# Spectral Distortion Anisotropy from Inflation for Primordial Black Holes

David Zegeye, Keisuke Inomata, Wayne Hu  Kavli Institute for Cosmological
Physics, Department of Astronomy & Astrophysics, Enrico Fermi Institute, The
University of Chicago, Chicago, IL 60637, USA

###### Abstract

Single field inflationary models that seek to greatly enhance small scale
power in order to form primordial black holes predict both a squeezed
bispectrum that is enhanced by this small scale power and a potentially
detectable enhancement of CMB spectral distortions. Despite this combination,
spectral distortion anisotropy on CMB scales remains small since the squeezed
bispectrum represents an unobservable modulation of the scale rather than
local amplitude for the short wavelength acoustic power that dissipates and
forms the $\mu$ spectral distortion. The leading order amplitude effect comes
from the local modulation of acoustic dissipation at the beginning of the
$\mu$ epoch at the end of thermalization by a long wavelength mode that is
correlated with CMB anisotropy itself. Compensating factors from the
suppression by the square of the ratio the comoving horizon at thermalization
to the smallest detectable primary CMB scales ($\sim 0.0005$) and maximal
allowed enhancement of $\mu$ ($\sim 5000$) leaves a signal in the $\mu T$
cross spectrum that is still well beyond the capabilities of PIXIE or LiteBIRD
due to sensitivity and resolution while remaining much larger than in single
field slow roll inflation and potentially observable.

##  I Introduction

There has been much recent interest in primordial black holes (PBHs) from an
extremely large enhancement of small scale fluctuations during inflation (e.g.
Kawasaki  _et al._ (2016); Garcia-Bellido and Ruiz Morales (2017); Ezquiaga
_et al._ (2018); Kannike  _et al._ (2017); Germani and Prokopec (2017);
Ballesteros and Taoso (2018); Inomata  _et al._ (2018); Hertzberg and Yamada
(2018); Motohashi and Hu (2017); Cheng  _et al._ (2019); Byrnes  _et al._
(2019); Passaglia  _et al._ (2019); Drees and Xu (2021); Carrilho  _et al._
(2019); Passaglia  _et al._ (2020); Ng and Wu (2021); Inomata  _et al._
(2021a, b)) given their potential to explain dark matter and the binary black
hole mergers detected by LIGO-Virgo Bird  _et al._ (2016); Clesse and García-
Bellido (2017); Sasaki  _et al._ (2016); Kashlinsky (2016); Kashlinsky  _et
al._ (2018); García-Bellido  _et al._ (2021). Moreover primordial fluctuations
on scales much smaller than those probed by the cosmic microwave background
(CMB) and large-scale structure are currently relatively poorly constrained
Josan  _et al._ (2009); Assadullahi and Wands (2010); Inomata and Nakama
(2019); Kawasaki  _et al._ (2021).

Spectral distortions in the CMB are one way to constrain primordial power on
small scales. After the thermalization epoch where photon number changing
processes in the plasma drop out of equilibrium, the energy dissipated in
small scale acoustic waves leaves observable distortions in the frequency
spectrum Sunyaev and Zeldovich (1970); Daly (1991); Hu  _et al._ (1994);
Chluba and Sunyaev (2012); Chluba (2013). The amplitude of these acoustic
waves is itself enhanced if the PBH scale is not much smaller than the
dissipation scale at the end of thermalization. The current bounds for
chemical potential or $\mu$ distortions from COBE-FIRAS already place strong
bounds on such models and rule out PBHs as a significant fraction of the dark
matter between $10^{4}-10^{13}M_{\odot}$ Nakama  _et al._ (2017); Carr  _et
al._ (2020).

While current spectral distortion limits can be greatly improved with future
space-based spectrometers Chluba  _et al._ (2021); Hazumi  _et al._ (2020), an
absolute measurement is limited by contamination from foregrounds and
systematics in addition to instrument sensitivity. Anisotropy in spectral
distortions, if they are large and correlated with CMB anisotropy in
temperature and polarization, provides a promising complementary approach that
uses cross correlation and differential measurement to mitigate these issues
Pajer and Zaldarriaga (2012); Ganc and Komatsu (2012); Pajer and Zaldarriaga
(2013); Chluba  _et al._ (2017); Cabass  _et al._ (2018); Remazeilles and
Chluba (2018); Remazeilles  _et al._ (2021).

These correlations can arise if the amplitude of small scale power is
modulated by long-wavelength fluctuations due to the squeezed bispectrum.
However in single-field slow roll inflation, the squeezed bispectrum obeys the
Maldacena consistency relation Maldacena (2003), where long wavelength
curvature fluctuations modulate the scale and not the amplitude of the small
scale power spectrum. Moreover this scale modulation is unobservable locally
since a coordinate system established using clocks and rulers locally cannot
reference the global coordinate system in which the long-wavelength mode is
embedded Pajer  _et al._ (2013). Since spectral distortions depend only on the
amount of power dissipated and not the globally referenced comoving scale from
which it originates, single-field slow roll inflation does not produce
spectral distortion anisotropy at leading order Pajer and Zaldarriaga (2012);
Cabass  _et al._ (2018).

Recently it has been suggested that single field PBH models may violate the
conditions that make the spectral distortion anisotropy vanish at leading
order in the squeezed bispectrum Özsoy and Tasinato (2021a, b). All such PBH
models must violate the slow roll assumption in order to enhance small scale
power sufficiently rapidly Motohashi and Hu (2017). Proposed PBH models
typically have a period of non-attractor behavior which causes a violation of
the Maldacena consistency relation for long wavelength modes that exit the
horizon sufficiently close to or in the non-attractor phase Namjoo  _et al._
(2013); Martin  _et al._ (2013); Cai  _et al._ (2018); Passaglia  _et al._
(2019).

In this work we show that for scales relevant for correlation with CMB
anisotropy, which are much, much larger than the horizon at the onset of the
non-attractor phase, the Maldacena consistency relation holds in this limit
and spectral distortion anisotropy is suppressed by the square of the ratio
between the long wavelength scale and a characteristic short wavelength scale.
Nevertheless, these suppressed effects can still be much larger than they are
in slow roll inflation Cabass  _et al._ (2018) and appear mainly due to the
modulation of the dissipation scale at the end of thermalization by the long
wavelength mode.

The outline of this paper is as follows. In §II, we review non-attractor
“ultra-slow roll” mechanism for enhancing small scale power during inflation
and discuss the modulation of short-wavelength modes by long-wavelength modes
in the form of the consistency relation. For the relevant long-wavelength
scales we show that the usual Maldacena consistency relation holds and, in
Appendix A, we explicitly verify that once combined with its impact on short-
wavelength acoustic evolution, the zeroth order effect of the long-wavelength
mode is a dilation of scales that is unobservable locally. In §IV and §V, we
calculate the leading order effect of the long wavelength density fluctuation
on the local amplitude of acoustic oscillations and their dissipation at the
end of thermalization respectively. We discuss the implications of this
modulation on spectral distortion anisotropy in §V. In §VII we assess the
prospects for detecting this signal and discuss these results in §VIII.

##  II USR Consistency Relation

A common aspect of many inflationary models that enhance the small scale
curvature power spectrum by the orders of magnitude that would be required to
later form PBHs in the radiation dominated epoch is a transient period of so-
called “ultra slow-roll” (USR) Kinney (2005). During USR, the inflaton $\phi$
rolls faster than can be sustained by the slope of its potential $V(\phi)$,

$\left|\frac{d\phi}{dN}\right|\gg\left|\frac{V^{\prime}(\phi)}{H^{2}}\right|,$
(1)

where $H=d\ln a/dt=dN/dt$ is the Hubble parameter and $N$ is the efold. This
excess kinetic energy can arise for example from a very flat potential around
an inflection point (e.g. Garcia-Bellido and Ruiz Morales (2017)) or from a
sudden increase in kinetic energy due to a downward feature in the potential
Inomata  _et al._ (2021a). In either case, the excess kinetic energy of the
field then redshifts away as $(d\phi/dN)^{2}\propto a^{-6}$ and the curvature
fluctuation in unitary gauge grows as $\zeta=-\delta\phi/(d\phi/dN)\propto
a^{3}$ once it crosses the horizon.

Modes that were already outside the horizon at the beginning of the USR phase
also experience growth to the extent that they have not completely frozen out.
Following the treatment in Ref. Passaglia  _et al._ (2019), during the
preceding slow roll (SR) phase superhorizon curvature fluctuations evolve to a
constant as $d\ln\zeta/dN\approx-(k/aH)^{2}$. Once USR begins
$d\zeta/dN\propto a^{3}$ and so the value of $(k/aH)_{\rm USR}^{2}$, where the
subscript denotes the beginning of the USR epoch, determines the growth of
these modes. Given the opposite sign of $d\ln\zeta/dN$ during SR, there is a
mode for which this growth is just sufficient to overcome the initially
constant SR piece and cause a near zero crossing in the curvature power
spectrum at $k=k_{\rm dip}$ where

$\left(\frac{k_{\rm dip}}{aH}\right)^{2}_{\rm USR}\sim e^{-3\Delta N_{\rm
USR}},$ (2)

