# CMB temperature trispectrum from accreting primordial black holes

Trey W. Jensen  Yacine Ali-Haïmoud  Center for Cosmology and Particle Physics
Department of Physics, New York University
New York, NY 10003, USA

###### Abstract

It is known that Primordial Black Holes (PBHs) can leave an imprint on Cosmic
Microwave Background (CMB) anisotropy power spectra, due to their accretion-
powered injection of energy into the recombining plasma. Here we study a
qualitatively new CMB observable sourced by accreting PBHs: the temperature
trispectrum or connected 4-point function. This non-Gaussian signature is due
to the strong spatial modulation of the PBH accretion luminosity, thus
ionization perturbations, by large-scale supersonic relative velocities
between PBHs and the accreted baryons. We first derive a factorizable
quadratic transfer function for free-electron fraction inhomogeneities induced
by accreting PBHs. We then compute the perturbation to the CMB temperature
anisotropy due to a general modification of recombination, and apply our
results to accreting PBHs. We calculate a new contribution to the temperature
power spectrum due to the spatial fluctuations of the ionization perturbation
induced by accreting PBHs, going beyond past studies which only accounted for
its homogeneous part. While these contributions are formally comparable, we
find the new part to be subdominant, due to the poor correlation of the
perturbed temperature field with the standard CMB anisotropy. For the first
time, we compute the temperature trispectrum due to accreting PBHs. This
trispectrum is weakly correlated with the local-type primordial non-
Gaussianity trispectrum, hence constraints on the latter do not lead to
competitive bounds on accreting PBHs. We also forecast Planck’s sensitivity to
the temperature trispectrum sourced by accreting PBHs. Excitingly, we find it
to be more sensitive to PBHs under $\sim 10^{3}M_{\odot}$ than current
temperature-only power spectrum constraints. This result motivates our future
work extending this study to temperature and polarization trispectra induced
by inhomogeneously-accreting PBHs.

##  I Introduction

Not only are primordial black holes (PBHs) a probe of very early-Universe
physics, but they could also be the culprit behind several cosmological and
astrophysical mysteries. For instance, even if they constitute only a small
fraction of cold dark matter (CDM), intermediate-mass PBHs (1-$10^{4}$
$M_{\odot}$) could be the seed for supermassive black holes (Bean and
Magueijo, 2002) or account for recent LIGO/Virgo gravitational wave
observations (Bird  _et al._ , 2016). Thus, even if the abundance of PBHs in
this mass range is heavily constrained (Carr  _et al._ , 2021), it proves
invaluable to inspect further.

The intermediate-mass range is where PBH accretion may leave a non-negligible
signature on the cosmic microwave background (CMB). The underlying physical
phenomena that lead to an indirect signal are the following. PBHs accrete
primordial plasma throughout cosmic time; some fraction of the in-falling
material is converted into radiation; this radiation propagates and deposits
energy into the background recombining plasma, heating and ionizing it;
finally, this change to the ionization history perturbs the last-scattering
surface, ultimately altering the observed CMB temperature and polarization
anisotropy. In fact, the strongest constraints on the abundance of PBHs in
this mass range come from this effect Carr  _et al._ (2021); however, previous
literature have only looked for a signal in 2-point CMB anisotropy statistics
(Ricotti  _et al._ , 2008; Ali-Haïmoud and Kamionkowski, 2017; Poulin  _et
al._ , 2017).

One avenue that has not been inspected is the non-Gaussianity that is induced
in the CMB by accreting PBHs. Although the PBH accretion rate and radiation
power are largely uncertain, they necessarily depend on the magnitude of the
local relative velocity between the accreted matter (baryons) and PBHs, which
behave as CDM on large scales (Ricotti  _et al._ , 2008; Ali-Haïmoud and
Kamionkowski, 2017; Poulin  _et al._ , 2017). This dependence implies a
spatial modulation of the luminosity of accreting PBHs and thus
inhomogeneities in their perturbation to recombination. It is known that
inhomogeneous recombination generates non-Gaussian signatures in CMB
anisotropies (Senatore  _et al._ , 2009; Khatri and Wandelt, 2009; Dvorkin
_et al._ , 2013). The goal of this paper is to quantify this _qualitatively
different_ CMB signature of accreting PBHs, for the first time.

The effect considered here is similar in spirit to that studied in Ref.
Dvorkin  _et al._ (2013) in the context of dark matter (DM) annihilation,
with, however, two major differences. First, since the PBH luminosity depends
on the relative velocity squared, the lowest-order non-Gaussian statistics
induced by accreting PBHs is the _trispectrum_ , or connected 4-point
function. This is to be contrasted with the bispectrum (3-point function)
sourced by energy injection from inhomogeneous annihilating DM Dvorkin  _et
al._ (2013). Second, in the case of annihilating DM, the inhomogeneity in
energy injection is of order the DM density fluctuation around recombination,
that is of order $\sim 10^{-3}$ on scales $k\sim 0.1$ Mpc-1. In contrast, the
PBH luminosity has _order-unity_ fluctuations on the same scales Jensen and
Ali-Haïmoud (2021), as it is strongly modulated by supersonic relative
velocities Tseliakhovich and Hirata (2010). This implies that the
inhomogeneities in the free-electron fraction sourced by accreting PBHs are
comparable to its mean enhancement, as we demonstrated explicitly in Ref.
Jensen and Ali-Haïmoud (2021), hereafter Paper I. We therefore expect the
amplitude of the non-Gaussian signature of accreting PBHs to be $\sim 10^{3}$
times larger than that of inhomogeneously annihilating DM, at equal amplitudes
of the 2-point function perturbation.

In Paper I, we found that, for a PBH abundance saturating CMB power-spectra
limits, the free-electron perturbation is of order $\delta_{e}\sim 10^{-3}$
around $z\sim 10^{3}$, both in mean and in root-mean-square (see Fig. 14 of
Paper I). This relatively large effect implies that the CMB trispectrum could
be significantly more sensitive to PBHs than CMB power spectra, as we now show
with two simple order-of-magnitude estimates. First, without any exotic energy
injection nor primordial non-Gaussianity, recombination is intrinsically
inhomogeneous, with perturbations $\delta_{e,\rm std}\sim 10^{-4}$ Novosyadlyj
(2006); Senatore  _et al._ (2009). This leads to non-Gaussianities with an
amplitude just below detectability threshold for Planck Senatore  _et al._
(2009); Huang and Vernizzi (2013). This suggests that an inhomogeneity of
order $\delta_{e}\sim 10^{-3}$ would lead to a non-Gaussian signal with a
signal-to-noise ratio (SNR) of order 10. Second, in the presence of a
perturbation $\delta_{e}$ to the free-electron fraction, the CMB temperature
anisotropy $\Theta=\Theta^{(0)}+\Theta^{(1)}$ is displaced from its standard
value $\Theta^{(0)}\sim\zeta$, where $\zeta\sim 10^{-4.5}$ is the primordial
curvature perturbation, by an amount $\Theta^{(1)}\sim\delta_{e}\zeta$. We
therefore expect the connected 4-point function to be of order
$\langle\Theta\Theta\Theta\Theta\rangle_{\rm
c}=\langle\Theta^{(0)}\Theta^{(0)}\Theta^{(0)}\Theta^{(1)}\rangle\sim
10^{-3}\langle\zeta^{2}\rangle^{2}$ for a PBH abundance saturating CMB power-
spectra limits. In comparison, primordial trispectra lead to 4-point functions
of order $\langle\Theta\Theta\Theta\Theta\rangle_{\rm c}\sim g_{\rm
NL}\langle\zeta^{2}\rangle^{3}\sim 10^{-9}g_{\rm
NL}\langle\zeta^{2}\rangle^{2}$. Planck’s upper limits on the amplitude of
local-type primordial non-Gaussianity is $|g_{\rm NL}|\lesssim 10^{5}$ Planck
Collaboration (2020a), implying that Planck is sensitive to a 4-point function
of order $\langle\Theta\Theta\Theta\Theta\rangle_{\rm c}\sim
10^{-4}\langle\zeta^{2}\rangle^{2}$. Here again, this estimates indicates that
PBHs saturating CMB power spectra limits could lead to a trispectrum
detectable with SNR $\sim 10$. Put differently, the trispectrum could be
sensitive to PBH abundances an order of magnitude below current CMB power-
spectra limits. As an ancillary effect, the perturbation of CMB power spectra
induced by accreting PBHs ought to be modified by order unity when properly
accounting for the inhomogeneities in $\delta_{e}$, which were neglected in
past works Ricotti  _et al._ (2008); Ali-Haïmoud and Kamionkowski (2017);
Poulin  _et al._ (2017).

These promising estimates warrant a detailed calculation of the effects of
inhomogeneously-accreting PBHs on CMB power spectra and trispectra. In this
work, we take the first step in this program by computing the temperature-only
2-point and 4-point functions. We moreover forecast Planck’s sensitivity to
PBHs from the temperature trispectrum. We find that the inhomogeneity in
recombination only leads to a $\lesssim 10\%$ correction to the effect of
accreting PBHs on the temperature power spectrum. We also find that the
temperature trispectrum is approximately as sensitive to accreting PBHs as the
temperature power spectrum is, and is thus not quite as powerful a probe as
our simple order-of-magnitude estimates indicated. This is likely due to the
imperfect correlation between the standard temperature anisotropy
$\Theta^{(0)}$ and the perturbation $\Theta^{(1)}$ sourced by inhomogeneous
ionization fluctuations. Still, we find that, for $M_{\rm pbh}\lesssim
10^{3}M_{\odot}$, the CMB temperature trispectrum would be a more sensitive
probe of accreting PBHs than the temperature power spectrum is. This result
motivates exploring the full temperature and polarization trispectrum, which
we take up in future work.

The remainder of this paper is organized as follows. In Section II we begin by
briefly reviewing accreting PBHs as a source of inhomogeneous recombination.
By assuming spherical accretion and taking the luminosity prescription from
Ref. Ali-Haïmoud and Kamionkowski (2017) (hereafter AK17), we derive a
quadratic transfer function for the perturbed free-electron fraction. This
transfer function incorporates the radiation transport simulation and
perturbed recombination calculation from Paper I. We are able to make the
transfer function factorizable with some justified approximations specific to
accreting PBHs, which tremendously reduces the computational cost of
calculating the high-dimensional trispectrum.

In Section III we derive general equations for the perturbed temperature
anisotropy at first order in free-electron fraction perturbations, starting
from the Boltzmann-Einstein system, and using the line-of-sight method (Seljak
and Zaldarriaga, 1996). The results of this section are general and not
limited to perturbations from accreting PBHs. As in previous works (Dvorkin
_et al._ , 2013; Khatri and Wandelt, 2010) we neglect “feedback” terms in the
first-order perturbation. However, for the first time we quantify the error
induced by this approximation in the case of the power-spectrum perturbation
induced by homogeneous free-electron perturbations.

In Section IV we apply these results to recombination perturbations due to
accreting PBHs. We compute the perturbation to the temperature anisotropy
power spectrum sourced by the _inhomogeneous part_ of free-electron fraction
perturbations, which we find to be more than an order of magnitude smaller
than its counterpart induced by the homogeneous effect on the ionization
history. We moreover compute the temperature trispectrum induced by accreting
PBHs, given in Eq. (94), which is one of the main results of this work.

In Section V we extract new limits on PBH abundance from Planck upper bounds
on the local-shape primordial trispectrum Planck Collaboration (2020a), which
indirectly constrains the PBH-induced trispectrum with which it partially
overlaps. But due to a poor correlation between the two trispectra, the
constraints are an order of magnitude weaker than the constraints from the
power spectra analysis. We also forecast the sensitivity of Planck to the
temperature 4-point function induced by accreting PBHs, based on the optimal
trispectrum estimator of Ref. Smith  _et al._ (2015). We are able to make
these computations efficiently by pre-computing purely geometric rotational-
invariant coefficients. We find that the temperature trispectrum could probe
PBH abundances lower than current temperature-only power-spectrum limits for
$M_{\rm pbh}\lesssim 10^{3}M_{\odot}$. We conclude and discuss future work in
Section VI.

We discuss a few points in more detail in the Appendices. In Appendix A, we
justify the approximation of general non-linear functions of $v_{\rm bc}$ by a
biased tracer of $v_{\rm bc}^{2}$. We describe our numerical resolution and
convergence tests in Appendix B. We review a few useful properties of spin-
weighted spherical harmonics in Appendix C, which we then use in Appendix D to
derive simple expressions for the rotational-invariant quantities involved in
the trispectrum sensitivity forecast calculation. In Appendix E we compute the
auto-power spectrum of the temperature perturbation induced by accreting PBHs
and its correlation coefficient with the standard temperature anisotropy.
Latsly, in Appendix F, we inspect the redshift dependence of the signal-to-
noise ratio of the PBH-induced trispectrum.

##  II Perturbed recombination from accreting PBHs

In this section we briefly review the effect of accreting PBHs on the
ionization history. We derive an approximate factorized form for the free-
electron fraction fluctuations, quadratic in the initial perturbations, which
will help simplify our trispectrum calculations later on.

###  II.1 Effect of accreting PBHs on the ionization history: general
expressions

If present in the early Universe, PBHs would accrete baryons which would power
some radiation—at minimum, the heated, compressed and eventually ionized
accreted gas would emit free-free radiation. The PBH luminosity $L$ is a
function of the baryon sound speed $c_{s}$ and of the magnitude of the local
relative velocity between baryons and dark matter $\bm{v}_{\rm bc}(\bm{r})$
(both evaluated far from the accretion region). The detailed dependence is
estimated in AK17, accounting for Compton heating and Compton drag, and in two
limiting regimes for the ionization structure of the accretion flow;
throughout this paper, and unless otherwise stated, we will assume the most
conservative “collisionally-ionized” limit. Following AK17, we approximate the
effect of relative velocities by adding them in quadrature to the baryon sound
speed $c_{s}$, i.e. approximating $L(c_{s};v_{\rm bc}\neq 0)\approx
L(\sqrt{c_{s}^{2}+v_{\rm bc}^{2}};0)$. While the baryon sound speed is very
nearly homogeneous near recombination, relative velocities have large-scale
fluctuations, with rms values of order five times the sound speed
Tseliakhovich and Hirata (2010); as a consequence, the PBH luminosity
$L(\bm{r})=\overline{L}(1+\delta_{L}(z,\bm{r}))$ is strongly inhomogeneous,
tracing the large-scale fluctuations of relative velocities.

Assuming, to simplify, that PBHs all have the same mass $M_{\rm pbh}$ and make
a fraction $f_{\rm pbh}$ of the dark matter, their accretion-powered
luminosity leads to a volumetric energy injection rate

$\displaystyle\dot{\rho}_{\rm inj}(z,\bm{r})$
$\displaystyle=\overline{\dot{\rho}}_{\rm
inj}(z)\left(1+\delta_{L}(z,\bm{r})\right),$
$\displaystyle\overline{\dot{\rho}}_{\rm inj}(z)$ $\displaystyle\equiv f_{\rm
pbh}\frac{\overline{\rho}_{c}(z)}{M_{\rm pbh}}\overline{L}(z)$ (1)

where $z$ is the redshift and $\overline{\rho}_{c}$ is the mean dark matter
mass density. Note that this equation is trivially generalizable to an
extended mass distribution. This inhomogeneously-injected energy is partially
deposited at some later time, and some distance away from the injection site.
Some of this energy is deposited in the form of extra ionizations, leading to
a perturbation $\Delta x_{e}(z,\bm{r})$ to the free-electron fraction. The
latter is a convolution of the volumetric energy injection rate with a
dimensionless injection-to-ionization Green’s function. In Fourier space, this
convolution is a simple product:

$\Delta
x_{e}(z,\bm{k})=\int_{z}^{\infty}\frac{dz^{\prime}}{1+z^{\prime}}G_{x_{e}}^{\rm
inj}(z,z^{\prime},k)\frac{\overline{\dot{\rho}}_{\rm inj}}{n_{\rm
H}HE_{I}}\Big{|}_{z^{\prime}}\delta_{L}(z^{\prime},\bm{k}),$ (2)

where $n_{\rm H}$ is the mean number density of hydrogen, $H$ is the Hubble
rate, and $E_{I}\equiv 13.6$ eV is hydrogen’s ionization energy. The
homogeneous part of the ionization-fraction perturbation is obtained from a
similar time integral, involving the homogeneous part of the Green’s function:

$\overline{\Delta
x_{e}}(z)=\int_{0}^{a}\frac{dz^{\prime}}{1+z^{\prime}}G_{x_{e}}^{\rm
inj}(z,z^{\prime},0)\frac{\overline{\dot{\rho}}_{\rm inj}}{n_{\rm
H}HE_{I}}\Big{|}_{z^{\prime}}.$ (3)

In Paper I, we computed the Green’s function $G_{x_{e}}^{\rm
inj}(z,z^{\prime},k)$ numerically, by convolving the injection-to-deposition
Green’s function obtained from a radiative transfer code with the deposition-
to-ionization Green’s functions computed with a modified HYREC-2 Lee and Ali-
Haïmoud (2020); Ali-Haïmoud and Hirata (2011, 2010).

###  II.2 Quadratic transfer function of ionization perturbations

The scale dependence of the luminosity perturbations $\delta_{L}$ is non-
trivial, as the PBH luminosity is a nonlinear function of $v_{\rm bc}^{2}$.
However, as we will see below, at lowest order $\Delta x_{e}$ affects CMB
anisotropy statistics only through cross-correlations with other fields. As we
demonstrate in Appendix A, to a good approximation these cross-correlations
can be obtained by approximating the full function by a biased tracer of
$v_{\rm bc}^{2}$, with the same first moment:

$\delta_{L}(z,\bm{r})\approx b(z)\left(\frac{v_{\rm bc}^{2}(z,\bm{r})}{\langle
v_{\rm bc}^{2}\rangle(z)}-1\right),\ \ \ b\equiv\frac{3}{2}\frac{\langle
v_{\rm bc}^{2}\delta_{L}\rangle}{\langle v_{\rm bc}^{2}\rangle}.$ (4)

This approximation is most accurate in both the large-scale and small-scale
regimes, and as a consequence is reasonably accurate at all scales. We show
the bias parameter $b$ as a function of redshift in Fig. 1 for several black
hole masses, for the AK17 accretion luminosity model. It is systematically
negative, reflecting the suppression of accretion rate and luminosity in
regions of large relative velocity, and its absolute value is roughly of order
unity across a broad range of masses and redshifts. Although the accretion
model is highly uncertain, we expect that these qualitative features should be
robust, and hold even for very different accretion models, such as disk-like
accretion Poulin  _et al._ (2017).

Figure 1:  Bias parameter $b(z)$ of the PBH accretion luminosity, approximated
as a biased tracer of $v_{\rm bc}^{2}$ (see precise definition in Eq. 4). The
bias is shown as a function of redshift and for several PBH masses
$1M_{\odot}\leq M_{\rm pbh}\leq 10^{4}M_{\odot}$.

Assuming scalar initial conditions and linear evolution, the relative velocity
field is purely longitudinal, and we denote its transfer function by
$\widetilde{v}_{\rm bc}(z,k)$ defined such that

$\bm{v}_{\rm bc}(z,\bm{k})=-i\hat{k}~{}\widetilde{v}_{\rm
bc}(z,k)\zeta(\bm{k}),$ (5)

where $\zeta(\bm{k})$ is the primordial curvature perturbation. We then have

$\displaystyle v_{bc}^{2}(z,\bm{k})=$
$\displaystyle(\bm{v}_{bc}\cdot\bm{v}_{bc})(z,\bm{k})$ $\displaystyle=$
$\displaystyle-\int
D(k_{1}k_{2})\cancel{\delta}(\bm{k}_{1}+\bm{k}_{2}-\bm{k})(\hat{k}_{1}\cdot\hat{k}_{2})$
$\displaystyle\quad\quad\times\widetilde{v}_{\rm
bc}(z,k_{1})\widetilde{v}_{\rm
bc}(z,k_{2})\zeta(\bm{k}_{1})\zeta(\bm{k}_{2}),$ (6)