and $\Delta N_{\rm USR}$ is the number of efolds in the USR phase. For example
for the $\sim 10^{7}$ enhancement of power typical to PBH models $\Delta
N_{\rm USR}\approx\ln(10^{7})/6$ and so $k_{\rm dip}\sim 10^{-2}(aH)_{\rm
USR}$. Note that these scalings also imply that the enhancement in power for
$k_{\rm dip}<k\lesssim(aH)_{\rm USR}$ scales as $k^{4}$ Byrnes  _et al._
(2019). Thus the curvature power spectrum for $k\lesssim(aH)_{\rm USR}$ can be
parameterized as

$\Delta^{2}_{\zeta}(k)\equiv\frac{k^{3}P_{\zeta}(k)}{2\pi^{2}}=\Delta^{2}_{{\rm
SR}}(k)\left[1-\left(\frac{k}{k_{\rm dip}}\right)^{2}\right]^{2},$ (3)

where $\Delta^{2}_{\rm SR}$ is the SR power spectrum which we take for
simplicity to be a pure power law

$\Delta^{2}_{\rm SR}(k)=A_{s}\left(\frac{k}{0.05\,{\rm
Mpc}^{-1}}\right)^{n_{s}-1},$ (4)

with $A_{s}=2.1\times 10^{-9}$ and $n_{s}=0.965$ consistent with CMB
measurements Aghanim  _et al._ (2020). Note that this form can be derived more
rigorously in specific models and it holds to leading order in the downward
step PBH model (see Inomata  _et al._ (2021b) Eq. (39)). For
$k\gtrsim(aH)_{\rm USR}$ the result becomes dependent on the specific PBH
model but for the modes relevant for spectral distortions that satisfy current
observational constraints Eq. (3) suffices.

In Fig. 1, we show the power spectrum (3) for $k_{\rm dip}=681$ Mpc-1, which
we shall see below in Fig. 2 is the largest scale allowed by current
constraints on spectral distortions which arise from the power near the
dissipation scale at end of thermalization $k_{D}(z_{\rm th})\approx 10^{4}$
Mpc-1. Notice the large hierarchy of scales between these values and the
smallest scale accessible to measurements of the primary CMB anisotropy
$k_{L}\sim 0.1$ Mpc-1 due to the dissipation scale at recombination.

Figure 1:  The PBH curvature power spectrum $\Delta^{2}_{\zeta}(k)$ for a USR
enhancement starting at $k_{\rm dip}=681$ Mpc-1 which provides the largest
spectral distortion still allowed (see Fig. 2). Vertical lines show
$k_{L}=0.1$ Mpc-1 (of order the largest $k$ mode probed by the CMB),
$(aH)_{\rm th}$ the comoving size of the horizon at the end of thermalization,
$k_{\rm dip}$, and the dissipation wavenumber at the end of thermalization
$k_{D}(z_{\rm th})$. The slow roll power law extrapolation $\Delta_{\rm
SR}^{2}$ is also shown for reference.

The wavenumber $k_{\rm dip}$ separates two very different regimes for the
impact of long wavelength fluctuations on much shorter wavelength power or
equivalently the squeezed bispectrum. For a long wavelength $k_{L}\ll k_{\rm
dip}$ the squeezed bispectrum obeys the usual Maldacena consistency relation
Maldacena (2003)

$\displaystyle\lim_{k_{L}/k_{S}\rightarrow 0}B_{\zeta}(k_{L},k_{S},k_{S})$
$\displaystyle=$
$\displaystyle-\left[\frac{d\ln\Delta^{2}_{\zeta}(k_{S})}{d\ln k_{S}}+{\cal
O}\left(\frac{k_{L}}{k_{\rm dip}}\right)^{2}\right]$ (5) $\displaystyle\times
P_{\zeta}(k_{L})P_{\zeta}(k_{S}),$

which can be demonstrated by explicit calculations using the in-in or $\delta
N$ formalisms Passaglia  _et al._ (2019). The interpretation of this relation
is that the scale of the short wavelength modes are dilated by the nearly
constant long wavelength curvature perturbation, which acts as a spatial
fluctuation in the local scale factor