where from here on we denote $D(k_{1}\cdots k_{N})\equiv
d^{3}k_{1}/(2\pi)^{3}\cdots~{}d^{3}k_{N}/(2\pi)^{3}$ and
$\cancel{\delta}(\bm{k})\equiv(2\pi)^{3}\delta_{\rm D}(\bm{k})$.

We denote by $\delta_{e}\equiv\Delta
x_{e}/x_{e}^{(0)}=\overline{\delta}_{e}+\delta_{e,\rm inh}$ the total
fractional perturbation to the standard (and homogeneous) ionization history
$x_{e}^{(0)}$. The first part, $\overline{\delta}_{e}$, is the homogeneous
contribution, and the second part, $\delta_{e,\rm inh}$, is the inhomogeneity,
which has zero mean, $\langle\delta_{e,\rm inh}\rangle=0$.

Inserting Eq. (4) into Eq. (2), we obtain the Fourier transform of
$\delta_{e,\rm inh}$ for $\bm{k}\neq 0$:

$\displaystyle\delta_{e,\rm inh}(z,\bm{k})\equiv\frac{\Delta
x_{e}(z,\bm{k})}{x_{e}^{(0)}(z)}$ $\displaystyle\approx f_{\rm
pbh}\int\\!\\!D(k_{1}k_{2})\cancel{\delta}(\bm{k}_{1}+\bm{k}_{2}-\bm{k})$
$\displaystyle\quad\times
T_{e}(z,\bm{k}_{1},\bm{k}_{2})\zeta(\bm{k}_{1})\zeta(\bm{k}_{2}),$ (7)

where, for $\bm{k}_{1}+\bm{k}_{2}\neq\bm{0}$, the ionization-perturbation
quadratic transfer function $T_{e}$ is defined as

$\displaystyle T_{e}(z,\bm{k}_{1},\bm{k}_{2})$
$\displaystyle\equiv-\frac{\hat{k}_{1}\cdot\hat{k}_{2}}{x_{e}^{(0)}(z)}\int_{z}^{\infty}\frac{dz^{\prime}}{1+z^{\prime}}G_{x_{e}}^{\rm
inj}(z,z^{\prime},|\bm{k}_{1}+\bm{k}_{2}|)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\times\frac{\overline{\rho}_{c}\overline{L}~{}b}{M_{\rm
pbh}n_{\rm H}HE_{I}}\Big{|}_{z^{\prime}}\frac{\widetilde{v}_{\rm
bc}(k_{1})\widetilde{v}_{\rm bc}(k_{2})}{\langle v_{\rm
bc}^{2}\rangle}\Big{|}_{z^{\prime}}.$ (8)

We moreover define

$\displaystyle T_{e}(z,\bm{k}_{1},-\bm{k}_{1})=0,$ (9)

so that we may use Eq. (7) even for $\bm{k}=\bm{0}$, in which case it gives
$\delta_{e,\rm inh}(\bm{k}=\bm{0})=0$, as it should since $\delta_{e,\rm inh}$
is defined to have a vanishing spatial average111A more rigorous approach
would be to keep track of the term proportional to $\cancel{\delta}(\bm{k})$
in $\delta_{e,\rm inh}(\bm{k})$; upon cross-correlating with other fields,
this approach would give the same results as using Eq. (7) for all $\bm{k}$
and imposing $T_{e}(\bm{k}_{1},-\bm{k}_{1})=0$..

###  II.3 Factorized approximation of the quadratic ionization transfer
function

We now derive an approximate, factorized form for $T_{e}$, which will
tremendously simplify our subsequent calculations of CMB power spectra and
trispectra. We do so by making two approximations.

$(i)$ For $z\gtrsim 10^{3}$, the Green’s function $G_{x_{e}}^{\rm
inj}(z,z^{\prime})$ is peaked at $z^{\prime}\approx z$ (see Fig. 9 in Paper
I). We may therefore approximate the last ratio in Eq. (8) by its value at
$z^{\prime}=z$. For $z\lesssim 10^{3}$, the Green’s function is increasingly
broad; however, after kinematic decoupling at $z_{\rm dec}\approx 1020$,
relative velocities redshift as $\widetilde{v}_{\rm bc}(z,k)\propto(1+z)$,
independently of scale Tseliakhovich and Hirata (2010). Therefore, the last
term in Eq. (8) is independent of redshift for $z^{\prime}\lesssim z_{\rm
dec}$. We therefore make the following approximation in Eq. (8), which we
expect to be accurate at all redshifts:

$\frac{\widetilde{v}_{\rm bc}(k_{1})\widetilde{v}_{\rm bc}(k_{2})}{\langle
v_{\rm bc}^{2}\rangle}\Big{|}_{z^{\prime}}\approx\frac{\widetilde{v}_{\rm
bc}(k_{1})\widetilde{v}_{\rm bc}(k_{2})}{\langle v_{\rm
bc}^{2}\rangle}\Big{|}_{z}.$ (10)

This approximation implies the following simplification:

$\displaystyle T_{e}(z,\bm{k}_{1},\bm{k}_{2})$
$\displaystyle=(\hat{k}_{1}\cdot\hat{k}_{2})G_{e}(z,|\bm{k}_{1}+\bm{k}_{2}|)$
$\displaystyle\quad\quad\quad\quad\times\frac{\widetilde{v}_{\rm
bc}(z,k_{1})\widetilde{v}_{\rm bc}(z,k_{2})}{\langle v_{\rm
bc}^{2}\rangle_{z}},$ (11) $\displaystyle G_{e}(z,k)\equiv$
$\displaystyle-\int_{z}^{\infty}\frac{dz^{\prime}}{1+z^{\prime}}\frac{G_{x_{e}}^{\rm
inj}(z,z^{\prime},k)}{x_{e}^{(0)}(z)}\frac{\overline{\rho}_{c}\overline{L}~{}b}{M_{\rm
pbh}n_{\rm H}HE_{I}}\Big{|}_{z^{\prime}}.$ (12)

Figure 2:  Normalized injection-integrated Green’s function defined in Eq.
(12), at various redshifts, for 100-$M_{\odot}$ PBHs. This function is
approximately Gaussian with a characteristic cutoff $k_{*}(z)$, beyond which
ionization inhomogeneities are suppressed due to finite propagation of
injected photons.

$(ii)$ As illustrated in Fig. 2, we find that $G_{e}(z,k)$ is an approximately
Gaussian function of wavenumber, with a characteristic cutoff at a redshift-
dependent scale $k_{*}(z)$:

$G_{e}(z,k)\approx G_{e}(z,0)e^{-k^{2}/k_{*}^{2}(z)}.$ (13)

Given that $(\bm{k}_{1}+\bm{k}_{2})^{2}\leq 2k_{1}^{2}+2k_{2}^{2}$, we may
therefore approximately bracket $G_{e}$ as follows:

$\frac{G_{e}(z,\sqrt{2}k_{1})~{}G_{e}(z,\sqrt{2}k_{2})}{G_{e}(z,0)}\leq
G_{e}(z,|\bm{k}_{1}+\bm{k}_{2}|)\leq G_{e}(z,0).$ (14)

By default, we will conservatively approximate $G_{e}$ by the lower bound of
this range. This approximation is accurate at large scales
$k_{1},k_{2}\lesssim k_{*}(z)$, at which propagation effects are not relevant
to energy deposition.

With these two approximations, the quadratic ionization transfer function
takes on the factorized form

$\displaystyle T_{e}(z,\bm{k}_{1},\bm{k}_{2})$
$\displaystyle\approx(\hat{k}_{1}\cdot\hat{k}_{2})\Delta_{e}(z,k_{1})\Delta_{e}(z,k_{2}),$
(15) $\displaystyle\Delta_{e}(z,k)$
$\displaystyle\equiv\frac{G_{e}(z,\sqrt{2}k)}{\sqrt{G_{e}(z,0)}}\frac{\widetilde{v}_{\rm
bc}(z,k)}{\langle v_{\rm bc}^{2}\rangle_{z}^{1/2}},$ (16)

where we recall that this expression holds for $\bm{k}_{1}+\bm{k}_{2}\neq 0$
only, and that $T_{e}(z,\bm{k}_{1},-\bm{k}_{1})=0$.

Our approximation for $G_{e}(z,|\bm{k}_{1}+\bm{k}_{2}|)$ can significantly
underestimate the true signal at small scales, in particular for
$\bm{k}_{1}\approx-\bm{k}_{2}$ or $k_{1}\ll k_{2}$. Moreover, it modifies the
geometric dependence of the signal. In order to estimate the error that this
approximation induces, we will also show our results in the spatially-on-the-
spot approximation $G_{e}(z,k)\approx G_{e}(z,0)$, which systematically over-
estimates the signal. In that case, the quadratic ionization transfer function
still takes the form (15), but with
$\Delta_{e}(z,k)=\sqrt{G_{e}(z,0)}\frac{\widetilde{v}_{\rm bc}(z,k)}{\langle
v_{\rm bc}^{2}\rangle_{z}^{1/2}}$. We show $\Delta_{e}(z,k)$ as a function of
wavenumber and redshift in Fig. 3, both for our default approximation, and in
the on-the-spot limit.

Figure 3:  $\Delta_{e}(z,k)$ as defined in Eq. (16) for PBHs of 100
$M_{\odot}$ plotted for various redshifts as a function of scale (left) and
for various scales as a function of redshift (right). We also show this
function if we ignore photon propagation and consider energy deposition as
spatially on-the-spot (dashed lines). These plots reveal the shape of the
free-electron perturbations induced by $v_{\rm bc}^{2}$ as well as the
amplitude suppression when considering nonlocal energy deposition from
accreting PBHs.

##  III Temperature anisotropy from perturbed recombination: general equations

We turn to computing the temperature anisotropy in the presence of a general
deviation from the standard free-electron fraction evolution, including
spatial variations. Because observed CMB anisotropies are consistent with the
standard $\Lambda$CDM prediction and canonical homogeneous recombination, this
deviation is necessarily small, allowing for a perturbative treatment.

###  III.1 Temperature Boltzmann equation

The evolution of the phase-space distribution of photons is governed by the
Boltzmann-Einstein differential system. CMB photons follow geodesics in an
expanding universe subject to Thomson scattering off free electrons. Provided
photons remain thermal, they are described entirely by their temperature
fluctuations $\Theta(\eta,{\bm{x}},\hat{n})$ and transverse, symmetric trace-
free $3\times 3$ polarization tensor $P_{ab}(\eta,\bm{x},\hat{n})$, where
$\hat{n}$ is the propagation direction and $\eta$ is the conformal time. In
the conformal Newtonian gauge, the Boltzmann-Einstein equation for the
temperature perturbation is (Ma and Bertschinger, 1995)

$\displaystyle\frac{\textrm{d}\Theta}{\textrm{d}\eta}\equiv\dot{\Theta}+\hat{n}\cdot\nabla\Theta+\hat{n}\cdot\nabla\psi-\dot{\phi}$
$\displaystyle=\dot{\tau}~{}\mathcal{C}[\Theta,P_{ab},{\bm{v}}_{b}],$ (17)

where overdots denote partial derivatives with respect to $\eta$. In this
equation, ${\bm{v}}_{b}$ is the baryon velocity, and $\dot{\tau}\equiv
an_{e}\sigma_{T}$ is the conformal scattering rate by free electrons with
number density $n_{e}$, where $a$ is the scale factor, and $\sigma_{\rm T}$ is
the Thomson cross section. $\mathcal{C}$ is a linear operator encapsulating
the geometry of Thomson scattering:

$\displaystyle\mathcal{C}[\Theta,P_{ab},{\bm{v}}_{b}](\hat{n})$
$\displaystyle\equiv\mathcal{L}[\Theta,P_{ab},{\bm{v}}_{b}](\hat{n})-\Theta(\hat{n}),$
(18) $\displaystyle\mathcal{L}[\Theta,P_{ab},{\bm{v}}_{b}](\hat{n})$
$\displaystyle\equiv\Theta_{0}+\hat{n}\cdot\bm{v}_{b}+\hat{n}_{a}\hat{n}_{b}~{}\Pi_{ab},$
(19)

where $\Theta_{0}$ is the photon temperature monopole,

$\displaystyle\Theta_{0}\equiv\int\frac{d^{2}\hat{n}}{4\pi}\Theta(\hat{n}),$
(20)

and the symmetric trace-free tensor $\Pi_{ab}$ is a linear combination of the
photon quadrupole moment and the angle-averaged polarization tensor:

$\displaystyle\Pi_{ab}\equiv\int\frac{d^{2}\hat{n}}{4\pi}\left[\frac{1}{4}(3\hat{n}_{a}\hat{n}_{b}-\delta_{ab})\Theta(\hat{n})+\frac{3}{2}P_{ab}(\hat{n})\right].$
(21)

To close the Boltzmann-Einstein differential system, the evolution equations
for baryons, cold dark matter, neutrinos, and photon polarization are needed
(Ma and Bertschinger, 1995), but we do not explicitly list them here.

###  III.2 Standard solution

We now briefly review the standard solution obtained with scalar initial
conditions and for linear evolution, and given the standard, homogeneous free-
electron fraction $x_{e}^{(0)}$. The notation and expressions derived here
will be useful in the following section dealing with the perturbation to CMB
anisotropies induced by modified recombination. We denote all standard
variables by a superscript $(0)$, e.g. $\Theta^{(0)},\psi^{(0)}$, etc… For
short, we also denote
$\mathcal{C}^{(0)}\equiv\mathcal{C}[\Theta^{(0)},P_{ab}^{(0)},\bm{v}_{b}^{(0)}]$,
and similarly for $\mathcal{L}^{(0)}$.

The Boltzmann equation (17) is most easily solved in terms of the variable
$\Theta_{\rm eff}\equiv\Theta+\psi$, and in Fourier space. For the standard
case, it takes the form

$\displaystyle\dot{\Theta}_{\rm eff}^{(0)}+i\bm{k}\cdot\hat{n}\Theta_{\rm
eff}^{(0)}+\dot{\tau}^{(0)}\Theta_{\rm eff}^{(0)}=\dot{\tau}^{(0)}S^{(0)},$
(22) $\displaystyle
S^{(0)}\equiv\mathcal{L}^{(0)}+\psi^{(0)}+\frac{1}{\dot{\tau}^{(0)}}(\dot{\psi}^{(0)}+\dot{\phi}^{(0)}),$
(23)

The solution at an arbitrary conformal time is given by the line-of-sight
solution (Seljak and Zaldarriaga, 1996),

$\displaystyle\Theta^{(0)}_{\rm eff}(\eta,\bm{k},\hat{n})$
$\displaystyle=\exp\left(\int_{\eta}^{\eta_{0}}d\eta^{\prime\prime}~{}\dot{\tau}^{(0)}(\eta^{\prime\prime})\right)$
$\displaystyle\times\int_{0}^{\eta}d\eta^{\prime}~{}g(\eta^{\prime})S^{(0)}(\eta^{\prime},{\bm{k}},\hat{n})e^{i\bm{k}\cdot\hat{n}(\eta^{\prime}-\eta)},$
(24)

where $\eta_{0}$ is the conformal time today, and $g(\eta)$ is the standard
visibility function,

$\displaystyle
g(\eta)\equiv\dot{\tau}^{(0)}(\eta)\exp\left(-\int_{\eta}^{\eta_{0}}d\eta^{\prime}~{}\dot{\tau}^{(0)}(\eta^{\prime})\right),$
(25)

In particular, the line-of-sight solution today, and at the spatial origin
${\bm{x}}=0$, is given by (Seljak and Zaldarriaga, 1996)

$\displaystyle\Theta^{(0)}_{\rm eff}(\eta_{0},{\bm{x}}=0,\hat{n})=\int
Dk\int_{0}^{\eta_{0}}\\!\\!d\eta\,g(\eta)$ $\displaystyle
S^{(0)}(\eta,\bm{k},\hat{n})~{}e^{-i\bm{k}\cdot\hat{n}\chi},$ (26)

where from here on we denote $\chi\equiv\eta_{0}-\eta$.

Under scalar adiabatic initial conditions, the baryon velocity is purely
longitudinal, i.e. in Fourier space,
$\bm{v}_{b}^{(0)}(\bm{k})=-i\bm{k}\theta_{b}^{(0)}/k^{2}$. Moreover, we have
$\hat{n}_{a}\hat{n}_{b}\Pi_{ab}^{(0)}(\bm{k})=-\Pi^{(0)}P_{2}(\hat{k}\cdot\hat{n})$,
where $\Pi^{(0)}\equiv(F_{\gamma 2}+G_{\gamma 0}+G_{\gamma 2})/8$ is a
combination of the photon temperature quadrupole and polarization monopole and
quadrupole moments (here we used the notation of Ref. Ma and Bertschinger
(1995)). We thus have

$\displaystyle
S^{(0)}(\bm{k},\hat{n})=\Theta_{0}^{(0)}+\psi^{(0)}+\frac{1}{\dot{\tau}^{(0)}}(\dot{\psi}^{(0)}+\dot{\phi}^{(0)})$
$\displaystyle-\frac{i}{k}(\hat{k}\cdot\hat{n})\theta_{b}^{(0)}-P_{2}(\hat{k}\cdot\hat{n})\Pi^{(0)},$
(27)

where the dependence of $\Theta_{0}^{(0)},\bm{v}_{b}^{(0)}$, etc… on $\bm{k}$
is implicit.

To simplify this expression, note that
$-i\hat{k}\cdot\hat{n}~{}e^{-i\bm{k}\cdot\hat{n}\chi}=\partial_{k\chi}e^{-i\bm{k}\cdot\hat{n}\chi}$,
where $\partial_{k\chi}\equiv\frac{1}{k}\frac{\partial}{\partial\chi}$. In Eq.
(III.2) we may then conveniently substitute $S^{(0)}$ by an angle-independent
differential operator acting on the geometric exponential term,

$\displaystyle S^{(0)}(\bm{k},\hat{n})\rightarrow
S_{\partial}^{(0)}(\bm{k})\equiv\Theta_{0}^{(0)}+\psi^{(0)}+\frac{1}{\dot{\tau}^{(0)}}(\dot{\psi}^{(0)}+\dot{\phi}^{(0)})$
$\displaystyle+\frac{\theta_{b}^{(0)}}{k}\partial_{k\chi}+\Pi^{(0)}\left(\frac{3}{2}\partial_{k\chi}^{2}+\frac{1}{2}\right).$
(28)

Using this substitution, and from the Rayleigh formula,

$\displaystyle e^{-i\bm{k}\cdot\hat{n}\chi}$
$\displaystyle=\sum_{\ell}(-i)^{\ell}(2\ell+1)j_{\ell}(k\chi)P_{\ell}(\hat{n}\cdot\hat{k})$
(29) $\displaystyle=4\pi\sum_{\ell m}(-i)^{\ell}j_{\ell}(k\chi)Y_{\ell
m}(\hat{n})Y_{\ell m}^{*}(\hat{k}),$ (30)

where $j_{\ell}$ are the spherical Bessel functions, we may directly read off
the harmonic multipoles (for $\ell>0$) of the standard temperature anisotropy
from Eq. (III.2):

$\displaystyle\Theta^{(0)}_{\ell m}$
$\displaystyle=4\pi(-i)^{\ell}\\!\int\\!Dk\,Y^{*}_{\ell
m}({\hat{k}})\int_{0}^{\eta_{0}}\\!\\!d\eta\,g(\eta)S^{(0)}_{\partial}(\bm{k},\eta)j_{\ell}(k\chi),$
(31)

where the operator $S^{(0)}_{\partial}$ now acts on the Bessel function.

Lastly, we denote by $\widetilde{S}^{(0)}_{\partial}(k,\eta)$ the transfer
function of $S^{(0)}_{\partial}$, defined such that

$\displaystyle
S^{(0)}_{\partial}(\bm{k},\eta)=\widetilde{S}^{(0)}_{\partial}(k,\eta)\zeta(\bm{k}),$
(32)

where $\zeta(\bm{k})$ is the primordial curvature perturbation. We thus obtain