$\displaystyle\Delta_{\zeta}^{2}(k_{S},{\bm{y}})$ $\displaystyle\approx$
$\displaystyle\Delta_{\zeta}^{2}(k_{S}(1-\zeta_{L}({\bm{y}}))$ (6)
$\displaystyle\approx$
$\displaystyle\Delta_{\zeta}^{2}(k_{S})\left(1-\frac{d\ln\Delta^{2}_{\zeta}(k_{S})}{d\ln
k_{S}}\zeta_{L}({\bm{y}})\right),$

where local spatial variations are denoted by the comoving coordinate
${\bm{y}}$. Notice that the amplitude of the power spectrum does not actually
change locally, just the comoving scale that the power is associated with.
Physically, the long wavelength mode just changes the efold at which the short
wavelength mode enters into USR, not any of the dynamics due to the USR
growth. As shown in Ref. Cabass  _et al._ (2018) this change in scale cannot
modulate spectral distortions, which depend only on the amplitude of the power
dissipated, not the scale. In the Appendix we explicitly verify that a
constant $\zeta_{L}$ generates a locally unobservable dilation for both the
primordial non-Gaussianity and the subsequent dynamics of short wavelength
modes to second order in perturbations.

For $k_{L}\gtrsim k_{\rm dip}$, the leading effect of the long wavelength mode
is not simply a dilation of scales and hence the consistency relation (5) no
longer holds Namjoo  _et al._ (2013); Martin  _et al._ (2013). Physically, the
inflaton field fluctuation of the long wavelength mode changes the number of
efolds that the short wavelength mode experiences USR since there is no longer
an attractor solution that makes it equivalent to a shift along the background
phase space trajectory Passaglia  _et al._ (2019). For example, a constant
backwards fluctuation of the field means that locally the short wavelength
modes see more efolds of USR growth and so the actual amplitude of the small
scale power spectrum is modulated by the long wavelength mode. This then would
produce a spatial modulation in spectral distortions.

We shall see that for modes relevant for correlation with CMB anisotropy and
models that satisfy spectral distortion constraints, $k_{L}\ll k_{\rm dip}$.
Even in this regime, there is always some dynamical effect of short wavelength
growth in the long wavelength field to the extent that $\zeta_{L}\neq$ const.
As we have seen, here the evolving part of $\zeta_{L}$ during USR is
suppressed by $k_{L}^{2}/k_{\rm dip}^{2}$. This expectation is consistent with
an explicit calculation of the leading order correction the Maldacena
consistency relation given in Eq. (5) in an inflection point PBH model
Passaglia  _et al._ (2019). However since again the relevant $k_{L}\ll k_{\rm
dip}$ this correction produces a very small $\mu T$ correlation (see Fig. 1
and §VI).

##  III Local Expansion

After inflation, there are $(k_{L}/aH)^{2}$ suppressed modulations of short
wavelength physics by the long wavelength curvature fluctuations $\zeta_{L}$,
analogous to the corrections to the primordial bispectrum from $k_{\rm
dip}\propto(aH)_{\rm USR}$ in the previous section. To the local observer, the
long wavelength mode appears as a change to the background cosmology induced
by the density perturbation that it carries. Since synchronous observers are
freely falling test particles, we can absorb these synchronous gauge adiabatic
density fluctuations $\delta_{L}$ Hu  _et al._ (2016); Hu and Joyce (2017),

$\delta_{L}=\frac{1}{3}\left(\frac{k}{aH}\right)^{2}\zeta_{L}\propto a^{2}$
(7)

during radiation domination into a new background or “separate universe”.

Following the construction of the separate universe for the late time growth
of structure Sirko (2005); Li  _et al._ (2014), we can define a local
background density $\bar{\rho}_{L}$ as

$\bar{\rho}(1+\delta_{L})=\bar{\rho}_{L},$ (8)

which implies that the local observer sees a scale factor $a_{L}$ that is
related to the global scale factor $a$ at equal times as

$a_{L}\approx a(1-\delta_{L}/4),$ (9)

where we chosen the normalization such that $a_{L}\rightarrow a$ at early
times. Notice that this normalization removes the dilation effect of
$\zeta_{L}$ in Eq. (6) by measuring scales locally so that they coincide when
$\delta_{L}\ll 1$ at the end of inflation.

By virtue of this normalization and conservation of particles, at the same
numerical value of the scale factors (or efolds) from the end of inflation in
the local and global universe, all particle number and energy densities are
the same. However the scale factors do not coincide at the same time in the
local and global universe.

To extract the cosmological parameters of the separate universe, we can
express the expansion rate as a function of the local scale factor. With the
definition $H_{L}\equiv d\ln a_{L}/dt$ and the radiation dominated growth of
$\delta_{L}\propto a^{2}$ we have at equal times

$H_{L}^{2}\approx H^{2}(1-\delta_{L}).$ (10)

In the global universe let us define a reference epoch $a_{r}$ in the
radiation dominated limit where $H^{2}(a)=H_{r}^{2}(a_{r}/a)^{4}$ so that we
can express the local expansion rate in terms of the local scale factor as

$H_{L}^{2}=H_{r}^{2}\left(\frac{a_{r}}{a}\right)^{4}(1-\delta_{L})\approx
H_{r}^{2}\left(\frac{a_{r}}{a_{L}}\right)^{4}(1-2\delta_{L}).$ (11)

In the local universe this takes the form of the Friedmann equation with
$\delta_{L}\propto a^{2}\approx a_{L}^{2}$ playing the role of spatial
curvature to linear order in $\delta_{L}$, specifically

$H_{0}^{2}\Omega_{KL}=-2{\delta_{L}}(a_{r})H_{r}^{2}a_{r}^{2}$ (12)

so that $\delta_{L}(a_{r})$ at the reference epoch defines the comoving
curvature scale in units of the comoving Hubble scale at that epoch. Finally
it is useful to express the local conformal time as

$\eta_{L}=\int\frac{d\ln
a_{L}}{a_{L}H_{L}}\approx\frac{a_{L}}{a_{r}}\frac{1}{a_{r}H_{r}}\left(1+\frac{\delta_{L}}{3}\right).$
(13)

These relations now determine the modulation of all short wavelength
observables by a long wavelength curvature fluctuation $\zeta_{L}$ which to
leading order scales as $(k_{L}/aH)^{2}\zeta_{L}$. The specific size of the
modulation will then depend on the epoch at which it influences the short
wavelength observable the most. We shall see that for $\mu$ distortions and
PBH models, this is the end of the thermalization epoch.

##  IV Modulated Acoustic Power

The first step in understanding the local modulation of spectral distortions
is to determine the impact of the long wavelength curvature perturbation
$\zeta_{L}$ on the amplitude of short wavelength acoustic oscillations in the
CMB during radiation domination. These oscillations then dissipate via
diffusion damping leaving a spectral distortion after the thermalization
epoch. We shall see that since acoustic oscillations are generated at horizon
crossing of the short wavelength mode $k_{S}$, the impact of the long
wavelength mode occurs at horizon crossing of the short wavelength mode
$aH=k_{S}$ and therefore scales as $(k_{L}/k_{S})^{2}\zeta_{L}$. The impact of
CMB scale wavenumbers $k_{L}$ is therefore highly suppressed for the modes
that contribute to spectral distortions (see Fig. 1).

We can analytically understand this scaling in the simple case where the
photons dominate the radiation density, i.e. neglecting the effect of
neutrinos which only change the numerical factors and not the overall scaling.
In this case we can solve the perturbation equations in terms of the
continuity and Euler equations for the photon fluid under self gravity in the
photon-dominated local universe as (see Hu and White (1996) Eq. (10))

$\displaystyle\Delta_{\gamma}^{\prime}-\frac{y^{\prime}}{y}\Delta_{\gamma}$
$\displaystyle=$
$\displaystyle-\frac{4}{\sqrt{3}}\frac{1-{y^{\prime\prime}}/{y}+2({y^{\prime}}/{y})^{2}}{1+{6}/{(f_{K}y^{2})}}V_{\gamma},$
$\displaystyle V_{\gamma}^{\prime}+\frac{y^{\prime}}{y}V_{\gamma}$
$\displaystyle=$
$\displaystyle\left(1-\frac{6}{f_{K}y^{2}}\right)\frac{\sqrt{3}}{4}\Delta_{\gamma},$
(14)

where $\Delta_{\gamma}$ is the comoving gauge photon density perturbation,
$y=(\Omega_{\gamma}H_{0}^{2})^{-1/2}a_{L}k_{S}$, ${}^{\prime}=d/dx$ with
$x=k_{S}\eta_{L}/\sqrt{3}$, and $f_{K}=1+3\Omega_{KL}H_{0}^{2}/k_{S}^{2}$.
Here we have again made use of the fact that at the same value for the scale
factor, the physical density of the photons is the same in the local and
global universe. The synchronous and comoving gauge differ in their density
perturbations outside the horizon due to radiation pressure so
$\lim_{x\rightarrow 0}\Delta_{\gamma}\neq\delta_{L}$ but they do approach each
other for $x\gg 1$ in the regime relevant for spectral distortions.

Without the curvature perturbation induced by the long wavelength mode,
$x=y/\sqrt{3}$ during photon domination and the solution is analytic

$\lim_{y\rightarrow{\sqrt{3}}x}\Delta_{\gamma}=4\left(\frac{\sin x}{x}-\cos
x\right)\zeta_{S},$ (15)

from which we can read off the usual transfer function111This transfer
function is reduced to $-4/(1+4R_{\nu}/15)\sim-3.61$ when neutrinos with
$R_{\nu}=\rho_{\nu}/(\rho_{\nu}+\rho_{\gamma})$ are included Hu and Sugiyama
(1996) for acoustic oscillations $-4\cos x$ at $x\gg 1$

We can now solve for the leading order correction from the small curvature
induced by the long wavelength mode. Since we neglect neutrinos
$H_{r}^{2}=\Omega_{\gamma}H_{0}^{2}a_{r}^{-4}$ and we have

$y\approx{\sqrt{3}}x\left(1-\frac{\delta_{L}}{3}\right)\approx{\sqrt{3}}x\left(1-\frac{x^{2}}{{3}}\alpha\right),$
(16)

so that Eq. (14) becomes to leading order:

$\Delta_{\gamma}^{\prime\prime}+\left(1-\frac{2}{x^{2}}-\frac{4}{3}\alpha\right)\Delta_{\gamma}=0$
(17)

with

$\alpha\equiv\left(\frac{k_{L}}{k_{S}}\right)^{2}\zeta_{L}$ (18)

and $f_{K}=1-2\alpha$ constant in time.

Notice that at $x\gg 1$ this takes the form of an oscillator equation with a
perturbed constant frequency

$\tilde{x}=\left(1-\frac{2}{3}\alpha\right)x.$ (19)

We can now iterate to solve Eq. (17) to first order in $\alpha$ using the
zeroth order solution (15) to determine $\alpha\Delta_{\gamma}$ as an external
source,

$\Delta_{\gamma}\approx
4\zeta_{S}\left(1+\frac{4}{3}\alpha\right)\left(\frac{\sin\tilde{x}}{\tilde{x}}-\cos\tilde{x}\right).$
(20)

We can explicitly verify that this form solves (17) and satisfies
$\Delta_{\gamma}\approx 4\zeta_{S}x^{2}/3$ in $x\ll 1$ to linear order in
$\alpha$.

Therefore there is an $\alpha=(k_{L}/k_{S})^{2}\zeta_{L}$ change in the
amplitude and frequency of the acoustic wave. Since $k_{L}/k_{S}\ll 1$, this
is a large suppression factor and produces a negligible change in the local
spectral distortion once the acoustic modes have dissipated. We can therefore
hereafter assume that the power in acoustic modes at $k_{S}$ is effectively
the same in the global and local universe at the same value of the scale
factors.

##  V Modulated Thermalization

We can now compute the local dissipation of energy from the acoustic waves
into $\mu$ spectral distortions in the separate universe. Changes in
thermalization due to the local background induce larger local variations in
$\mu$ for PBH models, than the primordial effect from $k_{\rm dip}$ or the
acoustic growth, since we shall see they scale with the comoving horizon size
at the end of thermalization $(k_{L}/aH)_{\rm th}^{2}\zeta_{L}$.

Let us now see how thermalization is altered by the long-wavelength modulation
in the local universe. The thermalization rate for the joint action of double
Compton scattering $e^{-}+\gamma\leftrightarrow e^{-}+2\gamma$, which changes
photon number, and Compton scattering $e^{-}+\gamma\leftrightarrow
e^{-}+\gamma$, which redistributes energy, to establish a blackbody scales as
(see e.g. Danese and de Zotti (1982); Hu and Silk (1993))

$\Gamma_{\rm th}\propto T^{3/2}n_{e}\propto a_{L}^{-9/2},$ (21)

where $n_{e}$ is the free electron density and $T$ is the plasma temperature
and recall that at the same numerical value of the scale factors all particle
densities are the same in the local and global universe.

Let us define the thermalization time by the condition

$\int dt\Gamma_{\rm th}=\int_{\ln a_{\rm th}}^{0}d\ln a_{L}\frac{\Gamma_{\rm
th}}{H_{L}}=1.$ (22)

Using Eq. (11), the change in the efold of thermalization is given by

$\Delta\ln a_{\rm th}={2}\delta_{L}(a_{{\rm th}})$ (23)

such that in an overdensity, thermalization continues to a later efold due to
a slower expansion rate.

The energy in acoustic waves is dissipated when the photon diffusion length
crosses the wavelength. If this occurs after the thermalization epoch but
before the Compton scattering becomes inefficient at redistributing energy,
then this energy is transferred into a $\mu$ spectral distortion (see Eq. (30)
below). The diffusion wavenumber in the radiation epoch is given by

$k_{D}^{-2}\propto\int\frac{d\ln a}{aH}\frac{1}{n_{e}\sigma_{T}a},$ (24)

where $\sigma_{T}$ is the Thomson cross section (see Eq. (28) below for the
general expression). Again let us make use of the fact that the particle
densities are the same in the local and global universe at the same efold. We
then get the change in the diffusion wavenumber as

$\displaystyle\Delta\ln k_{D}(a_{\rm th})$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\Delta\ln a_{\rm th}-\frac{3}{10}\delta_{L}(a_{{\rm
th}})$ $\displaystyle=$ $\displaystyle-\frac{33}{10}\delta_{L}(a_{{\rm
th}})=-\frac{11}{10}\left(\frac{k_{L}}{aH}\right)_{\rm th}^{2}\zeta_{L},$

where recall that in the global universe $k_{D}\propto a^{-3/2}$. Therefore in
an overdensity the diffusion wavenumber decreases. The net result is that the
maximal wavenumber that dissipates into $\mu$ distortions given the same small
scale power spectrum in local coordinates as global coordinates is modulated
by $(k_{L}/aH)_{\rm th}^{2}\zeta_{L}$. Since the power spectrum is strongly
blue tilted, it is this modulation that changes the local value of $\mu$. This
is in contrast to single field slow roll inflation where the larger horizon
size at the end of the $\mu$ epoch makes the modulation at that time the
dominant effect Cabass  _et al._ (2018).

Figure 2:  The average spectral distortion $\bar{\mu}$ for various PBH
enhancement scales $k_{\rm dip}$ (see Eq. (3)). Values of $k_{\rm dip}<k_{\rm
min}$ are excluded by the FIRAS bound of $\bar{\mu}\leq 9\times 10^{-5}$
whereas allowed values that lead to enhancement over the slow roll prediction
of $\bar{\mu}\approx 2\times 10^{-8}$ follow $\bar{\mu}\propto k_{\rm
dip}^{-4}$.

##  VI Spectral Distortion Anisotropy

With the local change in the thermalization and dissipation scales due to the
long wavelength curvature perturbation as calculated in the previous section,
we can now determine how it modulates the $\mu$ spectral distortion. Because
long wavelength curvature perturbations also generate CMB anisotropy, this
leads to a potentially observable $\mu T$ correlation in PBH models.

First we calculate the average spectral distortion in the global universe with
the PBH power spectrum of Eq. (3). Following Cabass  _et al._ (2018); Chluba
_et al._ (2017),

$\bar{\mu}=\int d\ln k\Delta_{\zeta}^{2}(k)W(k),$ (26)

where to good approximation

$\displaystyle W(k)$ $\displaystyle\approx$
$\displaystyle-4.54k^{2}\int_{0}^{\infty}{d}z\frac{dk_{D}^{-2}}{dz}\mathcal{J}_{\mu}(z)\mathrm{e}^{-{2k^{2}}/{k_{D}^{2}(z)}}.$
(27)

Here the diffusion wavenumber $k_{D}$ is given by

$k_{D}^{-2}(z)=\frac{1}{6}\int_{z}^{\infty}\frac{dz}{H}\frac{1}{n_{e}\sigma_{T}a}\frac{R^{2}+16(1+R)/15}{(1+R)^{2}},$
(28)

where $R=3\rho_{b}/4\rho_{\gamma}$. In the radiation dominated epoch and with
the best fit $\Lambda$CDM parameters

$k_{D}\approx 4.05\times 10^{-6}(1+z)^{3/2}\,{\rm Mpc}^{-1}.$ (29)

The Green function $\mathcal{J}_{\mu}$ for $\mu$ distortions is well
approximated by Chluba (2013)

$\displaystyle\mathcal{J}_{\mu}\left(z\right)$ $\displaystyle\approx$
$\displaystyle\left[1-\exp\left(-\left[\frac{1+z}{5.8\times
10^{4}}\right]^{1.88}\right)\right]\Theta(z-z_{\rm rec})$ (30)
$\displaystyle\times e^{-(z/z_{{\rm th}})^{5/2}},$

where we have included a $\Theta$ step function at recombination $z_{\rm rec}$
since below this redshift there are no acoustic waves to dissipate. Here
$z_{\rm th}\approx 2\times 10^{6}$ is the thermalization redshift in the
global universe and the quantity in brackets of Eq. (30) determines the
gradual transition from the $\mu$ epoch to the Compton $y$ epoch. Notice that
$\bar{\mu}$ receives contributions from the initial inflationary power
spectrum at wavenumber $k$ mainly when it crosses $k_{D}$ but these
contributions are sharply cut off by the thermalization process above
$k_{D}(z_{\rm th})$.

In Fig. 2 we show $\bar{\mu}$ as a function of $k_{\rm dip}$ in Eq. (3). The
COBE-FIRAS constraint $\bar{\mu}<9\times 10^{-5}$ (95% CL) Fixsen  _et al._
(1996) places a limit of

$k_{\rm dip}>k_{\rm min}\approx 681\,{\rm Mpc}^{-1}$ (31)

and

$\bar{\mu}\approx 9\times 10^{-5}\left(\frac{k_{\rm min}}{k_{\rm
dip}}\right)^{4}$ (32)

for models with $k_{\rm min}<k_{\rm dip}\lesssim 5000$ Mpc-1 whereas for even
larger $k_{\rm dip}$, $\bar{\mu}$ asymptotes to its slow-roll value of
$\bar{\mu}\approx 2\times 10^{-8}$ because enhanced scales have already
dissipated by the thermalization epoch. Notice that even this smallest allowed
value of $k_{\text{dip}}$ is much greater than the comoving horizon at the end
of thermalization $(aH)_{\rm th}\approx 4.3$ Mpc-1.

The modulation of $k_{D}(z_{\rm th})$ due to the long wavelength mode
therefore modulates the local value of $\mu$ from its background value
$\bar{\mu}$. Given that the PBH power spectrum rises as approximately
$\Delta_{\zeta}^{2}\propto k^{{4-(1-n_{s})}}$ in the region which can enhance
$\mu$,222For consistency with Eq. (3) we have retained $1-n_{s}$ here but note
that any actual ${\cal O}(1-n_{s})$ correction to the local slope will depend
on the details of the model. we get the fractional change in the local value
of $\mu$ as

$\delta\ln\mu=[4-(1-n_{s})]\Delta\ln k_{D}\equiv b_{\rm
th}\left(\frac{k_{L}}{aH}\right)_{\rm th}^{2}\zeta_{L},$ (33)

where from Eq. (V) $b_{\rm th}\approx-{22}/{5}+11(1-n_{s})/10$. Note that this
approximation can in the future be improved by numerically recalibrating the
Green function $\mathcal{J}_{\mu}$ in the separate universe so we parameterize
the result as a “bias” factor $b_{\rm th}$, similar to the slow-roll
calculation of Ref. Cabass  _et al._ (2018) but with respect to the horizon
scale at the end of thermalization rather than the end of the $\mu$ epoch. For
example if we assume that the functional form of $\mathcal{J}_{\mu}$ remains
the same and only $z_{\rm th}$ changes according to Eq. (23), then $b_{\rm
th}\approx-4.1$ for $k_{\rm dip}=k_{\rm min}$.

Given the level of precision in these estimates, we simply take

$b_{\rm th}=-22/5$ (34)

as our fiducial bias.

Figure 3:  Maximum $\mu$ white noise level $C_{\ell}^{\mu\mu}=A_{\rm noise}$
that allows for a $S/N=1$ measurement of the maximal $\bar{\mu}=9\times
10^{-5}$ and fiducial $b_{\rm th}=-22/5$ PBH signal as a function of largest
multipole measured $\ell_{\rm max}$. $A_{\rm noise}$ saturates at $10^{-23}$
around $\ell_{\rm max}\sim 10^{3}$ and other models can be scaled as $A_{\rm
noise}\propto b_{\rm th}^{2}\bar{\mu}^{2}$.

Following Ref. Cabass  _et al._ (2018), we can characterize the long
wavelength correlation between CMB temperature anisotropy and $\mu$ anisotropy
with the angular cross power spectrum

$C_{\ell}^{\mu T}\equiv\bar{\mu}b_{\rm th}C_{\ell}^{\mu T,b_{\rm th}},$ (35)

where

$\displaystyle C_{\ell}^{\mu T,b_{\rm th}}=\frac{4\pi}{(aH)_{\rm
th}^{2}}\int_{0}^{\infty}{d}kk\Delta^{2}_{\zeta}(k)\Delta_{\ell}^{\mu}(k)\Delta_{\ell}^{T}(k).$
(36)

Here $\Delta_{\ell}^{T}$ is the CMB temperature transfer function which we
take from CAMB333 https://camb.info and

$\Delta^{\mu}_{\ell}(k)=e^{-k^{2}/q_{\mu,D}^{2}}j_{\ell}(\eta_{0}-\eta_{\rm
rec}),$ (37)

where $q_{\mu,D}\approx 0.084$ Mpc-1 is the dissipation scale of $\mu$
inhomogeneities at recombination Pajer and Zaldarriaga (2013).444This
approximation can also be refined in the future with the Green function for
spectral spatial anisotropy (Chluba, private communication). Notice that this
damping factor is comparable to that of the temperature transfer function at
$k_{D}(z_{\rm rec})$ and in combination they limit the integral in Eq. (35) to
the long wavelength $k_{L}$ values of the CMB. It is straightforward to
generalize this result for the cross correlation with CMB polarization with
the polarization transfer function which we leave that to a future work.

Finally notice that the primordial deviation from the dilation or consistency
relation and the separate universe growth of acoustic oscillations takes the
same form as Eq. (35) but are suppressed by factors of $(aH/k_{\rm
dip})^{2}_{\rm th}\lesssim 4\times 10^{-5}$ and $(aH/k_{D})^{2}_{\rm
th}\approx 10^{-7}$ respectively and hence provide a negligible correction to
$C_{\ell}^{\mu T}$ as calculated in Eq. (35) from the thermalization bias.

##  VII Signal-to-Noise

The $\mu T$ cross power spectrum as calculated in the previous section is
enhanced in PBH models by the $k^{4}$ rise in the small scale curvature power
spectrum as long as $k_{\rm dip}\ll k_{D}(z_{\rm th})$ but also suppressed by
the smallness of the density perturbation associated with $\zeta_{L}$ at the
thermalization epoch $(k_{L}/aH)^{2}_{\rm th}\zeta_{L}$. Consequently unlike
$\mu T$ correlations from nearly scale invariant perturbations in multifield
inflation Pajer and Zaldarriaga (2012); Ganc and Komatsu (2012); Pajer and
Zaldarriaga (2013); Chluba  _et al._ (2017); Cabass  _et al._ (2018);
Remazeilles  _et al._ (2021), the signal-to-noise is dominated by the smallest
angular scales at which the correlation can be detected, namely the damping
scale of primary CMB anisotropy $k_{L}\sim 0.1$ Mpc-1.

To estimate the signal-to-noise in $\mu T$, we take the Gaussian approximation
for $\mu$ anisotropy

$\left(\frac{S}{N}\right)^{2}=\sum_{\ell=2}^{\ell_{\rm
max}}(2\ell+1)\frac{(C_{\ell}^{\mu T})^{2}}{(C_{\ell}^{\mu
T})^{2}+C_{\ell}^{\mu\mu}C_{\ell}^{TT}},$ (38)

where the $C_{\ell}^{\mu\mu}$ and $C_{\ell}^{TT}$ terms in the denominator
include both the sample variance of the signal and any instrumental noise from
their measurement.

Given the smallness of the $\mu$ anisotropy, $C_{\ell}^{\mu\mu}$ will be noise
dominated for the forseeable future.555The ultimate limit comes from the
residual fluctuations from the averaging across CMB scales of the patches that
dissipate Pajer and Zaldarriaga (2012, 2013); Cabass  _et al._ (2018) but note
that in the PBH context the patches that contribute most to $\bar{\mu}$ are
much smaller in size than the SR contributions at the end of the $\mu$ epoch
and do not appreciably enhance the noise. To assess the maximal white noise
level $C_{\ell}^{\mu\mu}=A_{\rm noise}$ at which the signal is barely
detectable at $S/N=1$, we can drop the $C_{\ell}^{\mu T}$ sample variance term
in the denominator of Eq. (38) and solve for white noise term

$\displaystyle{A_{\rm noise}}$ $\displaystyle=$
$\displaystyle\sum_{\ell=2}^{\ell_{\rm max}}(2\ell+1)\frac{(C_{\ell}^{\mu
T})^{2}}{C_{\ell}^{TT}}$ $\displaystyle=$ $\displaystyle{b_{\rm
th}^{2}\bar{\mu}^{2}}\sum_{\ell=2}^{\ell_{\rm
max}}(2\ell+1)\frac{(C_{\ell}^{\mu T,b_{\rm th}})^{2}}{C_{\ell}^{TT}}.$

To estimate the maximal noise level for detection, we assume that the $TT$
measurement is cosmic variance limited to $\ell_{\rm max}$. In Fig. 3, we show
this maximal noise amplitude as a function of $\ell_{\rm max}$ for the maximal
signal $\bar{\mu}=9\times 10^{-5}$ and $b_{\rm th}=-22/5$. Notice that the
result saturates at around $\ell_{\rm max}\approx 10^{3}$ at a level of $\sim
10^{-23}$ since both $T$ and $\mu$ anisotropies are damped by diffusion. To
constrain other PBH models, this result can be scaled as $A_{\rm
noise}\propto{b_{\rm th}^{2}\langle\mu\rangle^{2}}$ using Eq. (32).

As this represents the maximal signal, detecting this effect will be
challenging experimentally. First, to optimize detection, an experiment would
need at least several arcminute scale resolution since

$C^{\mu\mu}_{\ell}=A_{\rm noise}e^{\frac{\ell^{2}\theta_{b}^{2}}{8\ln 2}},$
(40)

where $\theta_{b}$ is the full width half max of the beam in radians.

Second, the required noise level $A_{\rm noise}<10^{-23}$ is quite stringent
even ignoring foregrounds and systematics. For reference, a PIXIE-like mission
Kogut  _et al._ (2011); Kogut and Fixsen (2020) which aims to measure
$\overline{\mu}$ in single field slow roll with $\theta_{b}\approx
1.6^{\circ}$, will at best achieve $A_{\rm noise}\sim 10^{-15}$ which would
suffice to constrain a squeezed bispectrum signal from multifield inflation to
$\left|f_{\mathrm{NL}}\right|\lesssim 3000$ at 68% CL Pajer and Zaldarriaga
(2012); Ganc and Komatsu (2012). To constrain the maximal PBH signal, one
would require the equivalent sensitivity to $|f_{\rm NL}|<0.3$ but with a much
higher angular resolution. LiteBIRD can potentially achieve a detector
sensitivity of $A_{\rm noise}\approx 10^{-18}$ but still with only $\sim
0.5^{\circ}$ resolution and subject to foreground contamination Remazeilles
_et al._ (2021).

##  VIII Discussion

We have shown that in spite of the large squeezed bispectrum due to the
enhancement of small scale power in single field inflationary PBH models, the
spectral distortion anisotropy is highly suppressed since for scales relevant
to CMB cross correlation, it represents an unobservable modulation of global
scales rather than of the local amplitude of the short wavelength modes.
Nonetheless the $\mu$ anisotropy can be larger than in single field slow roll
inflation because of the large enhancement of small scale power itself. The
largest effect comes from the local modulation of the expansion rate at a
given locally measured efold from the end of inflation and hence the end of
the thermalization epoch. This modulation provides a spatial variation in the
amount of power in acoustic waves dissipated near thermalization that causes a
$\mu$ distortion that is correlated with CMB anisotropy itself.

This leading order correlation is enhanced by the $k^{4}$ rise in small scale
power in PBH models but suppressed by the square of the ratio of the comoving
horizon at end of thermalization to observable CMB scales, $(k_{L}/aH)_{\rm
th}^{2}\lesssim 0.0005$. On the other hand, the enhancement of the average
$\mu$ itself can be up to $\sim 5000$ and still satisfy current COBE-FIRAS
bounds. These compensating factors leave the signal potentially observable but
still well beyond the capability of proposed space-based instruments like
PIXIE and LiteBIRD. Moreover to detect the correlation at the smallest
observable CMB scales, where the signal peaks, would require a telescope with
several arcminute scale resolution or better.

These properties of the PBH $\mu T$ correlation suggest that in the future,
ground based instruments may provide a competitive path forward, given the
rapid advance in the scale of detectors deployed on large telescopes into the
CMB-S4 Abazajian  _et al._ (2016) era and beyond. Unlike the absolute
measurement of $\bar{\mu}$, systematics and foregrounds can also be mitigated
by differential measurements and cross correlation Ganc and Komatsu (2012).
Furthermore the $\mu T$ correlation can be supplemented by polarization cross
correlation Ravenni  _et al._ (2017); Remazeilles  _et al._ (2021). We leave
these studies to a future work.

###### Acknowledgements.

We thank Giovanni Cabass, Jens Chluba, Tom Crawford, Austin Joyce, Hayden Lee,
Enrico Pajer, Sam Passaglia, and LianTao Wang for useful discussions. D.Z. was
supported by the National Science Foundation Graduate Research Fellowship
Program under Grant No. DGE1746045. K.I. was supported by the Kavli Institute
for Cosmological Physics at the University of Chicago through an endowment
from the Kavli Foundation and its founder Fred Kavli. W.H. was supported by
U.S. Dept. of Energy contract DE-FG02-13ER41958 and the Simons Foundation.
Portions of this work was performed at the Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-1607611.

##  Appendix A Dilation Consistency Relation at Second Order

In this Appendix, we show that at second order in Newtonian gauge, the
response of the short wavelength density perturbation to the long wavelength
$\zeta_{L}$ is a pure dilation and when combined with the Maldacena
consistency relation there is no locally observable effect of a constant
$\zeta_{L}$. Since we focus on the $\mu$-distortion, mainly produced during
the radiation dominated epoch, we concretely calculate the second-order
perturbations in that era. Note that the calculation and notation here, which
differs from the main text, is based on Ref. Inomata (2021).

First, let us summarize our notation in this appendix. In Newtonian gauge, we
can write the scalar parts of the metric perturbations as

$\displaystyle ds^{2}$
$\displaystyle=g_{\mu\nu}dx^{\mu}dx^{\nu}=a^{2}\left\\{-(1+2\Phi^{(1)}+\Phi^{(2)})d\eta^{2}+\left[(1-2\Psi^{(1)}-\Psi^{(2)})\delta_{ij}\right]{dy^{i}dy^{j}}\right\\},$
(41)

where the superscript denotes the order in perturbations. Here, we assume a
perfect fluid for simplicity, which leads to $\Phi^{(1)}=\Psi^{(1)}$, and thus
ignore the correction from neutrinos. Then, we can express the energy-momentum
tensor as

$\displaystyle T^{\mu}_{\ \nu}$
$\displaystyle=(\rho+P)u^{\mu}u_{\nu}+P\delta^{\mu}_{\ \nu},$ (42)

where $\rho$ is the energy density, $P$ is the pressure, and $u_{\mu}$ is the
4-velocity. We take the following notation for their perturbations up to the
second order:

$\displaystyle\rho$
$\displaystyle=\rho^{(0)}+\delta\rho^{(1)}+\frac{1}{2}\delta\rho^{(2)},$ (43)
$\displaystyle P$ $\displaystyle=P^{(0)}+\delta P^{(1)}+\frac{1}{2}\delta
P^{(2)},$ (44) $\displaystyle u^{i}$ $\displaystyle=\frac{1}{a}\left(\delta
v^{{(1)},i}+\frac{1}{2}\delta
v^{{(2)},i}+\frac{1}{2}\delta{v^{i}_{V}}^{(2)}\right),$ (45)

where $\delta v$ is the velocity potential and $\delta v_{V}^{i}$ is the
vector part of the velocity perturbation. We define the density perturbation,
$\delta^{(a)}$, as

$\displaystyle\delta^{(a)}\equiv\frac{\delta\rho^{(a)}}{\rho^{(0)}}.$ (46)

At the first order in perturbations, this is related to the curvature
perturbation as

$\displaystyle\delta^{(1)}_{\bm{k}}(\eta)$
$\displaystyle=\displaystyle-\frac{2}{3}\zeta^{(1)}_{\bm{k}}T_{\delta}(x),$
(47)

where the subscript $\bm{k}$ represents the Fourier mode whereas perturbations
are in real space otherwise in this Appendix, $x\equiv k\eta$ with
$k\equiv|\bm{k}|$, and $\zeta^{(1)}=-(3/2)\Phi^{(1)}$ in the superhorizon
limit. The transfer function is defined as

$\displaystyle T_{\delta}(x)$
$\displaystyle\equiv\frac{6x(-6+x^{2})\cos(\tfrac{x}{\sqrt{3}})-12\sqrt{3}(-3+x^{2})\sin(\tfrac{x}{\sqrt{3}})}{x^{3}}.$
(48)

Next, let us calculate the second-order density perturbations. The density
perturbation is related to the other perturbations as Inomata (2021)

$\displaystyle\delta^{(2)}=$ $\displaystyle-2\Psi^{(2)}+2N^{j}_{\ i}{B^{i}_{\
j}}^{(2)}-\frac{2}{\mathcal{H}}{\Psi^{(2)}}^{\prime}+\frac{2}{3\mathcal{H}^{2}}\Psi^{{(2)},i}_{\quad\
,i}+\frac{2}{\mathcal{H}^{2}}\left({\Phi^{(1)}}^{\prime}\right)^{2}+\frac{16}{3\mathcal{H}^{2}}\Phi^{(1)}\Phi^{{(1)},i}_{\quad\
,i}$
$\displaystyle+\frac{1}{\mathcal{H}^{2}}\left(2-\frac{8}{9(1+w)}\right)\Phi^{{(1)},i}\Phi^{(1)}_{\quad,i}-\frac{8}{9(1+w)\mathcal{H}^{3}}(\Phi^{{(1)},i}\Phi^{(1)}_{\quad,i})^{\prime}-\frac{8}{9(1+w)\mathcal{H}^{4}}{\Phi^{(1)}}^{\prime}_{,i}{\Phi^{{(1)},i}}^{\prime},$
(49)

where the prime here denotes the derivative with respect to $\eta$,
$\mathcal{H}(\equiv a^{\prime}/a)$ is the conformal Hubble parameter, and $w$
is the equation of state parameter. $B^{i}_{\ j}$ and $N^{i}_{\ j}$ are
defined as

$\displaystyle{B^{i}_{\
j}}^{(2)}\equiv\left[\frac{4(5+3w)}{3(1+w)}\Phi^{{(1)},i}\Phi^{(1)}_{\quad,j}+\frac{8}{3(1+w)\mathcal{H}}\left(\Phi^{{(1)},i}{\Phi^{(1)}_{\quad,j}}\right)^{\prime}+\frac{8}{3(1+w)\mathcal{H}^{2}}{\Phi^{{(1)},i}}^{\prime}{\Phi^{(1)}_{\quad,j}}^{\prime}\right],$
(50) $\displaystyle N^{j}_{\ i}A^{i}_{\ j}(\bm{x})$
$\displaystyle\equiv\frac{3}{2}\nabla^{-2}\left(\frac{\partial^{j}\partial_{i}}{\nabla^{2}}-\frac{1}{3}\delta^{j}_{\
i}\right)A^{i}_{\ j}(\bm{x})$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}\left(-\frac{3}{2k^{2}}\right)\left(\frac{k^{j}k_{i}}{k^{2}}-\frac{1}{3}\delta^{j}_{\
i}\right){A_{\bm{k}}}^{i}_{\ j},$ (51)