$\displaystyle\Theta^{(0)}_{\ell m}$ $\displaystyle=4\pi(-i)^{\ell}\int
Dk~{}Y^{*}_{\ell m}({\hat{k}})\Delta_{\ell}(k)\zeta(\bm{k}),$ (33)
$\displaystyle\Delta_{\ell}(k)$
$\displaystyle\equiv\int_{0}^{\eta_{0}}d\eta~{}g(\eta)~{}\widetilde{S}^{(0)}_{\partial}(k,\eta)j_{\ell}(k\chi).$
(34)

Assuming the primordial curvature perturbation is Gaussian, with power
spectrum $P_{\zeta}(k)$, i.e. such that

$\displaystyle\langle\zeta(\bm{k})\zeta(\bm{k}^{\prime})\rangle=(2\pi)^{3}\delta_{\rm
D}(\bm{k}+\bm{k}^{\prime})P_{\zeta}(k)\equiv\cancel{\delta}(\bm{k}+\bm{k}^{\prime})P_{\zeta}(k),$
(35)

the canonical temperature anisotropy angular power spectrum,
$\braket{\Theta^{(0)}_{\ell
m}\Theta^{*(0)}_{\ell^{\prime}m^{\prime}}}\equiv\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}C_{\ell}^{(0)}$,
is then given by

$\displaystyle C_{\ell}^{(0)}=4\pi\int
Dk~{}[\Delta_{\ell}(k)]^{2}P_{\zeta}(k).$ (36)

The $\Delta_{\ell}(k)$ are the temperature fluctuation multipole transfer
functions that can be extracted from cosmological codes such as CLASS (Blas
_et al._ , 2011). In practice, since we will need to compute similar integrals
later on, we compute the conformal time integral in Eq. (34) ourselves, using
only the source term transfer functions in Eq. (III.2) from CLASS. We also
compute the $k$-integral in Eq. (36) ourselves, and checked that our results
match those of CLASS to high accuracy. We discuss our numerical resolution and
convergence tests in Appendix B.

###  III.3 Temperature anisotropy due to perturbed recombination

We now suppose the free-electron fraction is perturbed,
$x_{e}=x_{e}^{(0)}(1+\delta_{e})$. Importantly, we make no assumption about
the spatial dependence of $\delta_{e}$, which in general has both a
homogeneous and an inhomogeneous piece. As a result of the modified Thomson
scattering rate
$\dot{\tau}=\dot{\tau}^{(0)}(1+\delta_{e})\equiv\dot{\tau}^{(0)}+\dot{\tau}^{(1)}$,
all matter and metric fields also get altered:
$\Theta=\Theta^{(0)}+\Theta^{(1)}$, $\psi=\psi^{(0)}+\psi^{(1)}$, etc… For
short, we again denote
$\mathcal{C}^{(1)}\equiv\mathcal{C}[\Theta^{(1)},P_{ab}^{(1)},\bm{v}_{b}^{(1)}]$,
and similarly for $\mathcal{L}^{(1)}$.

In general, matter and metric fields depend nonlinearly on $\delta_{e}$;
however, in the limit of small $\delta_{e}$, we may solve them with a
perturbative expansion in $\delta_{e}\ll 1$. The zero-th order equation is the
canonical Boltzmann-Einstein system discussed in Sec. III.2. At first order in
$\delta_{e}\ll 1$, the photon temperature Boltzmann equation is

$\displaystyle\dot{\Theta}^{(1)}+\hat{n}\cdot\nabla\Theta^{(1)}+\hat{n}\cdot\nabla$
$\displaystyle\psi^{(1)}-\dot{\phi}^{(1)}=\dot{\tau}^{(0)}\mathcal{C}^{(1)}+\dot{\tau}^{(1)}\mathcal{C}^{(0)}.$
(37)

It is convenient to rewrite this equations in terms of the variable
$\Theta_{\rm eff}^{(1)}\equiv\Theta^{(1)}+\psi^{(1)}$, as follows:

$\displaystyle\dot{\Theta}_{\rm eff}^{(1)}+\hat{n}\cdot\nabla\Theta_{\rm
eff}^{(1)}+\dot{\tau}^{(0)}\Theta_{\rm eff}^{(1)}$
$\displaystyle=\dot{\tau}^{(0)}S^{(1)},$ (38)

where the source term $S^{(1)}$ will be discussed shortly. Again, Eq. (38) can
be easily solved in Fourier space, with the familiar line-of-sight solution.
In particular, the order-one photon temperature perturbation at present time
$\eta_{0}$, and at the spatial origin, takes the form

$\displaystyle\Theta^{(1)}_{\rm eff}$
$\displaystyle(\eta_{0},{\bm{x}}=0,\hat{n})=$ $\displaystyle\int
Dk\int_{0}^{\eta_{0}}\\!\\!d\eta\,g(\eta)S^{(1)}(\bm{k},\eta,\hat{n})e^{-i{\bm{k}}\cdot\hat{n}\chi}.$
(39)

The first-order source term $S^{(1)}$ contains two pieces:

$\displaystyle S^{(1)}$ $\displaystyle=S^{(1)\rm d}+S^{(1)\rm f},$ (40)
$\displaystyle S^{(1)\rm d}$
$\displaystyle\equiv\delta_{e}*\mathcal{C}^{(0)},$ (41) $\displaystyle
S^{(1)\rm f}$
$\displaystyle\equiv\mathcal{L}^{(1)}+\psi^{(1)}+\frac{1}{\dot{\tau}^{(0)}}\left(\dot{\psi}^{(1)}+\dot{\phi}^{(1)}\right).$
(42)

The first piece $S^{(1)\rm d}$, we coin as the “direct” term, as it depends
directly on the perturbed free-electron fraction $\delta_{e}$, and otherwise
on zero-th order terms through $\mathcal{C}^{(0)}$, which can thus be
extracted in a relatively straightforward fashion from CLASS. Note that
$\delta_{e}*\mathcal{C}^{(0)}$ denotes a multiplication in real space, or a
convolution in Fourier space. The second piece $S^{(1)\rm f}$, we dub the
“feedback” term, as it depends on first-order terms; it thus requires solving
explicitly for the infinite Boltzmann hierarchy similar to that solved at
zeroth order, but with an additional source term, containing wavemode mixing
due to convolutions in Fourier space Khatri and Wandelt (2010).

As in previous studies Khatri and Wandelt (2009); Dvorkin  _et al._ (2013), we
will not solve for the feedback term in this work. However, we now quantify
its magnitude for the first time, in the limit of _homogeneous_ perturbations
to recombination.

###  III.4 Magnitude of the feedback term for homogeneous $\delta_{e}$

We consider the limiting case where
$\delta_{e}(\eta,\bm{x})=\overline{\delta}_{e}(\eta)$ is homogeneous. Our
perturbative expansion in $\delta_{e}$ applies just as well in this case, as
long as $\overline{\delta}_{e}\ll 1$. We shall only include the “direct”
source term, and then explicitly check our results against the exact output
from CLASS, which can handle arbitrary homogeneous perturbations to the
recombination history, thus effectively account for both “direct” and
“feedback” sources (although the calculation is not split this way in CLASS).

Let us rewrite the direct source term as

$\displaystyle S^{(1)\rm d}_{\rm
hom}=\overline{\delta}_{e}\mathcal{C}^{(0)}=\overline{\delta}_{e}\left(\mathcal{L}^{(0)}+\psi^{(0)}\right)-\overline{\delta}_{e}\Theta_{\rm
eff}^{(0)},$ (43)

where the subscript “hom” is there to remind the reader that we are
considering a homogeneous free-electron fraction in this section.

The contribution of the second term to the innermost integral of Eq. (III.3)
can be rewritten in the form

$\displaystyle\int_{0}^{\eta_{0}}d\eta~{}g(\eta)~{}\overline{\delta}_{e}(\eta)\Theta_{\rm
eff}^{(0)}(\eta,\bm{k},\hat{n})e^{-i\bm{k}\cdot\hat{n}\chi}$
$\displaystyle=\int_{0}^{\eta_{0}}d\eta~{}g(\eta)~{}\overline{\mathcal{D}}_{e}(\eta)S^{(0)}(\eta,\bm{k},\hat{n})e^{-i\bm{k}\cdot\hat{n}\chi},$
(44)

where

$\displaystyle\overline{\mathcal{D}}_{e}(\eta)\equiv\int_{\eta}^{\eta_{0}}d\eta^{\prime}~{}\dot{\tau}^{(0)}(\eta^{\prime})\overline{\delta}_{e}(\eta^{\prime}).$
(45)

To obtain this result, we inserted the arbitrary-time line-of-sight solution
(III.2) for $\Theta_{\rm eff}^{(0)}$, and switched the order of integration.
We therefore arrive at the following expression for the direct contribution to
the first-order temperature perturbation in the homogeneous case:

$\displaystyle{\Theta}^{(1)\rm d}_{\rm hom}(\eta_{0},{\bm{x}}=0,\hat{n})=\int
Dk\int_{0}^{\eta_{0}}\\!\\!d\eta\,g(\eta)$
$\displaystyle\left[\overline{\delta}_{e}(\eta)(\mathcal{L}^{(0)}_{\partial}+\psi^{(0)})-\overline{\mathcal{D}}_{e}(\eta)S^{(0)}_{\partial}\right]e^{-i{\bm{k}}\cdot\hat{n}\chi},$
(46)

where $\mathcal{L}^{(0)}_{\partial}(\bm{k},\eta)$ is the operator obtained
from $\mathcal{L}^{(0)}(\bm{k},\eta,\hat{n})$ in the same fashion as
$S^{(0)}_{\partial}$ is obtained from $S^{(0)}$ (c.f. Eq. (III.2)).

Using the same steps as in Sec. III.2, we thus arrive at the following
expression for the spherical-harmonic components of the direct-only part of
$\Theta^{(1)}_{\rm hom}$:

$\displaystyle{\Theta}_{\ell m,\rm hom}^{(1)\rm d}$
$\displaystyle=4\pi(-i)^{\ell}\int Dk~{}Y_{\ell
m}^{*}(\hat{k})\Delta_{\ell,\rm hom}^{(1)\rm d}(k)\zeta(\bm{k}),$ (47)
$\displaystyle\Delta_{\ell,\rm hom}^{(1)\rm d}(k)$
$\displaystyle\equiv\int_{0}^{\eta_{0}}d\eta~{}g(\eta)~{}\Big{[}\overline{\delta}_{e}(\eta)(\widetilde{\mathcal{L}}^{(0)}_{\partial}+\widetilde{\psi}^{(0)})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad-\overline{\mathcal{D}}_{e}(\eta)\widetilde{S}^{(0)}_{\partial}\Big{]}j_{\ell}(k\chi),$
(48)

where $\widetilde{\mathcal{L}}^{(0)}_{\partial}(k,\eta)$ and
$\widetilde{\psi}^{(0)}(k,\eta)$ are the transfer functions of
$\mathcal{L}^{(0)}_{\partial}(\bm{k},\eta)$ and $\psi^{(0)}(\bm{k},\eta)$.

We may now compute the perturbation to the angular power spectrum. To linear
order in $\overline{\delta}_{e}\ll 1$, we have
$C_{\ell}=C_{\ell}^{(0)}+{C}_{\ell,\rm hom}^{(1)}$, where we defined
$2\braket{{\Theta}^{(1)}_{\ell m,\rm
hom}\Theta^{*(0)}_{\ell^{\prime}m^{\prime}}}\equiv\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}{{C}^{(1)}_{\ell,\rm
hom}}$. We find that the direct contribution to ${C}_{\ell,\rm hom}^{(1)}$ is
then

$\displaystyle{C}_{\ell,\rm hom}^{(1)\rm d}=8\pi\int
Dk~{}P_{\zeta}(k)\Delta_{\ell}(k)\Delta_{\ell,\rm hom}^{(1)\rm d}(k).$ (49)

We computed ${C}_{\ell,\rm hom}^{(1)\rm d}$ using the homogeneous part of the
free-electron perturbation sourced by accreting PBHs, as calculated in AK17.
We compare this result against the exact ${C}_{\ell,\rm hom}^{(1)}$ obtained
from CLASS in Fig. 4. We see that neglecting the feedback term $S^{(1)\rm
f}_{\rm hom}$ leads to an order $\sim 10\%$ relative error on ${C}_{\ell,\rm
hom}^{(1)}$ for relevant black hole masses, indicating that the term is
subdominant. While there is no guarantee that this subdominance carries over
in general at higher-order statistics, it still gives us some confidence that
neglecting $S^{(1)\rm f}$ is a reasonable approximation, at least as a first
step, and especially considering the large theoretical uncertainty in the PBH
accretion model.

In what follows, we will therefore approximate $S^{(1)}\approx S^{(1)\rm
d}=\delta_{e}*\mathcal{C}^{(0)}$, and no longer indicate that we use the
direct-term only by a label “d”.

Figure 4:  Fractional change to the temperature anisotropy power spectrum from
the homogeneous perturbation to the free-electron fraction,
$\overline{\delta_{e}}(\eta)$, for various PBH masses and abundances. We
compare the exact non-perturbative effect extracted from CLASS to the
perturbative solution including only the “direct” source term discussed in
Sec. III.4. Our approximation of neglecting the “feedback” term is reasonably
accurate and we assume this carries over for the inhomogeneous free-electron
fraction case.

###  III.5 Alternative calculation of $C_{\ell,\rm hom}^{(1)\rm d}$

Before moving to the full calculation of $\Theta^{(1)}$, including ionization
fraction inhomogeneities, we present an alternative calculation of
$\Delta_{\ell,\rm hom}^{(1)\rm d}(k)$, required for the cross-power spectrum
$C_{\ell,\rm hom}^{(1)}$. This approach relies on intermediate quantities also
used for the trispectrum calculation, and provides a useful cross check of our
numerical methods.

For any quantity $X(k,\hat{n}\cdot\hat{k})$, we define its Legendre multipole
moments $X_{\ell}(k)$ as usual through

$\displaystyle
X_{\ell}(k)\equiv\frac{i^{\ell}}{2}\int_{-1}^{1}d\mu~{}P_{\ell}(\mu)X(k,\mu),$
(50)

such that

$\displaystyle X(k,\hat{k}\cdot\hat{n})$
$\displaystyle=\sum_{\ell}(-i)^{\ell}(2\ell+1)X_{\ell}(k)P_{\ell}(\hat{n}\cdot\hat{k})$
(51) $\displaystyle=4\pi\sum_{\ell m}(-i)^{\ell}X_{\ell}(k)Y_{\ell
m}(\hat{k})Y_{\ell m}^{*}(\hat{n}).$ (52)

We denote by $\widetilde{\mathcal{C}}^{(0)}(\eta,k,\hat{k}\cdot\hat{n})$ the
transfer function of $\mathcal{C}^{(0)}(\eta,\bm{k},\hat{n})$ (i.e. such that
$\mathcal{C}^{(0)}=\widetilde{\mathcal{C}}^{(0)}\zeta$), and define
$\mathcal{J}(\eta,k,\hat{n}\cdot\hat{k})\equiv
e^{-i\chi\bm{k}\cdot\hat{n}}\widetilde{\mathcal{C}}^{(0)}(\eta,k,\hat{k}\cdot\hat{n})$.

Substituting
$S^{(1)}(\bm{k},\eta,\hat{n})e^{-i\bm{k}\cdot\hat{n}\chi}=\overline{\delta}_{e}\mathcal{J}(\eta,k,\hat{n}\cdot\hat{k})\zeta(\bm{k})$
and inserting the spherical-harmonic expansion of $\mathcal{J}$ into Eq.
(III.3), we then arrive again at Eq. (47), with

$\displaystyle\Delta_{\ell,\rm hom}^{(1)\rm
d}(k)\equiv\int_{0}^{\eta_{0}}d\eta~{}g(\eta)~{}\overline{\delta}_{e}(\eta)\mathcal{J}_{\ell}(\eta,k).$
(53)

Using the plane-wave expansion (30) and the Legendre expansion of the product
of two Legendre polynomials, we may relate the coefficients
$\mathcal{J}_{\ell}$ to the Legendre coefficients of
$\widetilde{\mathcal{C}}^{(0)}$ as follows:

$\mathcal{J}_{\ell}(\eta,k)=\frac{4\pi}{2\ell+1}\sum_{\ell_{1}\ell_{2}}i^{\ell-\ell_{1}-\ell_{2}}(g_{\ell_{1}\ell_{2}\ell})^{2}j_{\ell_{1}}(k\chi)\widetilde{\mathcal{C}}^{(0)}_{\ell_{2}}(\eta,k),$
(54)

where $g_{\ell_{1}\ell_{2}\ell}$ is proportional to a three-J symbol, and is
defined in Eq. (C). Since this coefficient is nonvanishing only if
$\ell_{1}+\ell_{2}+\ell$ is even, we may substitute
$i^{\ell-\ell_{1}-\ell_{2}}=(-1)^{(\ell-\ell_{1}-\ell_{2})/2}=(-1)^{(\ell+3\ell_{1}+3\ell_{2})/2}$.
The Legendre coefficients of the collision operator are given explicitly by

$\displaystyle\widetilde{\mathcal{C}}^{(0)}_{\ell}=$
$\displaystyle\frac{1}{3}\widetilde{v}_{b\gamma}^{(0)}\delta_{\ell
1}+\frac{1}{5}\widetilde{\Pi}^{(0)}\delta_{\ell
2}-\widetilde{\Theta}_{\ell}^{(0)}(1-\delta_{\ell 0}-\delta_{\ell 1}).$ (55)

The sums over $\ell_{1}$ and $\ell_{2}$ in Eq. (54) are formally infinite, and
must be truncated in practice. Since the higher $\ell_{2}$-modes from the
collision term are induced after the peak of the visibility function, we
choose to truncate the $\ell_{2}$ sum at some finite $\ell_{\rm cut}$. This
automatically renders the double sum finite, since for a given $\ell_{2}$,
$\ell_{1}$ is bounded by the triangle condition,
$|\ell-\ell_{2}|\leq\ell_{1}\leq\ell+\ell_{2}$.

We compute $\Delta_{\ell,\rm hom}^{(1)d}$ as given by Eq. (53) and use it to
obtain $C_{\ell,\rm hom}^{(1)d}$ from Eq. (49). We show the results in Fig. 5,
for various $\ell_{\rm cut}$, and compare them to the result obtained with the
line-of-sight commutation method described in Sec. III.4. We see that the
former converges to the latter as $\ell_{\rm cut}$ is increased, as it should,
giving us confidence in the robustness of our numerical methods and results.

Figure 5:  Fractional change to the temperature anisotropy power spectrum from
the homogeneous perturbation to the free-electron fraction,
$\overline{\delta_{e}}(\eta)$, due to 100 $M_{\odot}$ accreting PBHs, computed
with the “commutation” method, using Eq. (48), or with the “direct summation
method”, using Eq. (53), where $\mathcal{J}_{\ell}$ is given by the sum (54),
truncated at $\ell_{2}\leq\ell_{\rm cut}$. We see that the direct summation
result converges to the “commutation” result as $\ell_{\rm cut}$ is increased,
as it should.

##  IV Perturbed temperature anisotropy statistics due to inhomogeneously-
accreting PBHs

###  IV.1 Temperature anisotropy transfer functions

####  IV.1.1 Definitions

Neglecting lensing and other nonlinearities, the standard temperature
perturbation is linearly related to the primordial curvature perturbation,
through (c.f. Eq (31))

$\displaystyle\Theta_{\ell m}^{(0)}$ $\displaystyle=\int Dk~{}T^{(0)}_{\ell
m}(\bm{k})~{}\zeta(\bm{k}),$ (56) $\displaystyle T_{\ell m}^{(0)}(\bm{k})$
$\displaystyle\equiv 4\pi(-i)^{\ell}\Delta_{\ell}(k)Y_{\ell m}^{*}(\hat{k}).$
(57)

Approximating the free-electron fraction perturbation due to accreting PBHs as
quadratic in the initial conditions (c.f. Eq. (7)), the corresponding
temperature anisotropy perturbation due to accreting PBHs is _cubic_ in the
initial curvature perturbation. The goal of this section is to derive an
explicit expression for the cubic transfer function $T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$, defined through