where $A^{i}_{\ j}$ is an arbitrary tensor. Since $\delta^{(2)}$ depends on
$\Psi^{(2)}$, we need to calculate $\Psi^{(2)}$ first. The equation of motion
for $\Psi^{(2)}$ is Inomata (2021)

$\displaystyle{\Psi^{(2)}}^{\prime\prime}+3(1+c_{s}^{2})\mathcal{H}{\Psi^{(2)}}^{\prime}+\left[2\mathcal{H}^{\prime}+(3c_{s}^{2}+1)\mathcal{H}^{2}\right]\Psi^{(2)}-c_{s}^{2}\Psi^{{(2)},i}_{\quad\
,i}=S^{(2)},$ (52) $\displaystyle S^{(2)}\equiv$
$\displaystyle\left(3c_{s}^{2}-\frac{1}{3}\right)\Phi^{{(1)},i}\Phi^{(1)}_{\quad,i}+8c_{s}^{2}\Phi^{(1)}\Phi^{{(1)},i}_{\quad\
,i}+(3c_{s}^{2}+1)\left({\Phi^{(1)}}^{\prime}\right)^{2}+\left[(3c_{s}^{2}+1)\mathcal{H}^{2}+2\mathcal{H}^{\prime}\right]N^{j}_{\
i}{B^{i}_{\ j}}^{(2)}$ $\displaystyle\ +\mathcal{H}N^{j}_{\ i}\left({B^{i}_{\
j}}^{(2)}\right)^{\prime}+\frac{1}{3}N^{j}_{\ i}\left({B^{i}_{\
j}}^{(2)}\right)^{,k}_{\
,k}+\left(\frac{1}{3}-c_{s}^{2}\right)\frac{4}{3(1+w)\mathcal{H}^{2}}\left(\mathcal{H}\Phi^{{(1)},i}+{\Phi^{{(1)},i}}^{\prime}\right)\left(\mathcal{H}\Phi^{(1)}_{,i}+{\Phi^{(1)}}^{\prime}_{,i}\right),$
(53)