$\displaystyle\Theta_{\ell m,\rm inh}^{(1)}=f_{\rm pbh}\int
D(k_{1}k_{2}k_{3})T_{\ell m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle\times\zeta(\bm{k}_{1})\zeta(\bm{k}_{2})\zeta(\bm{k}_{3}),$ (58)

where the label “inh” indicates that here we focus on the
inhomogeneous-$\delta_{e}$ contribution to $\Theta^{(1)}$, recalling that it
also has a piece $\Theta^{(1)}_{\rm hom}$ due to the homogeneous
$\overline{\delta}_{e}$, which we computed in Sec. III.4, so that the total
$\Theta^{(1)}=\Theta^{(1)}_{\rm hom}+\Theta^{(1)}_{\rm inh}$.

In addition, we shall derive the harmonic coefficients of this cubic transfer
function, defined as

$\displaystyle T_{\ell m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=(4\pi)^{3}\sum_{\ell_{1}\ell_{2}\ell_{3}}(-i)^{\ell_{1}+\ell_{2}+\ell_{3}}$
$\displaystyle\times\sum_{m_{1}m_{2}m_{3}}T_{\ell_{1}\ell_{2}\ell_{3};\ell}^{m_{1}m_{2}m_{3};m}(k_{1},k_{2},k_{3})$
$\displaystyle\times
Y_{\ell_{1}m_{1}}(\hat{k}_{1})Y_{\ell_{2}m_{2}}(\hat{k}_{2})Y_{\ell_{3}m_{3}}(\hat{k}_{3}).$
(59)

####  IV.1.2 Calculation

Neglecting the “feedback” term, the source term for the line-of-sight solution
of $\Theta^{(1)}_{\rm inh}$ is the convolution between the collision term and
inhomogeneous part of the free-electron fraction, $S^{(1)}=\delta_{e,\rm
inh}*\mathcal{C}^{(0)}$. Using Eqs. (7) and (15), we can write explicitly

$\displaystyle S^{(1)}(\eta,\bm{k},\hat{n})$ $\displaystyle=f_{\rm
pbh}\int~{}D(k_{1}k_{2}k_{3})\cancel{\delta}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}-\bm{k})$
$\displaystyle\quad\times(\hat{k}_{1}\cdot\hat{k}_{2})\Delta_{e}(\eta,k_{1})\Delta_{e}(\eta,k_{2})$
$\displaystyle\quad\times\widetilde{\mathcal{C}}^{(0)}(\eta,k_{3},\hat{k}_{3}\cdot\hat{n})\zeta(\bm{k}_{1})\zeta(\bm{k}_{2})\zeta(\bm{k}_{3}),$
(60)

where again $\widetilde{\mathcal{C}}^{(0)}(\eta,k,\hat{k}\cdot\hat{n})$ is the
transfer function of $\mathcal{C}^{(0)}(\eta,\bm{k},\hat{n})$. Taking the
harmonic transform of the line-of-sight solution for $\Theta^{(1)}$, Eq.
(III.3), we then find

$\displaystyle T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\int_{0}^{\eta_{0}}d\eta~{}g(\eta)\int
d^{2}\hat{n}~{}Y_{\ell m}^{*}(\hat{n})$
$\displaystyle\times(\hat{k}_{1}\cdot\hat{k}_{2})\Delta_{e}(\eta,k_{1})\Delta_{e}(\eta,k_{2})$
$\displaystyle\times\widetilde{\mathcal{C}}^{(0)}(\eta,k_{3},\hat{k}_{3}\cdot\hat{n})e^{-i\chi\hat{n}\cdot(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})}.$
(61)

Note that this function is symmetric under exchange of $\bm{k}_{1}$ and
$\bm{k}_{2}$. Let us recall, also, that the expression above only holds for
$\bm{k}_{1}+\bm{k}_{2}\neq 0$, and that $T_{\ell
m}^{(1)}(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{3})=0$, since we are only considering
the inhomogeneous part (with zero mean) of the free-electron fraction
perturbation.

To obtain the harmonic coefficients of $T_{\ell m}^{(1)}$, we first rewrite
(denoting $\bm{\chi}\equiv\chi\hat{n}$),

$\displaystyle(\bm{k}_{1}\cdot\bm{k}_{2})e^{-i\bm{\chi}\cdot(\bm{k}_{1}+\bm{k}_{2})}$
$\displaystyle=-\left[\bm{\nabla}_{\bm{\chi}}e^{-i\bm{\chi}\cdot\bm{k}_{1}}\right]\cdot\left[\bm{\nabla}_{\bm{\chi}}e^{-i\bm{\chi}\cdot\bm{k}_{2}}\right]$
$\displaystyle=-\partial_{\chi}\left(e^{-i\bm{\chi}\cdot\bm{k}_{1}}\right)\partial_{\chi}\left(e^{-i\bm{\chi}\cdot\bm{k}_{2}}\right)$
$\displaystyle\quad-\frac{1}{\chi^{2}}\left[\bm{\nabla}_{\hat{n}}e^{-i\bm{\chi}\cdot\bm{k}_{1}}\right]\cdot\left[\bm{\nabla}_{\hat{n}}e^{-i\bm{\chi}\cdot\bm{k}_{2}}\right],$
(62)

where $\bm{\nabla}_{\bm{\chi}}$ is the gradient with respect to $\bm{\chi}$,
which we have split into its radial part $\hat{n}\partial_{\chi}$ and its
angular part $\frac{1}{\chi}\bm{\nabla}_{\hat{n}}$. Using the plane-wave
expansion (30), we thus have

$\displaystyle(\hat{k}_{1}\cdot\hat{k}_{2})e^{-i\chi\hat{n}\cdot(\bm{k}_{1}+\bm{k}_{2})}=-(4\pi)^{2}\sum_{\ell_{1}\ell_{2}}(-i)^{\ell_{1}+\ell_{2}}$
$\displaystyle\sum_{m_{1}m_{2}}Y_{\ell_{1}m_{1}}(\hat{k}_{1})Y_{\ell_{2}m_{2}}(\hat{k}_{2})$
$\displaystyle\times\Big{[}j_{\ell_{1}}^{\prime}(\chi
k_{1})j_{\ell_{2}}^{\prime}(\chi
k_{2})Y_{\ell_{1}m_{1}}^{*}(\hat{n})Y_{\ell_{2}m_{2}}^{*}(\hat{n})$
$\displaystyle+\frac{j_{\ell_{1}}(\chi k_{1})}{\chi
k_{1}}\frac{j_{\ell_{1}}(\chi k_{2})}{\chi
k_{2}}\bm{\nabla}_{\hat{n}}Y_{\ell_{1}m_{1}}^{*}(\hat{n})\cdot\bm{\nabla}_{\hat{n}}Y_{\ell_{2}m_{2}}^{*}(\hat{n})\Big{]}.$
(63)

Combining this result with the Legendre-expansion of
$e^{-i\chi\hat{n}\cdot\bm{k}_{3}}\widetilde{\mathcal{C}}^{(0)}(\eta,k_{3},\hat{k}_{3}\cdot\hat{n})$
[Eq. (54)], we are now in the position to compute
$T_{\ell_{1}\ell_{2}\ell_{3};\ell}^{m_{1}m_{2}m_{3};m}$ defined in Eq.
(IV.1.1):

$\displaystyle
T_{\ell_{1}\ell_{2}\ell_{3};\ell}^{m_{1}m_{2}m_{3};m}(k_{1},k_{2},k_{3})$
$\displaystyle=A_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})Q_{\ell_{1}\ell_{2}\ell_{3}\ell}^{m_{1}m_{2}m_{3}m}$
$\displaystyle+B_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell}^{m_{1}m_{2},m_{3}m},$
(64)

where the rotationally-invariant coefficients $A_{\ell_{1}\ell_{2},\ell_{3}}$
and $B_{\ell_{1}\ell_{2},\ell_{3}}$ are given by

$\displaystyle A_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})\equiv$
$\displaystyle-\int\\!\\!d\eta~{}g(\eta)j_{\ell_{1}}^{\prime}(\chi
k_{1})\Delta_{e}(\eta,k_{1})$ $\displaystyle\quad\times
j_{\ell_{2}}^{\prime}(\chi
k_{2})\Delta_{e}(\eta,k_{2})\mathcal{J}_{\ell_{3}}(\eta,k_{3}),$ (65)
$\displaystyle B_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})\equiv$
$\displaystyle-\int\\!\\!d\eta~{}g(\eta)\frac{j_{\ell_{1}}(\chi k_{1})}{\chi
k_{1}}\Delta_{e}(\eta,k_{1})$ $\displaystyle\quad\times\frac{j_{\ell_{2}}(\chi
k_{2})}{\chi k_{2}}\Delta_{e}(\eta,k_{2})\mathcal{J}_{\ell_{3}}(\eta,k_{3}),$
(66)

and the purely geometric terms
$Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}$ and
$\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}$
are integrals of the product of four spherical harmonics or their gradients:

$\displaystyle
Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}\\!\equiv$
$\displaystyle\quad\quad\\!\int\\!d^{2}\hat{n}~{}Y^{*}_{\ell_{1}m_{1}}(\hat{n})Y^{*}_{\ell_{2}m_{2}}(\hat{n})Y^{*}_{\ell_{3}m_{3}}(\hat{n})Y_{\ell_{4}m_{4}}^{*}(\hat{n}),$
(67)
$\displaystyle\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}\equiv$
$\displaystyle\quad\quad\int
d^{2}\hat{n}~{}\bm{\nabla}_{\hat{n}}Y^{*}_{\ell_{1}m_{1}}(\hat{n})\cdot\bm{\nabla}_{\hat{n}}Y^{*}_{\ell_{2}m_{2}}(\hat{n})Y^{*}_{\ell_{3}m_{3}}(\hat{n})Y_{\ell_{4}m_{4}}^{*}(\hat{n}).$
(68)

Note that we have separated the groups of indices on which the functions
depend fully symmetrically: $A_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})$
and $B_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})$ are symmetric under
exchange of $(\ell_{1},k_{1})$ with $(\ell_{2},k_{2})$,
$Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}$ is symmetric
under exchange of any two $(\ell,m)$ pairs, and
$\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}$ is
symmetric under exchange of $(\ell_{1},m_{1})$ with $(\ell_{2},m_{2})$, as
well as under exchange of $(\ell_{3},m_{3})$ with $(\ell_{4},m_{4})$.

###  IV.2 Perturbed temperature angular power spectrum

We have derived all the required transfer functions and are now equipped to
compute statistical properties of $\Theta^{(1)}_{\rm inh}$. Because the
perturbed temperature anisotropy is cubic in the primordial curvature
perturbation, it has a non-vanishing cross-correlation with the standard
temperature anisotropy. Using Eqs. (31) and (58), we have

$\displaystyle\langle\Theta_{\ell m,\rm
inh}^{(1)}\Theta_{\ell^{\prime}m^{\prime}}^{*(0)}\rangle$
$\displaystyle=f_{\rm
pbh}\int\\!\\!D(k_{1}k_{2}k_{3}k^{\prime})~{}T^{(1)}_{\ell
m}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ $\displaystyle\quad\times
T^{*(0)}_{\ell^{\prime}m^{\prime}}(\bm{k}^{\prime})\langle\zeta(\bm{k}_{1})\zeta(\bm{k}_{2})\zeta(\bm{k}_{3})\zeta^{*}(\bm{k}^{\prime})\rangle.$
(69)

Using Wick’s theorem, and recalling that $T^{(1)}_{\ell
m}(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{3})=0$, and that $T^{(1)}_{\ell m}$ is
symmetric in its first two arguments, we get

$\displaystyle\langle\Theta_{\ell m,\rm
inh}^{(1)}\Theta_{\ell^{\prime}m^{\prime}}^{*(0)}\rangle=$ $\displaystyle
2f_{\rm pbh}\int D(kk^{\prime})T^{(1)}_{\ell
m}(\bm{k}^{\prime},-\bm{k},\bm{k})$ $\displaystyle\quad\times
T^{*(0)}_{\ell^{\prime}m^{\prime}}(\bm{k}^{\prime})P_{\zeta}(k)P(k^{\prime}),$
(70)

From Eq. (61), we have

$\displaystyle T^{(1)}_{\ell
m}(\bm{k}^{\prime},-\bm{k},\bm{k})=-(\hat{k}^{\prime}\cdot\hat{k})\int_{0}^{\eta_{0}}d\eta~{}g(\eta)\Delta_{e}(\eta,k^{\prime})\Delta_{e}(\eta,k)$
$\displaystyle\times\int d^{2}\hat{n}~{}Y_{\ell
m}^{*}(\hat{n})\widetilde{\mathcal{C}}^{(0)}(\eta,k,\hat{k}\cdot\hat{n})e^{-i\chi\hat{n}\cdot\bm{k}^{\prime}}.$
(71)

Averaging over the direction $\hat{k}$, we then obtain,

$\displaystyle\int\frac{d^{2}\hat{k}}{4\pi}T^{(1)}_{\ell
m}(\bm{k}^{\prime},-\bm{k},\bm{k})$
$\displaystyle=i\int_{0}^{\eta_{0}}d\eta~{}g(\eta)\Delta_{e}(\eta,k^{\prime})\Delta_{e}(\eta,k)\widetilde{\mathcal{C}}^{(0)}_{1}(\eta,k)$
$\displaystyle\times\int d^{2}\hat{n}~{}Y_{\ell
m}^{*}(\hat{n})(\hat{k}^{\prime}\cdot\hat{n})e^{-i\chi\hat{n}\cdot\bm{k}^{\prime}},$
(72)

where $\widetilde{\mathcal{C}}^{(0)}_{1}(\eta,k)\equiv
i\frac{1}{2}\int_{-1}^{1}d\mu
P_{1}(\mu)\widetilde{\mathcal{C}}^{(0)}(\eta,k,\mu)$ is the order-1 Legendre
coefficient of $\widetilde{\mathcal{C}}^{(0)}$, which is proportional to the
baryon-photon relative velocity (or baryon-photon slip):

$\displaystyle\widetilde{\mathcal{C}}^{(0)}_{1}(\eta,k)=\frac{1}{3}\widetilde{v}_{b\gamma}(\eta,k),$
(73)

where we defined
$\bm{v}_{b\gamma}(\bm{k})\equiv(\bm{v}_{b}-\bm{v}_{\gamma})(\bm{k})=-i\hat{k}\widetilde{v}_{b\gamma}(k)\zeta(\bm{k})$.

Using the plane-wave expansion Eq. (30), this expression further simplifies to

$\displaystyle\int\frac{d^{2}\hat{k}}{4\pi}T^{(1)}_{\ell
m}(\bm{k}^{\prime},-\bm{k},\bm{k})$
$\displaystyle=-\frac{4\pi}{3}(-i)^{\ell}\int_{0}^{\eta_{0}}d\eta~{}g(\eta)\Delta_{e}(\eta,k^{\prime})\Delta_{e}(\eta,k)$
$\displaystyle\times\widetilde{v}_{b\gamma}(\eta,k)j_{\ell}^{\prime}(k^{\prime}\chi)Y_{\ell
m}^{*}(\hat{k}^{\prime}).$ (74)

Finally inserting Eqs. (57) and (74) into Eq. (70), we arrive at the following
simple result

$\displaystyle\langle\Theta_{\ell m,\rm
inh}^{(1)}\Theta_{\ell^{\prime}m^{\prime}}^{*(0)}\rangle=\frac{1}{2}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}C_{\ell,\rm
inh}^{(1)},$ (75)

where the cross-power spectrum is given by the conformal time integral

$\displaystyle C_{\ell,\rm inh}^{(1)}$ $\displaystyle=-\frac{16\pi}{3}f_{\rm
pbh}\int_{0}^{\eta}d\eta~{}g(\eta)\gamma(\eta)\mu_{\ell}(\eta),$ (76)

where we have defined

$\displaystyle\gamma(\eta)$ $\displaystyle\equiv\int
Dk~{}P_{\zeta}(k)\Delta_{e}(\eta,k)\widetilde{v}_{b\gamma}(\eta,k),$ (77)
$\displaystyle\mu_{\ell}(\eta)$ $\displaystyle\equiv\int
Dk~{}P_{\zeta}(k)\Delta_{e}(\eta,k)\Delta_{\ell}(k)j_{\ell}^{\prime}(k\chi).$
(78)

We see that the factorization of the free-electron perturbation transfer
function has allowed us to obtain a very simple expression for $C_{\ell,\rm
inh}^{(1)}$: it only requires pre-computing tables of $\mu_{\ell}(\eta)$ and
$\gamma(\eta)$, and then computing a one-dimensional integral. It can
equivalently be rewritten in the same form as Eq. (49):

$\displaystyle C_{\ell,\rm inh}^{(1)}$ $\displaystyle=$ $\displaystyle
8\pi\int Dk~{}P_{\zeta}(k)\Delta_{\ell}(k)\Delta_{\ell,\rm inh}^{(1)}(k),$
(79) $\displaystyle\Delta_{\ell,\rm inh}^{(1)}(k)$ $\displaystyle\equiv$
$\displaystyle-\frac{2}{3}f_{\rm
pbh}\int_{0}^{\eta_{0}}d\eta~{}g(\eta)\gamma(\eta)\Delta_{e}(k,\eta)j_{\ell}^{\prime}(k\chi).~{}~{}~{}$
(80)

Equation (76), or the equivalent form (79), constitute one of the main results
of this work.

We computed $C_{\ell,\rm inh}^{(1)}$ from Eq. (76), and checked that Eq. (79)
gives the same result. We show the result in Fig. 6, where we compare this
term to its counterpart $C_{\ell,\rm hom}^{(1)}$ sourced by the homogeneous
part of the free-electron fraction perturbation, for 100-$M_{\odot}$ accreting
PBHs. Even though these two contributions should in principle be comparable,
given that $\langle\delta_{e}^{2}\rangle^{1/2}\sim\overline{\delta}_{e}$ (see
Paper I), we find that $C_{\ell,\rm inh}^{(1)}$ is suppressed by a factor
$\sim 10-100$, depending on scale, relative to $C_{\ell,\rm hom}^{(1)}$ for
all black hole masses. This turns out to be due to both a poor correlation
between $\Theta^{(1)}_{\rm inh}$ and $\Theta^{(0)}$, and a suppression of the
characteristic amplitude of $\Theta^{(1)}_{\rm inh}$ itself. We expound on
this point in Appendix E.

Figure 6:  Fractional change to the temperature anisotropy power spectrum due
to accreting PBHs of 100 $M_{\odot}$ comprising all the dark matter. _Left:_
Comparison of the contribution due to the inhomogeneous part of the ionization
fraction perturbations, $C_{\ell,\rm inh}^{(1)}$, calculated for the first
time in this work, with the one arising from the homogeneous part of the free-
electron fraction, $C_{\ell,\rm hom}^{(1)}$, previously computed in AK17; we
also overlay the total change to the temperature power spectrum from both.
_Right_ : the ratio between $C_{\ell,\rm inh}^{(1)}$ and $C_{\ell,\rm
hom}^{(1)}$. Although one would expect $C_{\ell,\rm inh}^{(1)}$ to be of the
same order of magnitude as $C_{\ell,\rm hom}^{(1)}$ _a priori_ , we find in
practice that the former is $\sim 10-100$ times smaller that the latter. In
both cases, dashed curves correspond to the spatial on-the-spot approximation,
which neglects the spatial smearing of energy deposition due to the finite
propagation of injected photons.

###  IV.3 Temperature trispectrum

We now compute the connected four-point correlation function of temperature
anisotropy,

$\displaystyle\langle\Theta_{1}\Theta_{2}\Theta_{3}\Theta_{4}\rangle_{c}$
$\displaystyle\equiv\langle\Theta_{1}\Theta_{2}\Theta_{3}\Theta_{4}\rangle-\langle\Theta_{1}\Theta_{2}\rangle\langle\Theta_{3}\Theta_{4}\rangle$
$\displaystyle-\langle\Theta_{1}\Theta_{3}\rangle\langle\Theta_{2}\Theta_{4}\rangle-\langle\Theta_{1}\Theta_{4}\rangle\langle\Theta_{2}\Theta_{3}\rangle,$
(81)