where $c_{s}$ is the sound speed.

Here, we focus on the second-order perturbations induced by a superhorizon
mode “$L$” on the much smaller mode $k_{L}/k\ll 1$. This case corresponds to
the $\mu$-distortion production in the presence of the superhorizon
perturbations. In this case, we can approximate Eq. (52) as

$\displaystyle{\Psi^{(2)}}^{\prime\prime}+4\mathcal{H}{\Psi^{(2)}}^{\prime}-\frac{1}{3}\Psi^{{(2)},i}_{\quad\
,i}\simeq\frac{8}{3}\Phi^{(1)}_{L}\Phi^{{(1)},i}_{\quad\ ,i},$ (54)

where we have substituted The other contributions in $S^{(2)}$ are sub-leading
because they include spatial derivative on the superhorizon mode “$L$”. In the
Fourier space, the equation becomes

$\displaystyle{\Psi^{(2)}}^{\prime\prime}_{\bm{k}}+\frac{4}{\eta}{\Psi^{(2)}}^{\prime}_{\bm{k}}+\frac{k^{2}}{3}\Psi_{\bm{k}}^{(2)}\simeq-\frac{8}{3}\Phi_{L}^{(1)}k^{2}\Phi_{\bm{k}}^{(1)},$
(55)

where we have assumed that the superhorizon perturbation as the constant
quantity, $\Phi_{L}$, because it does not evolve in the Newtonian gauge. Note
that $\Phi_{L}$ here and below is still a function of spatial comoving
coordinates but variations are small locally compared to the short wavelength
mode so we have ignored them in the Fourier transform. Then, we can solve
$\Psi^{(2)}_{\bm{k}}$ using the Green function method:

$\displaystyle\Psi^{(2)}_{\bm{k}}(\eta)$
$\displaystyle\simeq\Psi^{(2)}_{\bm{k}}(0)T(x)+\int^{\eta}_{0}d\bar{\eta}\left(\frac{a(\bar{\eta})}{a(\eta)}\right)^{2}G(k,\eta;\bar{\eta})\left(-\frac{8}{3}\Phi_{L}k^{2}\Phi^{(1)}_{\bm{k}}(\bar{\eta})\right),$
(56)

where the concrete expression of the Green function is given by

$\displaystyle
kG(k,\eta;\bar{\eta})=-\Theta(\eta-\bar{\eta})\frac{x\bar{x}}{\sqrt{3}}\left[j_{1}(x/\sqrt{3})y_{1}(\bar{x}/\sqrt{3})-j_{1}(\bar{x}/\sqrt{3})y_{1}(x/\sqrt{3})\right].$
(57)

The initial condition of $\Psi^{(2)}$ is given by

$\displaystyle\Psi^{(2)}_{\bm{k}}(0)=-\frac{16}{9}\zeta^{(1)}_{L}\zeta^{(1)}_{\bm{k}},$
(58)

where we have here assumed a Gaussian distribution of $\zeta$, which means
that the contribution from primordial non-Gaussianity from the Maldacena
consistency relation is not included in $\Psi^{(2)}$ (and $\delta^{(2)}$) in
this expression. We will independently take into account that contribution
later. Then, we obtain $\Psi^{(2)}$ as

$\displaystyle\Psi^{(2)}_{\bm{k}}(\eta)$
$\displaystyle\simeq\left[\frac{16\sqrt{3}}{3}\frac{\sin(x/\sqrt{3})}{x}+\frac{64\cos(x/\sqrt{3})}{x^{2}}\right]\zeta_{L}^{(1)}\zeta^{(1)}_{\bm{k}},$
(59)