where the numbered subscripts index both $\ell$ and $m$,
$\Theta_{1}\equiv\Theta_{\ell_{1}m_{1}}$, and $c$ denotes subtracting out the
unconnected parts of the trispectrum. Recalling that
$\Theta=\Theta^{(0)}+\Theta^{(1)}_{\rm hom}+\Theta^{(1)}_{\rm inh}$, and that
$\Theta^{(0)}$ and $\Theta^{(1)}_{\rm hom}$ are both linear in the initial
Gaussian curvature perturbation, to lowest order in electron density
perturbations, the trispectrum is given by

$\displaystyle\langle\Theta_{1}\Theta_{2}\Theta_{3}\Theta_{4}\rangle_{c}$
$\displaystyle=\langle\Theta_{1,\rm
inh}^{(1)}\Theta_{2}^{(0)}\Theta_{3}^{(0)}\Theta_{4}^{(0)}\rangle_{c}$
$\displaystyle+\langle\Theta_{1}^{(0)}\Theta_{2,\rm
inh}^{(1)}\Theta_{3}^{(0)}\Theta_{4}^{(0)}\rangle_{c}$
$\displaystyle+\langle\Theta_{1}^{(0)}\Theta_{2}^{(0)}\Theta_{3,\rm
inh}^{(1)}\Theta_{4}^{(0)}\rangle_{c}$
$\displaystyle+\langle\Theta_{1}^{(0)}\Theta_{2}^{(0)}\Theta_{3}^{(0)}\Theta_{4,\rm
inh}^{(1)}\rangle_{c}.$ (82)

We may now compute each term using Eq. (58) for $\Theta_{\ell m,\rm
inh}^{(1)}$. For instance, the last term is

$\displaystyle\langle\Theta_{1}^{(0)}\Theta_{2}^{(0)}\Theta_{3}^{(0)}\Theta_{4,\rm
inh}^{(1)}\rangle_{c}=f_{\rm pbh}\int
D(kk^{\prime}k^{\prime\prime})T_{4}^{(1)}(\bm{k},\bm{k}^{\prime},\bm{k}^{\prime\prime})$
$\displaystyle\times\Big{[}\langle\zeta(\bm{k})\zeta(\bm{k}^{\prime})\zeta(\bm{k}^{\prime\prime})\Theta^{(0)}_{1}\Theta^{(0)}_{2}\Theta^{(0)}_{3}\rangle$
$\displaystyle-\langle\zeta(\bm{k})\zeta(\bm{k}^{\prime})\zeta(\bm{k}^{\prime\prime})\Theta^{(0)}_{1}\rangle\langle\Theta_{2}^{(0)}\Theta_{3}^{(0)}\rangle$
$\displaystyle-\langle\zeta(\bm{k})\zeta(\bm{k}^{\prime})\zeta(\bm{k}^{\prime\prime})\Theta^{(0)}_{2}\rangle\langle\Theta_{1}^{(0)}\Theta_{3}^{(0)}\rangle$
$\displaystyle-\langle\zeta(\bm{k})\zeta(\bm{k}^{\prime})\zeta(\bm{k}^{\prime\prime})\Theta^{(0)}_{3}\rangle\langle\Theta_{1}^{(0)}\Theta_{2}^{(0)}\rangle\Big{]},$
(83)

where we used the explicit definition of the connected 4-point function, Eq.
(81). Using Wick’s theorem to compute the 6-point and 4-point functions of
Gaussian fields appearing in the integrand above, simplifying, and renaming
dummy integration variables, we arrive at

$\displaystyle\langle\Theta_{1}^{(0)}\Theta_{2}^{(0)}\Theta_{3}^{(0)}\Theta_{4,\rm
inh}^{(1)}\rangle_{c}=f_{\rm pbh}\int D(kk^{\prime}k^{\prime\prime})$
$\displaystyle\times\langle\zeta(\bm{k})\Theta_{1}^{(0)}\rangle\langle\zeta(\bm{k}^{\prime})\Theta_{2}^{(0)}\rangle\langle\zeta(\bm{k}^{\prime\prime})\Theta_{3}^{(0)}\rangle$
$\displaystyle\times\left[T_{4}^{(1)}(\bm{k},\bm{k}^{\prime},\bm{k}^{\prime\prime})+5\textrm{
perms.}\right]$ (84)

where the 5 permutations involve all other possible permutations of
$\bm{k},\bm{k}^{\prime},\bm{k}^{\prime\prime}$. The relevant two-point
functions are easily computed with the line-of-sight expression for
$\Theta_{\ell m}^{(0)}$, Eq. (33), and we obtain

$\displaystyle\langle\zeta(\bm{k})\Theta_{\ell
m}^{(0)}\rangle=4\pi(-i)^{\ell}Y_{\ell
m}^{*}(\hat{k})\Delta_{\ell}(k)P_{\zeta}(k).$ (85)

Integrating over the wavenumbers’ directions, and using the harmonic
decomposition of $T^{(1)}$ given in Eq. (IV.1.1), we thus arrive at

$\displaystyle\langle\Theta_{\ell_{1}m_{1}}^{(0)}\Theta_{\ell_{2}m_{2}}^{(0)}\Theta_{\ell_{3}m_{3}}^{(0)}\Theta_{\ell_{4}m_{4},\rm
inh}^{(1)}\rangle_{c}=(4\pi)^{3}f_{\rm pbh}\int D(k_{1}k_{2}k_{3})$
$\displaystyle\times
P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{3})\Delta_{\ell_{1}}(k_{1})\Delta_{\ell_{2}}(k_{2})\Delta_{\ell_{3}}(k_{3})$
$\displaystyle\times\Big{[}T_{\ell_{1}\ell_{2}\ell_{3};\ell_{4}}^{m_{1}m_{2}m_{3};m_{4}}(k_{1},k_{2},k_{3})+5\textrm{~{}perms.}\Big{]},$
(86)

where the 5 permutations involve all other possible permutations of
$k_{1},k_{2},k_{3}$ simultaneously with the corresponding permutation of the
indices $\ell_{i},m_{i},i=1,2,3$, i.e. such that the position of the index
$\ell_{i},m_{i}$ always corresponds to the position of $k_{i}$.

We now take advantage of the factorized form of
$T_{\ell_{1}\ell_{2}\ell_{3};\ell_{4}}^{m_{1}m_{2}m_{3};m_{4}}(k_{1},k_{2},k_{3})$,
given in Eqs. (64)-(66). In addition to the function $\mu_{\ell}(\eta)$
defined in Eq. (78), we define the following functions of time and multipole:

$\displaystyle\nu_{\ell}(\eta)$ $\displaystyle\equiv\int
Dk~{}P_{\zeta}(k)\Delta_{e}(\chi,k)\Delta_{\ell}(k)\frac{j_{\ell}(k\chi)}{k\chi},$
(87) $\displaystyle\lambda_{\ell}(\eta)$ $\displaystyle\equiv\int
Dk~{}P_{\zeta}(k)\Delta_{\ell}(k)\mathcal{J}_{\ell}(\eta,k).$ (88)

We then define the following one-dimensional integrals:

$\displaystyle\mathcal{A}_{\ell_{1}\ell_{2},\ell_{3}}\equiv$
$\displaystyle-2(4\pi)^{3}\int
d\eta~{}g(\eta)~{}\mu_{\ell_{1}}(\eta)\mu_{\ell_{2}}(\eta)\lambda_{\ell_{3}}(\eta),$
(89) $\displaystyle\mathcal{B}_{\ell_{1}\ell_{2},\ell_{3}}\equiv$
$\displaystyle-2(4\pi)^{3}\int
d\eta~{}g(\eta)~{}\nu_{\ell_{1}}(\eta)\nu_{\ell_{2}}(\eta)\lambda_{\ell_{3}}(\eta),$
(90)

which are symmetric in their first two arguments. We then find, using the
symmetry of $T^{(1)}$ in its first two arguments, and the symmetries of the
$Q$ and $\widetilde{Q}$ symbols defined in Eqs. (67), (68):

$\displaystyle\frac{1}{f_{\rm
pbh}}\langle\Theta_{\ell_{1}m_{1}}^{(0)}\Theta_{\ell_{2}m_{2}}^{(0)}\Theta_{\ell_{3}m_{3}}^{(0)}\Theta_{\ell_{4}m_{4},\rm
inh}^{(1)}\rangle_{c}$
$\displaystyle=\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3})}Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}+\mathcal{B}_{\ell_{1}\ell_{2},\ell_{3}}\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}$
$\displaystyle+\mathcal{B}_{\ell_{2}\ell_{3},\ell_{1}}\widetilde{Q}_{\ell_{2}\ell_{3},\ell_{1}\ell_{4}}^{m_{2}m_{3},m_{1}m_{4}}+\mathcal{B}_{\ell_{3}\ell_{1},\ell_{2}}\widetilde{Q}_{\ell_{3}\ell_{1},\ell_{2}\ell_{4}}^{m_{3}m_{1},m_{2}m_{4}},$
(91)

where222Note that we do not use the standard symmetrization notation, i.e. do
not divide by the number of terms, in order to avoid the proliferation of
numerical prefactors.

$\displaystyle\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3})}\equiv\mathcal{A}_{\ell_{1}\ell_{2},\ell_{3}}+\mathcal{A}_{\ell_{2}\ell_{3},\ell_{1}}+\mathcal{A}_{\ell_{3}\ell_{1},\ell_{2}}.$
(92)

Finally, summing over the four permutations in Eq. (IV.3), we arrive at the
main result of this work, which is the temperature trispectrum sourced by
accreting PBHs:

$\displaystyle\langle\Theta_{\ell_{1}m_{1}}\Theta_{\ell_{2}m_{2}}\Theta_{\ell_{3}m_{3}}\Theta_{\ell_{4}m_{4}}\rangle_{c}=f_{\rm
pbh}\left(\mathcal{T}_{\rm
pbh}\right)_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}},$ (93)
$\displaystyle\left(\mathcal{T}_{\rm
pbh}\right)_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}\equiv\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}$
$\displaystyle+\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}+\mathcal{B}_{\ell_{1}\ell_{3},(\ell_{2}\ell_{4})}\widetilde{Q}_{\ell_{1}\ell_{3},\ell_{2}\ell_{4}}^{m_{1}m_{3},m_{2}m_{4}}$
$\displaystyle+\mathcal{B}_{\ell_{1}\ell_{4},(\ell_{2}\ell_{3})}\widetilde{Q}_{\ell_{1}\ell_{4},\ell_{2}\ell_{3}}^{m_{1}m_{4},m_{2}m_{3}}+\mathcal{B}_{\ell_{2}\ell_{3},(\ell_{1}\ell_{4})}\widetilde{Q}_{\ell_{2}\ell_{3},\ell_{1}\ell_{4}}^{m_{2}m_{3},m_{1}m_{4}}$
$\displaystyle+\mathcal{B}_{\ell_{2}\ell_{4},(\ell_{1}\ell_{3})}\widetilde{Q}_{\ell_{2}\ell_{4},\ell_{1}\ell_{3}}^{m_{2}m_{4},m_{1}m_{3}}+\mathcal{B}_{\ell_{3}\ell_{4},(\ell_{1}\ell_{2})}\widetilde{Q}_{\ell_{3}\ell_{4},\ell_{1}\ell_{2}}^{m_{3}m_{4},m_{1}m_{2}},$
(94)

where we have defined the symmetrized coefficients

$\displaystyle\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}$
$\displaystyle\equiv\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3})}+\mathcal{A}_{(\ell_{2}\ell_{3}\ell_{4})}$
$\displaystyle\quad+\mathcal{A}_{(\ell_{3}\ell_{4}\ell_{1})}+\mathcal{A}_{(\ell_{4}\ell_{1}\ell_{2})},$
(95) $\displaystyle\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}$
$\displaystyle\equiv\mathcal{B}_{\ell_{1}\ell_{2},\ell_{3}}+\mathcal{B}_{\ell_{1}\ell_{2},\ell_{4}}.$
(96)

##  V Trispectrum constraints and sensitivity forecasts

In this section we compute and present trispectrum constraints and sensitivity
forecasts on the fraction of dark matter made of PBHs, $f_{\rm pbh}$. A full
trispectrum analysis of the _Planck_ satellite temperature data would be very
challenging and is well beyond the scope of this work. Instead, we compute the
overlap of the PBH-induced trispectrum with the local-type primordial non-
Gaussianity (PNG) trispectrum template, in order to extract an indirect limit
on the PBH abundance, given Planck’s limits on $g_{\rm NL}^{\rm loc}$ (Planck
Collaboration, 2020a). In addition, we forecast Planck’s sensitivity to the
trispectrum induced by accreting PBHs. For the scope of this paper, we ignore
biases that may arise due to lensing or other nonlinear effects, but they
should of course be accounted for in a full data analysis.

###  V.1 General equations

Given that the trispectrum induced by accreting PBHs is approximately linear
in $f_{\rm pbh}$, as given by Eq. (93), one can build an optimal quartic
estimator $\widehat{f}_{\rm pbh}$ for $f_{\rm pbh}$ Smith  _et al._ (2015);
Regan  _et al._ (2010). Its precise expression will not be needed here, and is
given in Eq. (24) of Ref. Smith  _et al._ (2015). The inverse variance of this
estimator is given by Eq. (25) in Ref. Smith  _et al._ (2015). Approximating
the noise covariance matrix as diagonal in $\ell$, the variance of the
estimator is given by

$\displaystyle\sigma_{f_{\rm pbh}}^{2}=\left\langle\mathcal{T}_{\rm
pbh}\cdot\mathcal{T}_{\rm pbh}\right\rangle^{-1},$ (97)

where for any two trispectra $\mathcal{T}_{A},\mathcal{T}_{B}$, we define
their inverse-noise weighted dot product as

$\displaystyle\left\langle\mathcal{T}_{A}\cdot\mathcal{T}_{B}\right\rangle$
$\displaystyle\equiv\frac{f_{\rm
sky}}{4!}\sum_{\ell^{\prime}s}\frac{1}{C^{\prime}_{\ell_{1}}C^{\prime}_{\ell_{2}}C^{\prime}_{\ell_{3}}C^{\prime}_{\ell_{4}}}$
$\displaystyle\times\sum_{m^{\prime}s}\left(\mathcal{T}_{A}\right)_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}\left(\mathcal{T}_{B}\right)_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}},$
(98)

where $C_{\ell}^{\prime}\equiv C_{\ell}+N_{\ell}$ is the total variance of the
observed CMB temperature, including both the cosmological signal $C_{\ell}$
and instrumental noise $N_{\ell}$, $f_{\rm sky}$ is the fraction of the sky
covered by the experiment, and the sums carry over all four indices.

Primordial non-Gaussianity also generates a CMB temperature trispectrum,
proportional to a non-Gaussianity parameter $g_{\rm NL}$:

$\langle\Theta_{\ell_{1}m_{1}}\Theta_{\ell_{2}m_{2}}\Theta_{\ell_{3}m_{3}}\Theta_{\ell_{4}m_{4}}\rangle_{c}=g_{\rm
NL}~{}\left(\mathcal{T}_{\rm
png}\right)_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}.$ (99)

One can build an optimal estimator $\widehat{g}_{\rm NL}$ for $g_{\rm NL}$ in
the same way as $f_{\rm pbh}$. The non-Gaussianity sourced by inhomogeneously
accreting PBHs would lead to a systematic bias in this estimator, even in the
absence of primordial non-Gaussianity. This bias is linear in $f_{\rm pbh}$:

$\displaystyle\langle\Delta\widehat{g}_{\rm NL}\rangle_{\rm pbh}$
$\displaystyle=f_{\rm pbh}~{}\mathcal{R},$ (100) $\displaystyle\mathcal{R}$
$\displaystyle\equiv\sigma_{g_{\rm NL}}^{2}\left\langle\mathcal{T}_{\rm
pbh}\cdot\mathcal{T}_{\rm png}\right\rangle,$ (101)

where $\sigma_{g_{\rm NL}}^{2}$ is the variance of the quadratic estimator
$\widehat{g}_{\rm NL}$, given by

$\displaystyle\sigma_{g_{\rm NL}}^{2}$
$\displaystyle\equiv\left\langle\mathcal{T}_{\rm png}\cdot\mathcal{T}_{\rm
png}\right\rangle^{-1}.$ (102)

Constraints on the amplitude $g_{\rm NL}$ of primordial non-Gaussianity
therefore directly translate into bounds on the PBH abundance $f_{\rm pbh}$.
In what follows we will specifically consider the local-type primordial
trispectrum, which is most tightly constrained by CMB anisotropy observations,
and whose shape is given in (Smith  _et al._ , 2015),

$\displaystyle\left(\mathcal{T}^{\rm loc}_{\rm
png}\right)_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}$
$\displaystyle=\mathcal{C}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}},$
(103) $\displaystyle\mathcal{C}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}$
$\displaystyle\equiv\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3},\ell_{4}}+3~{}\textrm{perm},$
(104) $\displaystyle\mathcal{C}_{\ell_{1}\ell_{2}\ell_{3},\ell_{4}}$
$\displaystyle\equiv 6\int
r^{2}dr~{}\beta_{\ell_{1}}(r)\beta_{\ell_{2}}(r)\beta_{\ell_{3}}(r)\alpha_{\ell_{4}}(r),$
(105)

where $Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}$ is given
by Eq. (67) and we used the standard notation of Refs. Smith  _et al._ (2015);
Komatsu  _et al._ (2005):

$\displaystyle\alpha_{\ell}(r)$ $\displaystyle\equiv\frac{5}{3}(4\pi)\int
Dk~{}\Delta_{\ell}(k)j_{\ell}(kr),$ (106) $\displaystyle\beta_{\ell}(r)$
$\displaystyle\equiv\frac{3}{5}(4\pi)\int
Dk~{}\Delta_{\ell}(k)j_{\ell}(kr)P_{\zeta}(k).$ (107)

###  V.2 Sums over $m$’s

Before proceeding with the numerical evaluation of Eqs. (97) and (100), we
first simplify the sums over $m$’s, which involve purely geometric quantities.
Specifically, we define

$\displaystyle(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}\equiv$
$\displaystyle\sum_{m^{\prime}s}\left(Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}\right)^{2},$
(108)
$\displaystyle(\mathcal{Q\widetilde{Q}})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}\equiv$
$\displaystyle\sum_{m^{\prime}s}{Q^{*}}_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}},$
(109)
$\displaystyle(\mathcal{\widetilde{Q}}^{2})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}\equiv$
$\displaystyle\sum_{m^{\prime}s}\left(\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}\right)^{2},$
(110) $\displaystyle(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
T})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}\equiv$
$\displaystyle\sum_{m^{\prime}s}\widetilde{Q^{*}}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}\widetilde{Q}_{\ell_{3}\ell_{4},\ell_{1}\ell_{2}}^{m_{3}m_{4},m_{1}m_{2}},$
(111) $\displaystyle(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
S})_{\ell_{1},\ell_{2}\ell_{3},\ell_{4}}\equiv$
$\displaystyle\sum_{m^{\prime}s}\widetilde{Q^{*}}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}\widetilde{Q}_{\ell_{1}\ell_{3},\ell_{2}\ell_{4}}^{m_{1}m_{3},m_{2}m_{4}},$
(112)

where “T” stands for transpose and “S” for “scrambled”. The same symmetry
rules apply where each set of indices divided by or surrounded by commas are
symmetric. We simplify these quantities in Appendix D, where we reduce them to
a single sum of products of 3-J symbols.

Inserting Eq. (94) into Eq. (97), carrying out the sums over $m$’s, and
simplifying, the inverse variance of the estimator $\widehat{f}_{\rm pbh}$
becomes

$\displaystyle\left(\sigma_{f_{\rm pbh}}^{2}\right)^{-1}$
$\displaystyle=\frac{f_{\rm
sky}}{4!}\sum_{\ell^{\prime}s}\frac{1}{C^{\prime}_{\ell_{1}}C^{\prime}_{\ell_{2}}C^{\prime}_{\ell_{3}}C^{\prime}_{\ell_{4}}}$
$\displaystyle\times\Big{[}~{}~{}~{}~{}~{}\left(\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}\right)^{2}(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}$
$\displaystyle~{}~{}~{}~{}+12~{}\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}(\mathcal{Q\widetilde{Q}})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}$
$\displaystyle~{}~{}~{}~{}+6~{}\left(\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}\right)^{2}(\mathcal{\widetilde{Q}}^{2})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}$
$\displaystyle~{}~{}~{}~{}+6~{}\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}\mathcal{B}_{\ell_{3}\ell_{4},(\ell_{1}\ell_{2})}(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
T})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}$
$\displaystyle~{}~{}~{}~{}+24~{}\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}\mathcal{B}_{\ell_{1}\ell_{3},(\ell_{2}\ell_{4})}(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
S})_{\ell_{1},\ell_{2}\ell_{3},\ell_{4}}\Big{]}.$ (113)

Similarly, the bias on local-type non-Gaussianity due to accreting PBHs
simplifies to

$\displaystyle\langle\Delta\widehat{g}^{\rm loc}_{\rm NL}\rangle_{\rm pbh}$
$\displaystyle=f_{\rm pbh}\times\frac{f_{\rm sky}}{4!}\sigma_{g_{\rm NL}^{\rm
loc}}^{2}\sum_{\ell^{\prime}s}\frac{\mathcal{C}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}}{C^{\prime}_{\ell_{1}}C^{\prime}_{\ell_{2}}C^{\prime}_{\ell_{3}}C^{\prime}_{\ell_{4}}}$
$\displaystyle\quad\times\Big{[}~{}~{}~{}~{}\mathcal{A}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}$
$\displaystyle\quad\quad+6~{}\mathcal{B}_{\ell_{1}\ell_{2},(\ell_{3}\ell_{4})}(\mathcal{Q\widetilde{Q}})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}\Big{]},$
(114)

where the inverse variance of $\widehat{g}_{\rm NL}^{\rm loc}$ is given by

$\displaystyle\left(\sigma_{g_{\rm NL}^{\rm loc}}^{2}\right)^{-1}$
$\displaystyle=\frac{f_{\rm
sky}}{4!}\sum_{\ell^{\prime}s}\frac{\left(\mathcal{C}_{(\ell_{1}\ell_{2}\ell_{3}\ell_{4})}\right)^{2}}{C^{\prime}_{\ell_{1}}C^{\prime}_{\ell_{2}}C^{\prime}_{\ell_{3}}C^{\prime}_{\ell_{4}}}(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}.$
(115)

###  V.3 Application to Planck data

We now apply the above results to the Planck experiment (Planck Collaboration,
2020b, c, a). The relevant fraction of sky coverage is $f_{\rm sky}=0.78$
(Planck Collaboration, 2020a), and the instrumental noise $N_{\ell}$ is
obtained from combining the noises of the 100, 143 and 217 GHz frequency
channels,

$\displaystyle N_{\ell}=\left[\sum_{c}N^{-1}_{\ell,c}\right]^{-1},$ (116)

where, for each channel $c$, the noise is modelled as a Gaussian with variance
per pixel $\sigma^{2}_{c}$ and beam size $\theta_{{\rm FWHM},c}$:

$\displaystyle N_{\ell,c}=\left(\frac{\sigma_{c}~{}\theta_{{\rm
FWHM},c}}{T_{0}}\right)\exp\left[\frac{\ell(\ell+1)\theta^{2}_{{\rm
FWHM},c}}{8\ln 2}\right],$ (117)

where $T_{0}=2.73$ K is the CMB monopole. The respective parameters for each
channel are333https://wiki.cosmos.esa.int/planckpla/index.php/Main_Page

$\nu_{c}$ | $\theta_{\rm FWHM,c}$ | $\sigma_{c}$
---|---|---
100 GHz | 9.66′ | 10.77 $\mu$K
143 GHz | 7.27′ | 6.40 $\mu$K
217 GHz | 5.01′ | 12.48 $\mu$K

The Planck 2018 limits on $g^{\rm loc}_{\rm NL}$ are given by (Planck
Collaboration, 2020a):

$\displaystyle g_{\rm NL}^{\rm loc}$ $\displaystyle=\left(-5.8\pm
6.5\right)\times 10^{4}\ \ \ (68\%~{}~{}\textrm{confidence})$ (118)
$\displaystyle\equiv\widehat{g_{\rm NL}^{\rm loc}}\pm\sigma_{g_{\rm NL}^{\rm
loc}}$ (119)

As a cross check of our numerical code, we compared the standard deviation of
the local-type trispectrum estimator that we obtain from Eq. (115) to the one
reported by the Planck collaboration, and given above. We find that they agree
within 5%.

To derive an indirect bound on $f_{\rm pbh}$ from the Planck constraint on
$g_{\rm NL}^{\rm loc}$, we proceed as follows. Using Bayes’ theorem, and
assuming the estimator for $g_{\rm NL}^{\rm loc}$ has a Gaussian distribution,
the un-normalized posterior probability distribution for $f_{\rm pbh}$ is
given by

$\displaystyle\mathcal{P}(f_{\rm
pbh})\propto\exp\left[-\frac{1}{2}\frac{\left(\mathcal{R}f_{\rm
pbh}-\widehat{g_{\rm NL}^{\rm loc}}\right)^{2}}{\sigma_{g_{\rm NL}^{\rm
loc}}^{2}}\right]H(f_{\rm pbh}),$ (120)

where in this context $H(x)$ designates the Heaviside function, enforcing a
positive prior on $f_{\rm pbh}$. The $(1-\epsilon)$-confidence upper limit on
$f_{\rm pbh}$, is then obtained from solving the implicit equation

$\int_{f_{\rm
pbh}}^{\infty}df~{}\mathcal{P}(f)=\epsilon~{}\int_{0}^{\infty}df~{}\mathcal{P}(f).$
(121)

###  V.4 Results and discussion

We are now fully equipped to compute upper limits on $f_{\rm pbh}$ indirectly
from the Planck constraint on $g_{\rm NL}^{\rm loc}$, and to forecast Planck’s
sensitivity to $f_{\rm pbh}$ from the temperature trispectrum. We shall
compare these limits and forecasts to Planck power-spectra limits on $f_{\rm
pbh}$. We obtain the latter with exactly the same procedure as in AK17, but
using Planck 2018 data Planck Collaboration (2020b) (instead of 2015).
Specifically, we use the foreground-marginalized Plik-lite log-likelihood for
$C_{\ell}$’s at $\ell\geq 30$, which we Taylor-expand near the Planck best-fit
cosmology, and account approximately for the low-$\ell$ data by imposing a
Gaussian prior on the optical depth to reionization. For the joint $TT,TE,EE$
limits, we use the modified version of HYREC and CLASS as implemented by AK17;
in particular we use their approximate homogeneous injection-to-deposition
Green’s function. In addition to the joint $TT,TE,EE$ limits, we also compute
a $TT$-only upper limit on $f_{\rm pbh}$ – we still retain the optical depth
prior, however, so our “$TT$-only” limits are technically temperature +
low-$\ell$ polarization limits. For a fair comparison with our $TTTT$
trispectrum limits and forecasts, for the $TT$ limit we compute the effect of
accreting PBHs at first order in $f_{\rm pbh}$, including the “direct” term
only, and using our more accurate injection-to-ionization Green’s function. We
also include the effect of inhomogeneous ionization perturbations on the
temperature power spectrum for completeness, but this makes a negligible
difference on the results.

We find that the indirect limit on $f_{\rm pbh}$ obtained from Planck’s bounds
on $g_{\rm NL}^{\rm loc}$ is systematically one order of magnitude weaker than
the $TT$-only power spectrum limit, for all PBH masses. This is due to the
weak overlap of the trispectrum induced by primordial non-Gaussianity with the
one induced by accreting PBHs: we find that the correlation coefficient of the
two shapes is less than $10\%$ across all black hole masses (using the dot
product defined in Eq. (98)). We therefore do not show this limit on our final
figure.

We show our forecasted 1-$\sigma$ sensitivity of Planck to the trispectrum of
accreting PBHs in Fig. 7, alongside current Planck power spectra upper limits
on $f_{\rm pbh}$. The upper set of curves correspond to the conservative
collisional ionization limit of AK17, while the lower set of curves correspond
to the photoionization limit (see AK17 for details). In both cases the
qualitative results are the same: we see that the temperature-only trispectrum
is not as sensitive as we had expected it to be a priori, as its sensitivity
is comparable to current $TT$ upper limits (rather than an order of magnitude
better than joint temperature and polarization limits). Nevertheless, the
temperature trispectrum is still more sensitive than temperature-only power
spectrum constraints for $M_{\rm pbh}\lesssim 10^{3}M_{\odot}$. In particular,
the temperature-only trispectrum has the potential to probe PBHs lighter by a
factor $\sim 2$ than the current reach of temperature-only power spectrum
limits.

Interestingly, the mass dependence of the trispectrum sensitivity forecast is
shallower than that of the power spectrum constraints. Moreover, we find that
making the spatial on-the-spot approximation (as described in Sec. II.3)
affects trispectrum forecasts by no more than $20\%$. Both of these features
can be explained qualitatively by the different redshift dependence of the
trispectrum and power spectrum signals, which we explore in Appendix F.

Fig. 7 also shows the updated Planck joint temperature and polarization power-
spectrum constraints ($TT,TE,EE$). We see that these constraints are tighter
than the $TT$-only constrains by about an order of magnitude. This stems from
the relatively larger effect of recombination perturbations on the
polarization signal (see e.g. Fig. 13 of AK17), indicating a stronger cross-
correlation of the perturbed CMB polarization with the unperturbed field. This
provides a strong motivation to extend our work to all temperature and
$E$-mode polarization trispectra, $TTTE,TTEE,TEEE,EEEE$, which may be
significantly more sensitive to accreting PBHs than the temperature-only
trispectrum. In addition, the inhomogeneity in the free-electron fraction
ought to induce $B$-mode polarization of magnitude comparable with the
corresponding $E$-mode polarization, $B^{(1)}_{\rm inh}\sim E^{(1)}_{\rm
inh}$. This means that, to linear order in $f_{\rm pbh}$, trispectra involving
one $B$ mode ($TTTB,TTEB,TEEB,EEEB$) ought to carry comparable signal to the
corresponding 4-point functions involving temperature and $E$-modes only.
Importantly, absent primordial tensor modes or accreting PBHs, the primary
(unlensed) CMB $B$-mode polarization vanishes. Therefore, after delensing, one
can effectively eliminate cosmic variance in the $B$-mode measurement. We thus
expect these $B$-mode trispectra to have a significantly enhanced signal-to-
noise ratio relative to their $E$-mode counterparts Meerburg  _et al._ (2016).
We defer to a future publication the extension of this work to polarization
trispectra.

Figure 7:  _Planck_ 2018 CMB power spectra constraints (solid lines) and
temperature trispectrum forecasted sensitivity (dashed red line) to the
fraction of dark matter in PBHs, as a function of PBH mass. Our forecasted
sensitivity from the temperature trispectrum is better than $TT$-only
constraints for $M_{\rm pbh}\lesssim 10^{3}M_{\odot}$ for both the collisional
ionization (thick lines) and photoionization (thin lines) limits (see AK17 for
details about these different regimes).

##  VI Conclusions

This work is the second part of a series of three papers studying the imprints
of inhomogeneously-accreting PBHs on CMB anisotropies, in particular their
higher-order statistics. The first part, Ref. (Jensen and Ali-Haïmoud, 2021),
inspected in detail how inhomogeneous energy injection from non-uniformly
accreting PBHs perturbs recombination. In the present analysis, we compute the
perturbed temperature anisotropy and its 2-point and 4-point functions. In the
upcoming third paper of this series, we will extend this work to polarization.

Our main results can be summarized as follows:
$(i)$ The inhomogeneous part of the free-electron perturbation leads to a
sub-10% effect on the perturbation to the CMB temperature power spectrum. In
other words, it is sufficient to only account for the average perturbation to
the free-electron fraction when computing the effect of accreting PBHs on the
CMB temperature power spectrum. This sub-dominant contribution was not
expected a priori and is due to the poor correlation of the perturbed CMB
temperature field with the standard temperature anisotropy. It is not
guaranteed that the same holds true for CMB polarization power spectra.
$(ii)$ We set new constraints on the PBH abundance, obtained indirectly from
Planck’s upper limits on local-type primordial non-Gaussianity. Indeed, the
shape of PBH-induced trispectrum overlaps with that of primoridal non-
Gaussianities, although weakly. This weak correlation implies that our new
constraints are not competitive with existing CMB temperature power spectrum
constraints. Still, they provide a qualitatively different probe of PBH
abundance, complementary to the usual 2-point function limits.
$(iii)$ We forecast the sensitivity of Planck to the temperature trispectrum
induced by inhomogeneously-accreting PBHs. Although our numerical results show
a weaker sensitivity than what could have been expected from simple order-of-
magnitude estimates, still we find that the temperature trispectrum would be
sensitive to PBH abundances lower than current bounds from the CMB
temperature-only 2-point function, for $M_{\rm pbh}\lesssim 10^{3}M_{\odot}$.
This is our most important result, which demonstrates that the CMB trispectrum
is indeed a useful probe of PBHs.

The calculation of higher-order CMB statistics is quite involved, and we
necessarily had to make several approximation to keep it tractable. First,
following previous studies of perturbed recombination, we only accounted for
the “direct” piece of the source term for the perturbation to CMB
anisotropies, and neglected the “feedback” piece. Unlike previous studies,
however, we explicitly quantified this approximation in the limiting case of
homogeneous ionization perturbations, and showed that it is accurate to better
than $\sim 20\%$ in that case. Still, a rigorous and definitive calculation of
the trispectrum should eventually include the “feedback” term self-
consistently. Second, we made several approximations in order to derive a
factorized quadratic transfer function for the free-electron fraction
perturbation. In particular, we conservatively approximated the injection-to-
ionization Green’s function by a factorized form that bounds it from below.
This approximation was needed to get a factorized trispectrum, much more
manageable computationally than the exact trispectrum would be. In order to
quantify the error induced by this approximation, we also considered the limit
of spatially-on-the-spot energy deposition, which bounds the injection-to-
ionization Green’s function from above. We found that all our results are
nearly unchanged when considering this limit, thus giving us confidence in
their robustness. Third, in our analysis of the primordial non-Gaussianity
bias and our trispectrum sensitivity forecast, we neglected non-Gaussianities
induced by CMB lensing. An actual analysis of CMB data should of course
correct for the lensing bias.

The most uncertain part of our calculation remains the physics of accretion
and radiation. All our numerical results rely on the semi-analytic model of
AK17 Ali-Haïmoud and Kamionkowski (2017), with a simple prescription for the
effect of relative velocities. While of course the quantitative results would
change with different assumptions about the accretion geometry and radiative
efficiency, it seems unavoidable that the PBH accretion luminosity should be
strongly modulated by large-scale supersonic relative velocities. We also
neglected entirely the effects of non-linear clustering post-recombination
Inman and Ali-Haïmoud (2019). We expect relative velocities would also
modulate the baryon content of the first halos, hence the accretion rate in
these environments. Hence, our results should still be robust qualitatively,
regardless of the details of the accretion model, or of the relevance of
accretion in non-linear halos. Moreover, the formalism we develop is quite
general and could be applied to arbitrary perturbations of recombination
spatially modulated by relative velocities, or even more generally quadratic
in initial conditions.

Even if the temperature trispectrum is not quite as sensitive to PBHs as we
had anticipated from the simple order of magnitude presented in the
introduction, our results are still very significant and promising. Indeed, we
uncovered a completely new CMB observable to probe PBHs, with a sensitivity
comparable to, and in some cases better than, current CMB temperature power-
spectrum constraints. Importantly, while several energy injection processes
could in principle mimick the effect of accreting PBHs in CMB power spectra,
to our knowledge the trispectrum signature studied in this work is unique to
them. These considerations provide strong motivation to extend this work and
study the polarization signal of inhomogeneously accreting PBHs. In addition
to trispectra involving $E$-mode polarization ($TTTE,TTEE,TEEE,EEEE$), we also
expect $B$-mode non-Gaussianity, in the form of $TTTB,TTEB,TEEB,EEEB$
trispectra at leading order in PBH abundance. These $B$-mode trispectra ought
to have amplitudes comparable to their $E$-mode counterparts, but much lower
noise. We defer the computation of these promising observables to the third
and last installment of this series of publications.

## Acknowledgements

YAH is a CIFAR-Azrieli Global scholar and acknowledges funding from CIFAR.
This work was supported in part by NSF grant No. 1820861. We thank Daan
Meerburg for suggesting to seek a factorized form for the trispectrum, and
Colin Hill for comments on this manuscript.

##  Appendix A Correlation functions involving a function of relative velocity

In what follows we denote $\bm{v}\equiv\bm{v}_{\rm bc}$ the relative velocity
of baryons and dark matter. We need to compute $(N+1)$-point functions of the
form

$\displaystyle\braket{F(\bm{k})\delta_{1}(\bm{k}_{1})\cdots\delta_{N}(\bm{k}_{N})}=P_{N}(\bm{k},\bm{k}_{1},\cdots,\bm{k}_{N})$
$\displaystyle\times(2\pi)^{3}\delta^{(3)}(\bm{k}+\bm{k}_{1}+\cdots+\bm{k}_{N}),$
(122)

where $F$ depends on position only through the magnitude $v$ of the relative
velocity field, i.e. $F(\bm{x})=F(v(\bm{x}))$, and has zero mean,
$\braket{F}=0$, and $\delta_{1},\cdots,\delta_{N}$ are scalars also with zero
mean linearly related to the primordial curvature perturbation. This
$(N+1)$-point function is non-zero only if $N$ is even, given that $F$ is an
even function of relative velocity. The $(N+1)$-spectrum $P_{N}$ is the
Fourier transform of the $(N+1)$-point correlation function

$\xi_{N}(\bm{x}_{1},\cdots,\bm{x}_{N})\equiv\braket{F(v(\bm{0}))\delta_{1}(\bm{x}_{1})\cdots\delta_{N}(\bm{x}_{N})}.$
(123)

The goal of this appendix is to derive an approximate expression for
$\xi_{N}$, from which one can also approximate $P_{N}$.

In full generality, provided $\bm{v},\delta_{1},\cdots,\delta_{N}$ are
Gaussian-distributed, we have

$\displaystyle\xi_{N}(\bm{x}_{1},\cdots,\bm{x}_{N})=\int
d^{3}v~{}d\delta_{1}\cdots d\delta_{N}~{}F(v)\delta_{1}\cdots\delta_{N}$
$\displaystyle\times\frac{1}{\sqrt{2\pi\det(\bm{C})}}\exp\left[-\frac{1}{2}\bm{X}^{\rm
T}\cdot\bm{C}^{-1}\cdot\bm{X}\right],$ (124)

where

$\displaystyle\bm{X}^{\rm
T}\equiv(\widetilde{\bm{v}},\widetilde{\bm{\delta}}^{\rm
T})\equiv\left(\frac{\bm{v}}{\sigma_{1d}},\frac{\delta_{1}}{\sigma_{\delta_{1}}},\cdots,\frac{\delta_{N}}{\sigma_{\delta_{N}}}\right),$
(125)

with $\sigma_{1d}^{2}\equiv\braket{v^{2}}/3$ and
$\sigma_{\delta_{i}}^{2}\equiv\braket{\delta_{i}^{2}}$. $\bm{C}$ is the
$(N+3)$ by $(N+3)$ normalized covariance matrix of
$\widetilde{\bm{v}}(\bm{0}),\widetilde{\delta_{1}}(\bm{x}_{1}),\cdots,\widetilde{\delta_{N}}(\bm{x}_{N})$.
Explicitly, this matrix is given by $\bm{C}=\bm{C}_{0}+\bm{\Delta}$, with

$\displaystyle\bm{C}_{0}\equiv\begin{pmatrix}\mbox{\Large 1}_{3\times
3}&\vline&\mbox{\Large 0}_{3\times N}\\\\[8.0pt] \hline\cr\mbox{\Large
0}_{N\times 3}&\vline&\bm{C}_{\widetilde{\delta}}\end{pmatrix},$ (126)
$\displaystyle\bm{\Delta}\equiv\begin{pmatrix}\mbox{\Large 0}_{3\times
3}&\vline&\begin{matrix}\bm{\Xi}_{1}^{\rm T}\\\ \bm{\Xi}_{2}^{\rm T}\\\
\bm{\Xi}_{3}^{\rm T}\end{matrix}\\\
\hline\cr\begin{matrix}\bm{\Xi}_{1}&\bm{\Xi}_{2}&\bm{\Xi}_{3}\end{matrix}&\vline&\mbox{\Large
0}_{N\times N}\end{pmatrix},$ (127)