where we have neglected the contribution of $\mathcal{O}(x^{-3})$ because we
want to know the evolution of the second-order density perturbations on
subhorizon scales, which is related to the $\mu$-distortion production. In the
large $x$ limit, Eq. (49) can be approximated in the Fourier space as

$\displaystyle\delta^{(2)}_{\bm{k}}$
$\displaystyle\simeq-\frac{2x^{2}}{3}\Psi^{(2)}_{\bm{k}}-2x\frac{d\Psi^{(2)}_{\bm{k}}}{dx}-\frac{16x^{2}}{3}\Phi_{L}^{(1)}\Phi^{(1)}_{\bm{k}}$
$\displaystyle\simeq\left(-\frac{32\sqrt{3}}{9}x\sin(x/\sqrt{3})-32\cos(x/\sqrt{3})\right)\zeta_{L}^{(1)}\zeta^{(1)}_{\bm{k}}.$
(60)

Then, we finally get the following expression for the total energy density
perturbation:

$\displaystyle\delta_{\bm{k}}(\eta)$
$\displaystyle=\delta^{(1)}_{\bm{k}}(\eta)+\frac{1}{2}\delta^{(2)}_{\bm{k}}(\eta)$
$\displaystyle\simeq-\frac{2}{3}\left(1+\frac{4}{3}\zeta^{(1)}_{L}\right)T_{\delta}\left(x\left(1-\frac{4}{3}\zeta^{(1)}_{L}\right)\right)\zeta^{(1)}_{\bm{k}},$
(61)

where the approximate equality is valid in $x\gg 1$.

In the following, we relate this result to the dilation transformation. The
dilation consistency relation in the Newtonian gauge is given by Horn  _et
al._ (2014)

$\displaystyle\lim_{\bm{q}\rightarrow
0}\frac{\left<{\pi_{\bm{q}}\delta_{\bm{k}_{1}}\cdots\delta_{\bm{k}_{N}}}\right>^{c}}{P_{\pi}(q)}\epsilon(\eta)=\sum^{N}_{a=1}\left[-\epsilon(\eta_{a})\left(\left.\frac{\bar{\rho}^{\prime}}{\bar{\rho}}\right|_{\eta_{a}}+\partial_{\eta_{a}}\right)+\lambda\left(3+\bm{k}_{a}\cdot\partial_{\bm{k}_{a}}\right)\right]\left<{\delta_{\bm{k}_{1}}\cdots\delta_{\bm{k}_{N}}}\right>^{c},$
(62)

where
$\left<{\pi_{\bm{q}}\pi_{\bm{p}}}\right>=(2\pi^{3})\delta_{D}(\bm{q}+\bm{p})P_{\pi}(q)$,
$\pi$ is the velocity potential ($\pi=\delta v$), and the superscript $c$ of
the braket means the connected part. The $\epsilon$ and $\lambda$ are the
coordinate transformation parameters to go from a gauge where the effect is a
pure dilation and are given by
$y^{\mu}\rightarrow{y^{\prime}}^{\mu}=y^{\mu}+\xi^{\mu}$ with
$\xi^{0}=\epsilon$ and $\xi^{i}=\lambda y^{i}$. During a radiation dominated
epoch, we can derive the following relations: $\lambda=-3\epsilon/\eta$,
$\bar{\rho}^{\prime}/\bar{\rho}=-4/\eta$, and $\pi\simeq\eta\zeta/3$ in the
superhorizon limit. Here, we take $\lambda=\zeta_{L}$ and $N=2$ to determine
the modulation of short wavelength power by $\zeta_{L}$, such that the
transformation is from comoving gauge to Newtonian gauge and rewrite the
consistency relation as

$\displaystyle\lim_{\bm{q}\rightarrow
0}\frac{\left<{\zeta_{\bm{q}}\delta_{\bm{k}_{1}}\delta_{\bm{k}_{2}}}\right>^{c}}{P_{\zeta}(q)}\zeta_{L}^{(1)}=\sum^{2}_{a=1}\left[\left(\frac{4}{3}\zeta_{L}^{(1)}-\frac{1}{3}\zeta_{L}^{(1)}\eta\,\partial_{\eta}\right)-\zeta_{L}^{(1)}\left(3+\bm{k}_{a}\cdot\partial_{\bm{k}_{a}}\right)\right]\left<{\delta_{\bm{k}_{1}}\delta_{\bm{k}_{2}}}\right>^{c}.$
(63)

Note again that $\zeta_{L}$ is in real space as $\zeta_{L}(\bm{y})$, but its
perturbation scale is much larger than the scales of $\bm{k}_{1}$ and
$\bm{k}_{2}$ so that we can approximately consider it constant except when
considering the correlation to the large scale mode $\bm{q}$. Using Eq. (61),
we can rewrite the left-hand side as

$\displaystyle\lim_{\bm{q}\rightarrow
0}\frac{\left<{\zeta_{\bm{q}}\delta_{\bm{k}_{1}}\delta_{\bm{k}_{2}}}\right>^{c}}{P_{\zeta}(q)}\zeta_{L}^{(1)}=$
$\displaystyle\frac{4}{9}\left[\left(1+\frac{8}{3}\zeta^{(1)}_{L}\right)T_{\delta}^{2}\left(x_{1}\left(1-\frac{4}{3}\zeta^{(1)}_{L}\right)\right)-T_{\delta}^{2}\left(x_{1}\right)\right](2\pi^{3})\delta_{D}(\bm{k}_{1}+\bm{k}_{2})P_{\zeta}(k_{1})$
$\displaystyle-\zeta_{L}^{(1)}\frac{d\ln\Delta_{\zeta}^{2}}{d\ln
k}\Big{|}_{k=k_{1}}(2\pi^{3})\delta_{D}(\bm{k}_{1}+\bm{k}_{2})P_{\zeta}(k_{1}),$
(64)

where we have used

$\displaystyle\left<{\zeta_{\bm{q}}^{(1)}\zeta_{L}^{(1)}(\bm{y})}\right>^{c}$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot\bm{y}}\left<{\zeta^{(1)}_{\bm{q}}\zeta^{(1)}_{\bm{k}}}\right>^{c}$
$\displaystyle=P_{\zeta}(q).$ (65)

The first line in Eq. (64) corresponds to the contribution from
$\delta^{(2)}$. On the other hand, the second line corresponds to the
Maldacena consistency relation Maldacena (2003), that is, it just comes from
the contribution proportional to
$\left<{\zeta_{\bm{q}}\zeta_{\bm{k}_{1}}\zeta_{\bm{k}_{2}}}\right>$, not
related to the evolution of the second-order scalar perturbations. On the
other hand, once we substitute

$\displaystyle\left<{\delta_{\bm{k}_{1}}\delta_{\bm{k}_{2}}}\right>^{c}=\frac{4}{9}T_{\delta}^{2}(x_{1})(2\pi^{3})\delta_{D}(\bm{k}_{1}+\bm{k}_{2})P_{\zeta}(k_{1})$
(66)

into the right-hand side of Eq. (63), we can see that the right-hand side is
the same as Eq. (64) at least at the lowest order in the perturbations,
$\mathcal{O}(\zeta_{L}P_{\zeta})$. From this, we can see that the dilation
consistency relation is satisfied, once both the second-order evolution and
primordial non-Gaussianity of the Maldacena relation are included, which
indicates that the constant $\zeta_{L}$ does not give locally observable
effects. Although our calculation is based on the perfect fluid assumption, we
can expect that the dilation consistency relation would be satisfied even in
imperfect fluid.