where $\bm{C}_{\widetilde{\delta}}$ is the $N\times N$ normalized covariance
matrix of the $\widetilde{\delta}$’s, and $\bm{\Xi}_{i},i=1,2,3$, are the
$N$-dimensional column vectors

$\bm{\Xi}_{i}\equiv\begin{pmatrix}\braket{{\widetilde{v}_{i}(0)\widetilde{\delta}_{1}(\bm{x}_{1})}}\\\
\vdots\\\
\braket{\widetilde{v}_{i}(0)\widetilde{\delta}_{N}(\bm{x}_{N})}\end{pmatrix}.$
(128)

In words, the matrix $\bm{C}_{0}$ includes all correlations except for the
velocity-$\delta$ correlations, which are included in $\bm{\Delta}$.

So far, these expressions are exact. We expect that, in general, $\bm{\Delta}$
is small for _any_ separation. Indeed, this is always true in the large-
separation limit. Moreover, statistical isotropy implies that
$\braket{\bm{v}(0)\delta(\bm{x})}\rightarrow 0$ when $\bm{x}\rightarrow 0$,
since there is no non-null isotropic rank-1 tensor. This can be seen in Fig. 8
where we correlate $v_{i}$ with the canonical monopoles of the $\Theta^{(0)}$
line-of-sight source transfer functions of ${S}^{(0)}$ (c.f. the first line of
Eq. (III.2)) for example.

Figure 8:  Correlation function of the monopole terms of the line-of-sight
source for $\Theta^{(0)}$ near recombination with the relative velocity
between CDM and baryons at various redshifts. Namely we plot the variance
normalized correlation function $\widetilde{c}_{r}$, where $\langle{\bf
v}(z){S^{(0)}_{0}}(z^{\prime}=1100)\rangle\equiv c_{r}(r,z)\hat{r}$, and the
subscript implies the monopole terms only. Even at intermediate scales it is
less than unity and justifies expanding the covariance matrix discussed in
Appendix A.

We may therefore expand $\bm{C}^{-1}$ around $\bm{C}_{0}^{-1}$ to compute
$\xi_{N}$. We’ll see that it is required to include terms at second order in
$\bm{\Delta}$:

$\bm{C}^{-1}=\bm{C}_{0}^{-1}-\bm{C}_{0}^{-1}\bm{\Delta}\bm{C}_{0}^{-1}+\bm{C}_{0}^{-1}\bm{\Delta}\bm{C}_{0}^{-1}\bm{\Delta}\bm{C}_{0}^{-1}+\mathcal{O}(\Delta^{3}),$
(129)

from which we get, at second order in $\Delta$,

$\displaystyle\exp\left[-\frac{1}{2}\bm{X}^{\rm
T}\bm{C}^{-1}\bm{X}\right]=\Lambda\times\exp\left[-\frac{1}{2}\bm{X}^{\rm
T}\bm{C}_{0}^{-1}\bm{X}\right],$ $\displaystyle\Lambda\equiv
1+\frac{1}{2}\bm{\tilde{X}}^{\rm T}\bm{\Delta}\bm{\tilde{X}}$
$\displaystyle~{}~{}+\frac{1}{8}\left(\bm{\tilde{X}}^{\rm
T}\bm{\Delta}\bm{\tilde{X}}\right)^{2}-\frac{1}{2}\bm{\tilde{X}}^{\rm
T}\bm{\Delta}\bm{C}_{0}^{-1}\bm{\Delta}\bm{\tilde{X}},~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(130)

with $\bm{\tilde{X}}\equiv\bm{C}_{0}^{-1}\bm{X}$. With this approximation, we
thus have

$\xi_{N}\approx\sqrt{\frac{\det(\bf{C}_{0})}{\det(\bf{C})}}\langle
F(v)\delta_{1}\cdots\delta_{N}\Lambda\rangle_{0}\approx\langle
F(v)\delta_{1}\cdots\delta_{N}\Lambda\rangle_{0},$ (131)

where the average $\langle\cdots\rangle_{0}$ is over the “unperturbed”
$(N+3)$-D Gaussian distribution with covariance matrix $\bf{C}_{0}$, which is
the product of two uncorrelated Gaussian distributions: an isotorpic Gaussian
distribution for $\bm{v}(\bm{0})$ and a $N$-dimensional Gaussian for
$(\delta_{1}(\bm{x}_{1}),\cdots,\delta_{N}(\bm{x}_{N}))$, with covariance
matrix $\bm{C}_{\delta}$. The second equality is valid to lowest order in
$\Delta$.

Upon integrating over velocities, the contribution of the first term in
$\Lambda$ (i.e. 1) vanishes, since $\braket{F(v)}=0$. Let us now compute the
other terms. First, let us compute

$\tilde{\bm{X}}^{\rm T}\bm{\Delta}=\left(\bm{\delta}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{1},\bm{\delta}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{2},\bm{\delta}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{3},\frac{v_{i}}{\sigma_{1\rm
d}^{2}}\bm{\Xi}_{i}^{\rm T}\right),$ (132)

where the first three terms are scalars, and the last term contains an
implicit sum over $i$, and is a $N$-dimensional vector. We therefore have

$\Lambda_{1}\equiv\bm{\tilde{X}}^{\rm
T}\bm{\Delta}\bm{\tilde{X}}=\frac{2v_{i}}{\sigma_{\rm
1d}^{2}}\bm{\Xi}_{i}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\delta}=\frac{2v_{i}}{\sigma_{\rm
1d}^{2}}\bm{\delta}^{\rm T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{i}.$ (133)

The second term $\Lambda_{1}$ is therefore linear in $v_{i}$. Therefore
$\langle F(v)\delta_{1}...\delta_{N}\Lambda_{1}\rangle_{0}=0$ since $\langle
v_{i}F(v)\rangle=0$, by isotropy.

We thus need to only include the third and last terms in $\Lambda$, quadratic
in $\Delta$.

Let’s start with the third term, proportional to $\Lambda_{1}^{2}$. From our
previous results, we have

$\Lambda_{1}^{2}=\frac{4}{\sigma_{\rm 1d}^{4}}v_{i}v_{j}\bm{\Xi}_{i}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\delta}\bm{\delta}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{j},$ (134)

where repeated indices are summed over. Using $\langle
F(v)v_{i}v_{j}\rangle_{0}=\frac{1}{3}\delta_{ij}\langle v^{2}F(v)\rangle_{0}$,
we thus find

$\displaystyle\Big{\langle}F(v)\delta_{1}\cdots\delta_{N}\Lambda_{1}^{2}\Big{\rangle}_{0}=\frac{4}{3\sigma_{1d}^{4}}\langle
v^{2}F(v)\rangle$ $\displaystyle\times\bm{\Xi}_{i}^{\rm
T}\bm{C}_{\delta}^{-1}\langle\delta_{1}\cdots\delta_{N}\bm{\delta}\bm{\delta}^{\rm
T}\rangle_{0}\bm{C}_{\delta}^{-1}\bm{\Xi}_{i}.$ (135)

On the other hand, we have

$\displaystyle\Lambda_{2}$ $\displaystyle\equiv$
$\displaystyle\bm{\tilde{X}}^{\rm
T}\bm{\Delta}\bm{C}_{0}^{-1}\bm{\Delta}\bm{\tilde{X}}=\frac{v_{i}v_{j}}{\sigma_{1d}^{4}}\bm{\Xi}_{i}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{j}$ (136) $\displaystyle+$
$\displaystyle\textrm{terms independent of $v$}.$

This implies

$\displaystyle\Big{\langle}F(v)\delta_{1}\cdots\delta_{N}\Lambda_{2}\Big{\rangle}_{0}=\frac{\braket{v^{2}F(v)}}{3\sigma_{1d}^{4}}\langle\delta_{1}\cdots\delta_{N}\rangle\bm{\Xi}_{i}^{\rm
T}\bm{C}_{\delta}^{-1}\bm{\Xi}_{i}.$ (137)

Therefore, combining terms, we obtain

$\displaystyle\xi_{N}(\bm{x}_{1},\cdots,\bm{x}_{N})$ $\displaystyle\approx$
$\displaystyle\frac{1}{8}\Big{\langle}F(v)\delta_{1}\cdots\delta_{N}\Lambda_{1}^{2}\Big{\rangle}_{0}$
(138) $\displaystyle-$
$\displaystyle\frac{1}{2}\Big{\langle}F(v)\delta_{1}\cdots\delta_{N}\Lambda_{2}\Big{\rangle}_{0}$
$\displaystyle=$
$\displaystyle\braket{v^{2}F(v)}S(\bm{x}_{1},\cdots,\bm{x}_{N}),$

where we have defined

$S\equiv\frac{1}{6\sigma_{1d}^{4}}\bm{\Xi}_{i}^{\rm
T}\bm{C}_{\delta}^{-1}\Big{\langle}\delta_{1}\cdots\delta_{N}\left(\bm{\delta\delta}^{\rm
T}-\bm{C}_{\delta}\right)\Big{\rangle}\bm{C}_{\delta}^{-1}\bm{\Xi}_{i}.$ (139)

We see that in this approximation, the shape of the $N$-point correlation is
entirely determined by $S(\bm{x}_{1},\cdots,\bm{x}_{N})$, regardless of the
function $F(v)$. The latter only affects the overall amplitude of the
correlation function, and only through its moment $\langle v^{2}F(v)\rangle$.

Therefore, to compute the $N$-point correlation function, one may substitute
$F(v)$ with a simpler function $\tilde{F}(v)$, as long as $\langle
v^{2}\tilde{F}(v)\rangle=\langle v^{2}F(v)\rangle$. The simplest such function
is $\tilde{F}(v)\equiv b_{F}\left(\frac{v^{2}}{3\sigma_{1d}^{2}}-1\right)$. It
is such that

$\langle v^{2}\tilde{F}(v)\rangle=2b_{F}\sigma_{1\rm d}^{2}.$ (140)

Hence, we may use $\tilde{F}(v)$ instead of $F(v)$ provided the parameter
$b_{F}$ is given by

$b_{F}=\frac{1}{2\sigma_{1d}^{2}}\braket{v^{2}F(v)}.$ (141)

This result was proven in configuration space but also holds in the Fourier
domain, where it is most useful: we have proven that, for any $N$-point
function involving $N$ scalar functions (provided the $\delta v_{i}$
correlations are sufficiently small at all separations), we may use
$\tilde{F}(v)=b_{F}\left(\frac{v^{2}}{3\sigma_{1d}^{2}}-1\right)$ in order to
compute $N$-point functions. Importantly, this means that the shape we derive
for the trispectrum should be relatively insensitive to the details of
accretion physics – the shape still has some dependence on it, as in practice
the bias parameter $b_{F}$ is redshift-dependent, in a way that depends on the
details of accretion.

##  Appendix B Numerical resolution and convergence

In this appendix we describe our sampling of $\eta,k$ and $\ell$ integrals and
sums.

Because each conformal-time integral relevant to the PBH-induced trispectrum
includes the visibility function $g(\eta)$, we sample $\eta$ more finely
during recombination. Starting from $z_{\rm max}=1400$, we sample $\eta$ with
logarithmic step size $\Delta\ln\eta=10^{-3}$ until $z_{\rm rec}=900$, after
which we increase the step size to $\Delta\ln\eta=2\times 10^{-2}$ until
$z_{\rm re}=10$, and then finally we sample linearly in $\eta$ until $z=0$
with step size $\Delta\eta=50$ Mpc.

For $k$ integrals, we compute quantities on a grid from $k_{\rm min}=10^{-5}$
Mpc-1 up to a maximum wave number $k_{\rm max}=5000\eta_{0}^{-1}$ with a step
size $\Delta k={\rm min}(\epsilon k,\kappa_{0})$, where $\epsilon=0.006$ and
$\kappa_{0}=10^{-4}$ Mpc-1, i.e. use logarithm spacing for low-$k$ to linear
spacing at high-$k$.

Finally, our $\ell$ sampling consists of the floors of an array of real $\ell$
values spaced logarithmically in $2\leq\ell<400$ with $\Delta\ln\ell=0.0225$,
and linearly in $400\leq\ell<\ell_{\max}=3000$ with $\Delta\ell=19.5$. Note
that these values were chosen to produce an $\ell$ sampling similar to the
standard output of CLASS, with almost double the resolution.

We reproduce the standard CMB temperature angular power spectrum, Eq. (36),
and compare to the output of CLASS (Blas  _et al._ , 2011). We find a sub-
percent fractional difference for all $\ell<\ell_{\rm max}$. We also recompute
all results with increased resolution prescribed via,

$\displaystyle(\Delta\ln\eta,\,\Delta\eta,\,k_{\rm max}$
$\displaystyle,\epsilon,\kappa_{0})\rightarrow$
$\displaystyle\left(\frac{2}{3}\Delta\ln\eta,\,\frac{2}{3}\Delta\eta,\,\frac{3}{2}k_{\rm
max},\frac{2}{3}\epsilon,\frac{2}{3}\kappa_{0}\right).$ (142)

We find, for both the power spectrum and trispectrum calculations, there is a
fractional change in the results only at sub-percent level, far below the
theoretical uncertainty of the problem at hand.

Lastly, the trispectrum results depend on the intermediate quantity
$\mathcal{J}_{\ell}$, given in Eq. (54) as an infinite double sum. We truncate
this sum at a maximum $\ell_{2}=\ell_{\rm cut}$ (which automatically truncates
the $\ell_{1}$ sum due to the triangle inequality). We find that our
trispectrum results are converged within $0.1\%$ by $\ell_{\rm cut}=50$.

##  Appendix C Spin-weighted spherical harmonics

Spin-weighted spherical harmonics are related to Wigner $D$-matrices. They
become regular spherical harmonics when their spin is zero, $(_{0}Y_{\ell
m})=Y_{\ell m}$ and inherit similar orthogonal and completeness relations.
They have the familiar property $(_{s}Y_{\ell
m})^{*}=(-1)^{s+m}(_{-s}Y_{\ell-m})$, as well as a product rule similar to the
Gaunt relation involving Wigner 3$j$ symbols (Varshalovich  _et al._ , 1988),

$\displaystyle~{}_{s_{1}}Y_{\ell_{1}m_{1}}(\hat{n})~{}_{s_{2}}Y_{\ell_{2}m_{2}}(\hat{n})=\sum_{s_{3},\ell_{3}m_{3}}g^{-s_{1}(-s_{2})(-s_{3})}_{\ell_{1}\ell_{2}\ell_{3}}$
$\displaystyle\times{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}}~{}_{s_{3}}Y^{*}_{\ell_{3}m_{3}}(\hat{n}),$
(143)

where the $g$-symbols are defined by

$\displaystyle g_{\ell_{1}\ell_{2}\ell_{3}}^{s_{1}s_{2}s_{3}}$
$\displaystyle\equiv\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell_{3}+1)}{4\pi}}{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
s_{1}&s_{2}&s_{3}\end{pmatrix}},$ $\displaystyle g_{\ell_{1}\ell_{2}\ell_{3}}$
$\displaystyle\equiv g_{\ell_{1}\ell_{2}\ell_{3}}^{0\,0\,0},$ (144)

For shorthand we also define the Gaunt coefficient,

$\displaystyle\mathcal{G}^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}\equiv\,$
$\displaystyle
g_{\ell_{1}\ell_{2}\ell_{3}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}.$ (145)

For the Wigner $3j$ symbols to be nonzero, the $\ell$’s in the first row must
be positive and obey the triangle inequality. Likewise, the sum of the bottom
row of azimuthal modes ($m_{1},m_{2},m_{3}$) must equate to zero, and each
must satisfy $-\ell_{i}\leq m_{i}\leq\ell_{i}$. The Wigner 3$j$ symbols also
have an orthogonality condition that we utilize,

$\displaystyle\sum_{m_{1}m_{2}}{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}}$
$\displaystyle{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell_{3}^{\prime}\\\
m_{1}&m_{2}&m_{3}^{\prime}\end{pmatrix}}=$
$\displaystyle\quad\quad\frac{\delta_{\ell_{3}\ell_{3}^{\prime}}\delta_{m_{3}m_{3}^{\prime}}}{2\ell_{3}+1}\\{\ell_{1}\ell_{2}\ell_{3}\\},$
(146)

where $\\{\ell_{1}\ell_{2}\ell_{3}\\}$ is 1 if the three $\ell$’s satisfy the
triangle inequality and 0 otherwise.

When summing over the azimuthal modes of the product of spin-weighted
spherical harmonics, we have,

$\displaystyle\sum_{m}(_{s}Y_{\ell m}(\hat{n}))(_{s^{\prime}}Y_{\ell
m}(\hat{n}^{\prime}))^{*}=(-1)^{s}\frac{2\ell+1}{4\pi}d^{\ell}_{ss^{\prime}}\left(\mu\right),$
(147)

where we have introduced the Wigner small $d$-functions and
$\mu\equiv\hat{n}\cdot\hat{n}^{\prime}$ (Smith  _et al._ , 2015). If
$s=s^{\prime}=0$, then the $d$-functions reduce to normal Legendre
polynomials. These $d$-functions themselves satisfy the orthogonality
condition,

$\displaystyle\int^{1}_{-1}d\mu~{}d^{\ell_{1}}_{ss^{\prime}}(\mu)d^{\ell_{2}}_{ss^{\prime}}(\mu)=\frac{2}{2\ell+1}\delta_{\ell_{1}\ell_{2}},$
(148)

and are equipped with the identity,

$\displaystyle
d^{\ell}_{-s(-s^{\prime})}(\mu)=d^{\ell}_{s^{\prime}s}(\mu)=(-1)^{s+s^{\prime}}d^{\ell}_{ss^{\prime}}(\mu).$
(149)

They also have the property that their product can be expanded via
Varshalovich  _et al._ (1988),

$\displaystyle
d^{\ell_{1}}_{s_{1}s_{2}^{\prime}}(\mu)d^{\ell_{2}}_{s_{2}s_{2}^{\prime}}(\mu)=$
$\displaystyle\quad\quad\sum_{\ell,s,s^{\prime}}(2\ell+1){\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
s_{1}&s_{2}&s\end{pmatrix}}d^{\ell}_{ss^{\prime}}(\mu){{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
s_{1}^{\prime}&s_{2}^{\prime}&s^{\prime}\end{pmatrix}}}.$ (150)

##  Appendix D Sums of products of $Q$ and $\widetilde{Q}$ symbols

Computing the multipoles of the nonlinear perturbation of temperature
anisotropy introduces integrals of products of four (spin-weighted) spherical
harmonics, denoted as the $Q$-symbols in Eq. (67) and (68). In this appendix
we lay out the math to simplify the sums and products of these $Q$-symbols
necessary for the first-order trispectrum calculations. We borrow the tools
introduced in Appendix C.

Lest us start with $(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}$
defined in Eq. (108). Given the definitions of
$Q_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}^{m_{1}m_{2}m_{3}m_{4}}$, it is given by

$(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}=\sum_{m^{\prime}s}\int
d^{2}\hat{n}\int
d^{2}\hat{n}^{\prime}\prod_{i}Y_{\ell_{i}m_{i}}(\hat{n})Y^{*}_{\ell_{i}m_{i}}(\hat{n}^{\prime}).$
(151)

Let us now use Eq. (147), which reduces to Legendre polynomials in this case:

$\sum_{m_{i}}Y_{\ell_{i}m_{i}}(\hat{n})Y^{*}_{\ell_{i}m_{i}}(\hat{n}^{\prime})=\frac{1}{4\pi}(2\ell_{i}+1)P_{\ell_{i}}(\mu),$
(152)

where $\mu\equiv\hat{n}\cdot\hat{n}^{\prime}$. We may carry out one of the
angular integrals and get

$\displaystyle(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}$
$\displaystyle=\frac{1}{(4\pi)^{2}}\frac{1}{2}\prod_{i}(2\ell_{i}+1)$
$\displaystyle\times\int_{-1}^{1}d\mu~{}P_{\ell_{1}}(\mu)P_{\ell_{2}}(\mu)P_{\ell_{3}}(\mu)P_{\ell_{4}}(\mu).$
(153)

Now, recall the product rule for Wigner $d$-functions (or Legendre polynomials
in this case), Eq (C),

$P_{\ell_{1}}P_{\ell_{2}}=\sum_{\ell}(2\ell+1){{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
0&0&0\end{pmatrix}}}^{2}P_{\ell}.$ (154)

Therefore, using the orthogonality relation Eq. (148), we obtain

$\displaystyle(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell_{4}}=$
$\displaystyle\frac{1}{(4\pi)^{2}}\prod_{i}(2\ell_{i}+1)\sum_{\ell}(2\ell+1)$
$\displaystyle\quad\quad\times{{\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
0&0&0\end{pmatrix}}}^{2}{{\begin{pmatrix}\ell&\ell_{3}&\ell_{4}\\\
0&0&0\end{pmatrix}}}^{2}$ $\displaystyle=$
$\displaystyle\sum_{\ell}\frac{1}{2\ell+1}(g_{\ell_{1}\ell_{2}\ell})^{2}(g_{\ell\ell_{3}\ell_{4}})^{2},$
(155)

with a result that should be independent of the grouping of the two pairs of
$\ell$’s.

To generalize this to sums involving $\widetilde{Q}$ is relatively
straightforward, but requires general Wigner small $d$-functions. We note that
the angular derivatives present in $\widetilde{Q}$ can be expressed in terms
of spin-1 spin-weighted spherical harmonics (Smith  _et al._ , 2015). That is,

$\displaystyle\bm{\nabla}_{\hat{n}}Y_{\ell_{1}m_{1}}\cdot\bm{\nabla}_{\hat{n}}Y_{\ell_{2}m_{2}}$
$\displaystyle=-\frac{1}{2}\sqrt{\ell_{1}(\ell_{1}+1)\ell_{2}(\ell_{2}+1)}$
$\displaystyle\times\sum_{s=\pm
1}~{}_{s}Y_{\ell_{1}m_{1}}~{}_{-s}Y_{\ell_{2}m_{2}},$ (156)

such that,

$\displaystyle\widetilde{Q}_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}^{m_{1}m_{2},m_{3}m_{4}}=-\frac{1}{2}\sqrt{\ell_{1}(\ell_{1}+1)\ell_{2}(\ell_{2}+1)}$
$\displaystyle\quad\sum_{s=\pm 1}\int
d^{2}\hat{n}~{}_{s}Y^{*}_{\ell_{1}m_{1}}(\hat{n})~{}_{-s}Y^{*}_{\ell_{2}m_{2}}(\hat{n})Y^{*}_{\ell_{3}m_{3}}(\hat{n})Y_{\ell_{4}m_{4}}^{*}(\hat{n}).$
(157)

For short we define

$\widetilde{g}_{\ell_{1}\ell_{2},\ell}\equiv\sum_{s=\pm
1}g_{\ell_{1},\ell_{2},\ell}^{s,-s,0},$ (158)

which is symmetric in its first two indices. Using the properties of
$d$-functions outlined in Appendix C, we obtain for Eq. (108) – (112),

$\displaystyle(\mathcal{Q\widetilde{Q}})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}=-\frac{1}{2}\sqrt{\ell_{1}(\ell_{1}+1)\ell_{2}(\ell_{2}+1)}$
$\displaystyle\quad\quad\times\sum_{\ell}\frac{1}{2\ell+1}g_{\ell_{1}\ell_{2}\ell}~{}\widetilde{g}_{\ell_{1}\ell_{2},\ell}(g_{\ell_{3}\ell_{4}\ell})^{2},$
(159)
$\displaystyle(\mathcal{\widetilde{Q}}^{2})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}=\frac{1}{4}\ell_{1}(\ell_{1}+1)\ell_{2}(\ell_{2}+1)$
$\displaystyle\quad\quad\times\sum_{\ell}\frac{1}{2\ell+1}(\widetilde{g}_{\ell_{1}\ell_{2},\ell})^{2}(g_{\ell\ell_{3}\ell_{4}})^{2},$
(160) $\displaystyle(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
T})_{\ell_{1}\ell_{2},\ell_{3}\ell_{4}}=\frac{1}{4}\sqrt{\ell_{1}(\ell_{1}+1)\ell_{2}(\ell_{2}+1)}$
$\displaystyle\quad\quad\times\sqrt{\ell_{3}(\ell_{3}+1)\ell_{4}(\ell_{4}+1)}$
$\displaystyle\quad\quad\times\sum_{\ell}\frac{1}{2\ell+1}(g_{\ell_{1}\ell_{2}\ell}~{}\widetilde{g}_{\ell_{1}\ell_{2},\ell})(g_{\ell_{3}\ell_{4}\ell}~{}\widetilde{g}_{\ell_{3}\ell_{4},\ell}),$
(161) $\displaystyle(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
S})_{\ell_{1},\ell_{2}\ell_{3},\ell_{4}}=-\frac{1}{4}\ell_{1}(\ell_{1}+1)\sqrt{\ell_{2}(\ell_{2}+1)\ell_{3}(\ell_{3}+1)}$
$\displaystyle\quad\times\sum_{s=\pm
1}\sum_{\ell}\frac{1}{2\ell+1}\widetilde{g}_{\ell_{1}\ell_{2},\ell}g_{\ell,\ell_{1},\ell_{2}}^{s,-s,0}~{}g_{\ell,\ell_{3},\ell_{4}}^{-s,s,0}~{}g_{\ell\ell_{3}\ell_{4}}.$
(162)

##  Appendix E Perturbed temperature anisotropy auto-power-spectrum

In this paper we found that the temperature-only trispectrum induced by
accreting PBHs was not as sensitive as we expected a priori. Additionally, the
amplitude of the power spectrum perturbation sourced by inhomogeneities in the
free-electron fraction, $C^{(1)}_{\ell,\rm inh}\equiv
2\langle\Theta^{(1)}_{\ell m,\rm inh}\Theta^{(0)*}_{\ell m}\rangle$ is up to
two orders of magnitude smaller than its counterpart $C^{(1)}_{\ell,\rm
hom}\equiv 2\langle\Theta^{(1)}_{\ell m,\rm hom}\Theta^{(0)*}_{\ell
m}\rangle$, as revealed in Fig. 6. In this appendix, we show this is due to a
combination of a poor correlation between $\Theta^{(1)}_{\rm inh}$ and the
standard CMB temperature anisotropy $\Theta^{(0)}$, and a suppression of the
characteristic amplitude of $\Theta^{(1)}_{\rm inh}$ itself, relative to its
counterpart $\Theta^{(1)}_{\rm hom}$. We do so by computing and comparing the
auto-power-spectra of $\Theta^{(1)}_{\rm hom}$ and $\Theta^{(1)}_{\rm inh}$.
The results are shown in Fig. 9.

From Eq. (47), the auto-power-spectrum of $\Theta^{(1)}_{\rm hom}$ is
trivially

$\displaystyle C_{\ell,\rm hom}^{(11)}=4\pi\int
Dk~{}P_{\zeta}(k)\left[\Delta_{\ell,\rm hom}^{(1)\rm d}(k)\right]^{2},$ (163)

where $\langle\Theta^{(1)}_{\ell m,\rm
hom}\Theta^{*(1)}_{\ell^{\prime}m^{\prime},\rm
hom}\rangle\equiv\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}C_{\ell,\rm
hom}^{(11)}$.

The auto-power-spectrum of the inhomogeneous-ionization counterpart,
$\langle\Theta^{(1)}_{\ell m,\rm
inh}\Theta^{*(1)}_{\ell^{\prime}m^{\prime},\rm
inh}\rangle\equiv\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}C_{\ell,\rm
inh}^{(11)}$, is much more involved. In what follows, we denote the integral
operator,

$\displaystyle\int
D(k_{1}k_{2}k_{3})P_{\zeta}(k_{1})P_{\zeta}(k_{2})P_{\zeta}(k_{3})\equiv\int\mathcal{D}^{3}\\!P.$
(164)

We begin similarly as we did for the trispectrum calculation in Sec. IV.3.
Starting with Eq. (58), using Wick’s theorem, and exploiting the fact that
$T^{(1)}_{\ell m}(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{3})=0$ and $T^{(1)}_{\ell m}$
is symmetric in its first two ${\bm{k}}$ arguments, we find,

$\displaystyle\langle$ $\displaystyle\Theta_{\ell m,\rm
inh}^{(1)}\Theta_{\ell^{\prime}m^{\prime},\rm inh}^{*(1)}\rangle=$
$\displaystyle\quad\quad f_{\rm
pbh}^{2}\int\mathcal{D}^{3}\\!P\left\\{4T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},-\bm{k}_{2})T_{\ell^{\prime}m^{\prime}}^{*(1)}(\bm{k}_{1},\bm{k}_{3},-\bm{k}_{3})\right.$
$\displaystyle\quad\quad\quad\quad+2T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})T_{\ell^{\prime}m^{\prime}}^{*(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle\quad\quad\quad\quad+4\left.T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})T_{\ell^{\prime}m^{\prime}}^{*(1)}(\bm{k}_{3},\bm{k}_{2},\bm{k}_{1})\right\\}.$
(165)

The first term can be solved with the same method as for the inhomogeneous
power spectrum in Sec. IV.2. That is,

$\displaystyle\int\mathcal{D}^{3}\\!P$ $\displaystyle~{}T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},-\bm{k}_{2})T_{\ell^{\prime}m^{\prime}}^{*(1)}(\bm{k}_{1},\bm{k}_{3},-\bm{k}_{3})=$
$\displaystyle\quad\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\frac{16\pi}{9}\int_{0}^{\eta_{0}}d\eta\int_{0}^{\eta_{0}}d\eta^{\prime}g(\eta)g(\eta^{\prime})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\times\mathcal{A}_{\ell}(\eta,\eta^{\prime})\gamma(\eta)\gamma(\eta^{\prime}),$
(166)

where $\gamma(\eta)$ is defined in Eq. (77) and,

$\displaystyle\mathcal{A}_{\ell}(\eta,\eta^{\prime})$
$\displaystyle\equiv\int\\!DkP_{\zeta}(k)\Delta_{e}(\eta,k)j_{\ell}^{\prime}(k\chi)\Delta_{e}(\eta^{\prime},k)j_{\ell}^{\prime}(k\chi^{\prime}).$
(167)

The remaining two terms are not as simple. Using Eq. (IV.1.1) and integrating
over all three $\hat{k}$’s (but restoring them for notational convenience by
absorbing factors of $4\pi$), the second term can be written as,

$\displaystyle\int\mathcal{D}^{3}\\!P~{}T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})T_{\ell^{\prime}m^{\prime}}^{*(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=$
$\displaystyle\quad\quad(4\pi)^{3}\sum_{m_{i},\ell_{i}}\int\mathcal{D}^{3}\\!P~{}T_{\ell_{1}\ell_{2}\ell_{2};\ell}^{m_{1}m_{2}m_{3};m}T_{\ell_{1}\ell_{2}\ell_{3};\ell^{\prime}}^{*m_{1}m_{2}m_{3};m^{\prime}},$
(168)

where we have suppressed the $k$ dependence in
$T_{\ell_{1}\ell_{2}\ell_{3};\ell_{4}}^{m_{1}m_{2}m_{3};m_{4}}$, defined in
Eq. (64), and the sum is over $\ell_{i}$ and $m_{i}$, with $i=1,2,3$. Without
having to expand the terms with spherical harmonics, we can exploit the fact
that the Universe is statistically isotropic and instead compute,

$\displaystyle\langle\Theta_{\ell m,\rm
inh}^{(1)}\Theta_{\ell^{\prime}m^{\prime},\rm
inh}^{*(1)}\rangle=\frac{\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}}{2\ell+1}\sum_{m^{\prime\prime}}\langle\Theta_{\ell
m^{\prime\prime},\rm inh}^{(1)}\Theta_{\ell m^{\prime\prime},\rm
inh}^{*(1)}\rangle.$ (169)

This enables us to use the machinery we derived in Appendix D and write the
second term as,

$\displaystyle\int\mathcal{D}^{3}\\!P~{}T_{\ell
m}^{(1)}(\bm{k}_{1},\bm{k}_{2},$
$\displaystyle\bm{k}_{3})T_{\ell^{\prime}m^{\prime}}^{*(1)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=$
$\displaystyle(4\pi)^{3}\frac{\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}}{2\ell+1}\\!\\!\sum_{m_{i},\ell_{i}}\\!$
$\displaystyle\int\\!\\!\mathcal{D}^{3}\\!P\\!\left\\{\vphantom{\widetilde{Q}^{2}}A^{2}_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell}\right.$
$\displaystyle+2AB_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})(\mathcal{Q\widetilde{Q}})_{\ell_{1}\ell_{2},\ell_{3}\ell}$
$\displaystyle\left.+B^{2}_{\ell_{1}\ell_{2},\ell_{3}}(k_{1},k_{2},k_{3})(\mathcal{\widetilde{Q}}^{2})_{\ell_{1}\ell_{2},\ell_{3}\ell}\right\\}.$

We apply the same logic to the third term. We then take advantage of the
factorized forms of Eqs. (64)-(66) to write a computationally manageable final
solution as,

$\displaystyle C_{\ell,\rm inh}^{(11)}=f_{\rm
pbh}^{2}\Bigr{[}4\mathfrak{A}_{\ell}+\sum_{\ell_{1}\ell_{2}\ell_{3}}\left(2\mathfrak{B}_{\ell_{1}\ell_{2},\ell_{3};\ell}+4\mathfrak{C}_{\ell_{1},\ell_{2},\ell_{3};\ell}\right)\Bigr{]},$
(170)

where

$\displaystyle\mathfrak{A}_{\ell}\equiv\\!\frac{4\pi}{9}\int_{0}^{\eta_{0}}\\!\\!\\!d\eta\int_{0}^{\eta_{0}}\\!\\!\\!d\eta^{\prime}g(\eta)g(\eta^{\prime})\mathcal{A}_{\ell}(\eta,\eta^{\prime})\beta(\eta)\beta(\eta^{\prime}),$
(171)
$\displaystyle\mathfrak{B}_{\ell_{1}\ell_{2},\ell_{3};\ell}\equiv\frac{(4\pi)^{3}}{2\ell+1}\int
d\eta\int d\eta^{\prime}g(\eta)g(\eta^{\prime})$
$\displaystyle\times\left\\{\mathcal{A}_{\ell_{1}}(\eta,\eta^{\prime})\mathcal{A}_{\ell_{2}}(\eta,\eta^{\prime})J_{\ell_{3}}(\eta,\eta^{\prime})(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell}\right.$
$\displaystyle\quad+2\mathcal{K}_{\ell_{1}}(\eta,\eta^{\prime})\mathcal{K}_{\ell_{2}}(\eta,\eta^{\prime})J_{\ell_{3}}(\eta,\eta^{\prime})(\mathcal{Q\widetilde{Q}})_{\ell_{1}\ell_{2},\ell_{3}\ell}$
$\displaystyle\quad\left.+\mathcal{B}_{\ell_{1}}(\eta,\eta^{\prime})\mathcal{B}_{\ell_{2}}(\eta,\eta^{\prime})J_{\ell_{3}}(\eta,\eta^{\prime})(\mathcal{\widetilde{Q}}^{2})_{\ell_{1}\ell_{2},\ell_{3}\ell}\right\\},$
(172)
$\displaystyle\mathfrak{C}_{\ell_{1},\ell_{2},\ell_{3};\ell}\equiv\frac{(4\pi)^{3}}{2\ell+1}\int
d\eta\int d\eta^{\prime}g(\eta)g(\eta^{\prime})$
$\displaystyle\times\left\\{\widetilde{\mathcal{A}}_{\ell_{1}}(\eta,\eta^{\prime})\widetilde{\mathcal{A}}_{\ell_{2}}(\eta^{\prime},\eta)\mathcal{A}_{\ell_{3}}(\eta,\eta^{\prime})(\mathcal{Q}^{2})_{\ell_{1}\ell_{2}\ell_{3}\ell}\right.$
$\displaystyle\quad\\!\\!\\!+2\widetilde{\mathcal{A}}_{\ell_{1}}(\eta,\eta^{\prime}){\mathcal{K}}_{\ell_{2}}(\eta,\eta^{\prime})\widetilde{\mathcal{B}}_{\ell_{3}}(\eta^{\prime},\eta)(\mathcal{Q\widetilde{Q}})_{\ell_{3}\ell_{2},\ell_{1}\ell}$
$\displaystyle\quad\\!\\!\\!\left.+\widetilde{\mathcal{B}}_{\ell_{1}}(\eta,\eta^{\prime})\widetilde{\mathcal{B}}_{\ell_{2}}(\eta^{\prime},\eta)\mathcal{B}_{\ell_{3}}(\eta,\eta^{\prime})(\mathcal{\widetilde{Q}\widetilde{Q}}^{\rm
S})_{\ell_{1},\ell_{2}\ell_{3},\ell}\right\\}\\!,$ (173)

with

$\displaystyle\mathcal{B}_{\ell}(\eta,\eta^{\prime})$
$\displaystyle\equiv\\!\int\\!\\!D(k)P_{\zeta}(k)\frac{j_{\ell}(\chi k)}{\chi
k}\Delta_{e}(\eta,k)\frac{j_{\ell}(\chi^{\prime}k)}{\chi^{\prime}k}\Delta_{e}(\eta^{\prime},k),$
$\displaystyle\mathcal{K}_{\ell}(\eta,\eta^{\prime})$
$\displaystyle\equiv\\!\int\\!\\!D(k)P_{\zeta}(k)j_{\ell}^{\prime}(\chi
k)\Delta_{e}(\eta,k)\frac{j_{\ell}(\chi^{\prime}k)}{\chi^{\prime}k}\Delta_{e}(\eta^{\prime},k),$
$\displaystyle J_{\ell}(\eta,\eta^{\prime})$
$\displaystyle\equiv\\!\int\\!\\!D(k)P_{\zeta}(k)\mathcal{J}_{\ell}(\eta,k)\mathcal{J}_{\ell}(\eta^{\prime},k),$
$\displaystyle\widetilde{\mathcal{A}}_{\ell}(\eta,\eta^{\prime})$
$\displaystyle\equiv\\!\int\\!\\!D(k)P_{\zeta}(k)j_{\ell}^{\prime}(\chi
k)\Delta_{e}(\eta,k)\mathcal{J}_{\ell}(\eta^{\prime},k),$
$\displaystyle\widetilde{\mathcal{B}}_{\ell}(\eta,\eta^{\prime})$
$\displaystyle\equiv\\!\int\\!\\!D(k)P_{\zeta}(k)\frac{j_{\ell}(\chi k)}{\chi
k}\Delta_{e}(\eta,k)\mathcal{J}_{\ell}(\eta^{\prime},k).$ (174)

We plot the results, Eq. (170) and Eq. (163), in Fig. 9. We see there is
ultimately an order of magnitude difference between the amplitudes of the
inhomogeneous and homogeneous temperature perturbation auto-power-spectrum. We
also see that the correlation between our newly computed inhomogeneous
temperature perturbation and the standard CMB temperature anisotropy is very
poor. Both these facts are likely the culprits behind both the unexpected
sensitivity from the forecast on the trispectrum and the two orders of
magnitude difference in the power spectra amplitudes we observe in Fig. 6.
Additionally, it can be seen that, if it were not for the very small
correlation at large $\ell$, the scale suppression due to photon propagation
would have a much bigger effect in both the inhomogeneous power spectrum and
trispectrum results.

Figure 9:  Top: auto-power-spectrum of the perturbed temperature anisotropy
due to accreting PBHs defined as $\langle\Theta^{(1)}_{\ell
m}\Theta^{*(1)}_{\ell^{\prime}m^{\prime}}\rangle\equiv\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}C_{\ell}^{(11)}$,
normalized by the standard angular power spectrum. Bottom: correlation
coefficients between $\Theta^{(1)}$ and $\Theta^{(0)}$. In both cases we