## References

  * Kawasaki  _et al._ (2016) M. Kawasaki, A. Kusenko, Y. Tada, and T. T. Yanagida, Phys. Rev. D94, 083523 (2016), arXiv:1606.07631 [astro-ph.CO] .
  * Garcia-Bellido and Ruiz Morales (2017) J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18, 47 (2017), arXiv:1702.03901 [astro-ph.CO] .
  * Ezquiaga  _et al._ (2018) J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales, Phys. Lett. B 776, 345 (2018), arXiv:1705.04861 [astro-ph.CO] .
  * Kannike  _et al._ (2017) K. Kannike, L. Marzola, M. Raidal, and H. Veermäe, JCAP 1709, 020 (2017), arXiv:1705.06225 [astro-ph.CO] .
  * Germani and Prokopec (2017) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6 (2017), arXiv:1706.04226 [astro-ph.CO] .
  * Ballesteros and Taoso (2018) G. Ballesteros and M. Taoso, Phys. Rev. D97, 023501 (2018), arXiv:1709.05565 [hep-ph] .
  * Inomata  _et al._ (2018) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 97, 043514 (2018), arXiv:1711.06129 [astro-ph.CO] .
  * Hertzberg and Yamada (2018) M. P. Hertzberg and M. Yamada, Phys. Rev. D 97, 083509 (2018), arXiv:1712.09750 [astro-ph.CO] .
  * Motohashi and Hu (2017) H. Motohashi and W. Hu, Phys. Rev. D 96, 063503 (2017), arXiv:1706.06784 [astro-ph.CO] .
  * Cheng  _et al._ (2019) S.-L. Cheng, W. Lee, and K.-W. Ng, Phys. Rev. D 99, 063524 (2019), arXiv:1811.10108 [astro-ph.CO] .
  * Byrnes  _et al._ (2019) C. T. Byrnes, P. S. Cole, and S. P. Patil, JCAP 06, 028 (2019), arXiv:1811.11158 [astro-ph.CO] .
  * Passaglia  _et al._ (2019) S. Passaglia, W. Hu, and H. Motohashi, Phys. Rev. D 99, 043536 (2019), arXiv:1812.08243 [astro-ph.CO] .
  * Drees and Xu (2021) M. Drees and Y. Xu, Eur. Phys. J. C 81, 182 (2021), arXiv:1905.13581 [hep-ph] .
  * Carrilho  _et al._ (2019) P. Carrilho, K. A. Malik, and D. J. Mulryne, Phys. Rev. D 100, 103529 (2019), arXiv:1907.05237 [astro-ph.CO] .
  * Passaglia  _et al._ (2020) S. Passaglia, W. Hu, and H. Motohashi, Phys. Rev. D 101, 123523 (2020), arXiv:1912.02682 [astro-ph.CO] .
  * Ng and Wu (2021) K.-W. Ng and Y.-P. Wu, (2021), arXiv:2102.05620 [astro-ph.CO] .
  * Inomata  _et al._ (2021a) K. Inomata, E. Mcdonough, and W. Hu, (2021a), arXiv:2104.03972 [astro-ph.CO] .
  * Inomata  _et al._ (2021b) K. Inomata, E. McDonough, and W. Hu, (2021b), arXiv:2110.14641 [astro-ph.CO] .
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Clesse and García-Bellido (2017) S. Clesse and J. García-Bellido, Phys. Dark Univ. 15, 142 (2017), arXiv:1603.05234 [astro-ph.CO] .
  * Sasaki  _et al._ (2016) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016), arXiv:1603.08338 [astro-ph.CO] .
  * Kashlinsky (2016) A. Kashlinsky, Astrophys. J. Lett. 823, L25 (2016), arXiv:1605.04023 [astro-ph.CO] .
  * Kashlinsky  _et al._ (2018) A. Kashlinsky, R. G. Arendt, F. Atrio-Barandela, N. Cappelluti, A. Ferrara, and G. Hasinger, Rev. Mod. Phys. 90, 025006 (2018), arXiv:1802.07774 [astro-ph.CO] .
  * García-Bellido  _et al._ (2021) J. García-Bellido, J. F. Nuño Siles, and E. Ruiz Morales, Phys. Dark Univ. 31, 100791 (2021), arXiv:2010.13811 [astro-ph.CO] .
  * Josan  _et al._ (2009) A. S. Josan, A. M. Green, and K. A. Malik, Phys. Rev. D79, 103520 (2009), arXiv:0903.3184 [astro-ph.CO] .
  * Assadullahi and Wands (2010) H. Assadullahi and D. Wands, Phys. Rev. D81, 023527 (2010), arXiv:0907.4073 [astro-ph.CO] .
  * Inomata and Nakama (2019) K. Inomata and T. Nakama, Phys. Rev. D99, 043511 (2019), arXiv:1812.00674 [astro-ph.CO] .
  * Kawasaki  _et al._ (2021) M. Kawasaki, H. Nakatsuka, and K. Nakayama, (2021), arXiv:2110.12620 [astro-ph.CO] .
  * Sunyaev and Zeldovich (1970) R. A. Sunyaev and Y. B. Zeldovich, Astrophys. Space Sci. 7, 3 (1970).
  * Daly (1991) R. A. Daly, Astrophys. J. 371, 14 (1991).
  * Hu  _et al._ (1994) W. Hu, D. Scott, and J. Silk, Astrophys. J. Lett. 430, L5 (1994), arXiv:astro-ph/9402045 .
  * Chluba and Sunyaev (2012) J. Chluba and R. A. Sunyaev, Mon. Not. Roy. Astron. Soc. 419, 1294 (2012), arXiv:1109.6552 [astro-ph.CO] .
  * Chluba (2013) J. Chluba, Mon. Not. Roy. Astron. Soc. 436, 2232 (2013), arXiv:1304.6121 [astro-ph.CO] .
  * Nakama  _et al._ (2017) T. Nakama, B. Carr, and J. Silk, (2017), arXiv:1710.06945 [astro-ph.CO] .
  * Carr  _et al._ (2020) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778 [astro-ph.CO] .
  * Chluba  _et al._ (2021) J. Chluba _et al._ , Exper. Astron. 51, 1515 (2021), arXiv:1909.01593 [astro-ph.CO] .
  * Hazumi  _et al._ (2020) M. Hazumi _et al._ (LiteBIRD), Proc. SPIE Int. Soc. Opt. Eng. 11443, 114432F (2020), arXiv:2101.12449 [astro-ph.IM] .
  * Pajer and Zaldarriaga (2012) E. Pajer and M. Zaldarriaga, Phys. Rev. Lett. 109, 021302 (2012), arXiv:1201.5375 [astro-ph.CO] .
  * Ganc and Komatsu (2012) J. Ganc and E. Komatsu, Phys. Rev. D 86, 023518 (2012), arXiv:1204.4241 [astro-ph.CO] .
  * Pajer and Zaldarriaga (2013) E. Pajer and M. Zaldarriaga, JCAP 02, 036 (2013), arXiv:1206.4479 [astro-ph.CO] .
  * Chluba  _et al._ (2017) J. Chluba, E. Dimastrogiovanni, M. A. Amin, and M. Kamionkowski, Mon. Not. Roy. Astron. Soc. 466, 2390 (2017), arXiv:1610.08711 [astro-ph.CO] .
  * Cabass  _et al._ (2018) G. Cabass, E. Pajer, and D. van der Woude, JCAP 08, 050 (2018), arXiv:1805.08775 [astro-ph.CO] .
  * Remazeilles and Chluba (2018) M. Remazeilles and J. Chluba, Mon. Not. Roy. Astron. Soc. 478, 807 (2018), arXiv:1802.10101 [astro-ph.CO] .
  * Remazeilles  _et al._ (2021) M. Remazeilles, A. Ravenni, and J. Chluba, (2021), arXiv:2110.14664 [astro-ph.CO] .
  * Maldacena (2003) J. M. Maldacena, JHEP 05, 013 (2003), arXiv:astro-ph/0210603 .
  * Pajer  _et al._ (2013) E. Pajer, F. Schmidt, and M. Zaldarriaga, Phys. Rev. D 88, 083502 (2013), arXiv:1305.0824 [astro-ph.CO] .
  * Özsoy and Tasinato (2021a) O. Özsoy and G. Tasinato, (2021a), arXiv:2111.02432 [astro-ph.CO] .
  * Özsoy and Tasinato (2021b) O. Özsoy and G. Tasinato, Phys. Rev. D 104, 043526 (2021b), arXiv:2104.12792 [astro-ph.CO] .
  * Namjoo  _et al._ (2013) M. H. Namjoo, H. Firouzjahi, and M. Sasaki, EPL 101, 39001 (2013), arXiv:1210.3692 [astro-ph.CO] .
  * Martin  _et al._ (2013) J. Martin, H. Motohashi, and T. Suyama, Phys. Rev. D 87, 023514 (2013), arXiv:1211.0083 [astro-ph.CO] .
  * Cai  _et al._ (2018) Y.-F. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D.-G. Wang, and Z. Wang, JCAP 05, 012 (2018), arXiv:1712.09998 [astro-ph.CO] .
  * Kinney (2005) W. H. Kinney, Phys. Rev. D 72, 023515 (2005), arXiv:gr-qc/0503017 .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Hu  _et al._ (2016) W. Hu, C.-T. Chiang, Y. Li, and M. LoVerde, Phys. Rev. D 94, 023002 (2016), arXiv:1605.01412 [astro-ph.CO] .
  * Hu and Joyce (2017) W. Hu and A. Joyce, Phys. Rev. D 95, 043529 (2017), arXiv:1612.02454 [astro-ph.CO] .
  * Sirko (2005) E. Sirko, Astrophys. J. 634, 728 (2005), arXiv:astro-ph/0503106 .
  * Li  _et al._ (2014) Y. Li, W. Hu, and M. Takada, Phys. Rev. D 89, 083519 (2014), arXiv:1401.0385 [astro-ph.CO] .
  * Hu and White (1996) W. Hu and M. J. White, Astrophys. J. 471, 30 (1996), arXiv:astro-ph/9602019 .
  * Hu and Sugiyama (1996) W. Hu and N. Sugiyama, Astrophys. J. 471, 542 (1996), arXiv:astro-ph/9510117 .
  * Danese and de Zotti (1982) L. Danese and G. de Zotti, Astronomy and Astrophysics 107, 39 (1982).
  * Hu and Silk (1993) W. Hu and J. Silk, Phys. Rev. D 48, 485 (1993).
  * Fixsen  _et al._ (1996) D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, Astrophys. J. 473, 576 (1996), arXiv:astro-ph/9605054 .
  * Kogut  _et al._ (2011) A. Kogut _et al._ , JCAP 07, 025 (2011), arXiv:1105.2044 [astro-ph.CO] .
  * Kogut and Fixsen (2020) A. Kogut and D. J. Fixsen, JCAP 05, 041 (2020), arXiv:2002.00976 [astro-ph.IM] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Ravenni  _et al._ (2017) A. Ravenni, M. Liguori, N. Bartolo, and M. Shiraishi, JCAP 09, 042 (2017), arXiv:1707.04759 [astro-ph.CO] .
  * Inomata (2021) K. Inomata, JCAP 03, 013 (2021), arXiv:2008.12300 [gr-qc] .
  * Horn  _et al._ (2014) B. Horn, L. Hui, and X. Xiao, JCAP 09, 044 (2014), arXiv:1406.0842 [hep-th] .