# The Atacama Cosmology Telescope: Galactic Dust Structure and the Cosmic PAH
Background in Cross-correlation with WISE

Rodrigo Córdova Rosado Department of Astrophysical Sciences, Peyton Hall,
Princeton University, 4 Ivy Lane, Princeton, NJ, 08544  Brandon S. Hensley
Department of Astrophysical Sciences, Peyton Hall, Princeton University, 4 Ivy
Lane, Princeton, NJ, 08544  Susan E. Clark Department of Physics, Stanford
University, Stanford, CA, 94305  Kavli Institute for Particle Astrophysics &
Cosmology, Stanford University, PO Box 2450, Stanford, CA, 94305  Adriaan J.
Duivenvoorden Center for Computational Astrophysics, Flatiron Institute, 162
5th Avenue, New York, NY 10010 USA  Zachary Atkins Joseph Henry Laboratories
of Physics, Jadwin Hall, Princeton University, Princeton, NJ, USA 08544  Elia
Stefano Battistelli Sapienza University of Rome, Physics Department, Piazzale
Aldo Moro 5, 00185 Rome, Italy  Steve K. Choi Department of Physics, Cornell
University, Ithaca, NY 14853, USA  Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Jo Dunkley Joseph Henry Laboratories of
Physics, Jadwin Hall, Princeton University, Princeton, NJ, USA 08544
Department of Astrophysical Sciences, Peyton Hall, Princeton University, 4 Ivy
Lane, Princeton, NJ, 08544  Carlos Hervías-Caimapo Instituto de Astrofísica
and Centro de Astro-Ingeniería, Facultad de Fìsica, Pontificia Universidad
Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile
Zack Li Canadian Institute for Theoretical Astrophysics, University of
Toronto, Toronto, ON, Canada M5S 3H8  Thibaut Louis  Université Paris-Saclay,
CNRS/IN2P3, IJCLab, 91405 Orsay, France  Sigurd Naess Institute of Theoretical
Astrophysics, University of Oslo, Norway  Lyman A. Page Joseph Henry
Laboratories of Physics, Jadwin Hall, Princeton University, Princeton, NJ, USA
08544  Bruce Partridge Department of Physics and Astronomy, Haverford College,
Haverford, PA, USA 19041  Cristóbal Sifón Instituto de Física, Pontificia
Universidad Católica de Valparaíso, Casilla 4059, Valparaíso, Chile  Suzanne
T. Staggs Joseph Henry Laboratories of Physics, Jadwin Hall, Princeton
University, Princeton, NJ, USA 08544  Cristian Vargas Instituto de Astrofísica
and Centro de Astro-Ingeniería, Facultad de Fìsica, Pontificia Universidad
Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile
Edward J. Wollack NASA/Goddard Space Flight Center, Greenbelt, MD, USA 20771

###### Abstract

We present a cross-correlation analysis between $1\arcmin$ resolution total
intensity and polarization observations from the Atacama Cosmology Telescope
(ACT) at 150 and 220 GHz and $15\arcsec$ mid-infrared photometry from the
Wide-field Infrared Survey Explorer (WISE) over 107 $12.5^{\circ}\times
12.5^{\circ}$ patches of sky. We detect a spatially isotropic signal in the
WISE$\times$ACT $TT$ cross power spectrum at $30\sigma$ significance that we
interpret as the correlation between the cosmic infrared background at ACT
frequencies and polycyclic aromatic hydrocarbon (PAH) emission from galaxies
in WISE, i.e., the cosmic PAH background. Within the Milky Way, the Galactic
dust $TT$ spectra are generally well-described by power laws in $\ell$ over
the range $10^{3}<\ell<10^{4}$, but there is evidence both for variability in
the power law index and for non-power law behavior in some regions. We measure
a positive correlation between WISE total intensity and ACT $E$-mode
polarization at $1000<\ell\lesssim 6000$ at $>3\sigma$ in each of 35 distinct
$\sim$100 deg2 regions of the sky, suggesting alignment between Galactic
density structures and the local magnetic field persists to sub-parsec
physical scales in these regions. The distribution of $TE$ amplitudes in this
$\ell$ range across all 107 regions is biased to positive values, while there
is no evidence for such a bias in the $TB$ spectra. This work constitutes the
highest-$\ell$ measurements of the Galactic dust $TE$ spectrum to date and
indicates that cross-correlation with high-resolution mid-infrared
measurements of dust emission is a promising tool for constraining the spatial
statistics of dust emission at millimeter wavelengths.

††facilities: Planck, WISE††software: Astropy (Astropy Collaboration et al.,
2013, 2018), emcee (Foreman-Mackey et al., 2013, 2019), Matplotlib (Hunter,
2007), NaMaster (Alonso et al., 2019), NumPy (van der Walt et al., 2011),
pixell (Naess et al., 2021), SciPy (Virtanen et al., 2020)

##  1 Introduction

The interstellar medium (ISM) is a turbulent environment. Energy is injected
at large physical scales by processes like stellar feedback, and a complex
turbulent energy cascade shapes the ISM over a vast range of physical scales
(e.g., Ferrière, 2001; Elmegreen & Scalo, 2004). High-dynamic range
observations of interstellar emission are critical for understanding the flow
of mass and energy in the ISM (e.g., Fissel et al., 2019a; Stinebring et al.,
2019). One common approach is to measure the power spectrum of ISM emission,
using tracers like neutral hydrogen (H I) emission or interstellar dust. These
power spectra are often found to be well-described by a power law, with a
power spectral index that can be compared to theoretical predictions (e.g.,
Crovisier & Dickey, 1983; Miville-Deschênes et al., 2003a; Miville-Deschênes &
Lagache, 2005; Miville-Deschênes et al., 2007; Martin et al., 2010, 2015;
Blagrave et al., 2017; Pingel et al., 2022).

Combining observations of dust emission from Planck and the Wide-field
Infrared Survey Explorer (WISE; Wright et al., 2010) with MegaCam measurements
of optical scattering from dust, Miville-Deschênes et al. (2016) demonstrated
that the dust power spectrum in total intensity (i.e., $TT$) is well-fit by a
power law $k^{-2.9\pm 0.1}$ from scales of degrees to $\sim 1^{\prime\prime}$,
corresponding to physical scales of $\sim 0.01$ pc. The power spectral index
of the observable column density is related to the statistics of the 3D
density field, which are in turn affected by turbulence and the phase
distribution of the gas, but are not directly measurable due to projection
effects (Miville-Deschênes et al., 2003b; Clark et al., 2019; Kalberla & Haud,
2019; Marchal & Miville-Deschênes, 2021). On scales greater than $5\arcmin$,
this index is consistent with measurements across the sky from the Planck
satellite (Planck Collaboration XXX, 2014; Planck Collaboration Int. XLVIII,
2016). In this work, we investigate the variability of this power law index at
smaller scales ($10^{3}<\ell<10^{4}$, corresponding to
$10\arcmin>\theta>1\arcmin$).

Planck observations established a robust positive correlation between dust
total intensity and dust $E$-mode polarization (i.e., $TE$) for multipoles
$\ell\lesssim 600$ over much of the sky (Planck Collaboration XI, 2020). Such
a correlation is expected if dust-bearing ISM structures are elongated along
magnetic field lines (Zaldarriaga, 2001; Huffenberger et al., 2020; Clark et
al., 2021). Indeed, H I filaments are ubiquitous across the sky with
orientations that are strongly correlated with the measured dust polarization
angles (Clark et al., 2015; Clark & Hensley, 2019). H I structure can thus be
used to measure properties of Galactic dust polarization in cross-correlation
(BICEP/Keck Collaboration et al., 2023; Halal et al., 2023). Filament-based
models successfully reproduce the observed $TE$ correlation (Hervías-Caimapo &
Huffenberger, 2022). However, in dense regions, it is observed that ISM
structures are preferentially oriented perpendicular to magnetic field lines
(i.e., negative $TE$; Planck Collaboration Int. XXXV, 2016; Bracco et al.,
2019). In this work, we extend the characterization of the $TE$ correlation to
smaller scales ($1000<\ell\lesssim 6000$).

To probe the small-scale $TT$ and $TE$ spectra of Galactic dust emission, we
employ new maps of millimeter dust emission from the Atacama Cosmology
Telescope (ACT) in both total and polarized intensity. The combination of
sensitivity, angular resolution, and sky coverage (sky fraction $f_{\rm
sky}\simeq 40\%$) afforded by ACT observations enable characterization of dust
at arcminute scales over a large sky area, and thus our investigation of the
universality of the small-scale dust power spectrum across a range of Galactic
environments. We employ ACT maps from two broad bands centered roughly at 150
and 220 GHz, with angular resolutions of $1.4\arcmin$ and $1.0\arcmin$,
respectively.

We complement the ACT data with full-sky observations of mid-infrared emission
from WISE. In particular, we use the maps of diffuse emission extracted from
observations in the W3 passband by Meisner & Finkbeiner (2014) (hereafter
MF14). Their custom processing of the WISE data removes compact sources and
associated data artifacts. The resulting map covers the full sky at
$15^{\prime\prime}$ resolution with noise properties independent of ACT,
enabling correlation analysis down to the ACT resolution limit.

The broad W3 passband is centered at 12 $\mu$m but has appreciable spectral
response from $\simeq$8–16 $\mu$m (Wright et al., 2010). For regions typical
of the diffuse ISM of the Galaxy, emission at these wavelengths is dominated
by the 7.7, 11.3, and 12.7 $\mu$m mid-infrared emission features associated
with polycyclic aromatic hydrocarbons (PAHs; e.g., Tielens, 2008). Continuum
emission is also present from PAHs and likely other nanoparticles undergoing
single photon heating. Other expected sources of diffuse emission in this map
include the Zodiacal light and the extragalactic background light.

As a ubiquitous component of the Galactic dust population, PAHs are generally
well-coupled to the larger grains responsible for far-infrared and millimeter
emission (e.g., Onaka et al., 1996; Mattila et al., 1996; Draine & Li, 2007).
Empirically, strong correlations have been observed between the MF14 map and
dust emission in the Planck bands (Hensley et al., 2016). However, the mass
fraction of dust in the form of PAHs varies over the sky (Planck Collaboration
Int. XXIX, 2016; Hensley et al., 2016) with an apparent dependence on ISM
phase (Hensley et al., 2022). Further, the PAH emission spectrum is sensitive
to the intensity and spectrum of the interstellar radiation field (Draine et
al., 2021). Thus, PAH emission and millimeter dust emission will not be
perfectly correlated in detail.

In this work, we assess the ability of high resolution observations of mid-
infrared dust emission to correlate with—and thus predict—dust emission
properties at millimeter wavelengths. We then use this combination of
independent datasets to characterize the Galactic dust power spectrum at small
scales both in intensity and polarization, as well as its variation across the
sky.

The MF14 WISE W3 map isolates all diffuse emission falling into the W3
passband. In principle, this includes emission from unresolved galaxies across
cosmic time. Chiang & Ménard (2019) found that the MF14 map has statistically
significant correlation with optical measurements of galaxies and active
galactic nuclei from redshifts $z\lesssim 2$. They interpret this as
redshifted PAH emission from galaxies, i.e., the cosmic PAH background.
Indeed, the 7.7 $\mu$m PAH feature is the strongest of the PAH features and
remains in the W3 passband until $z\sim$ 1\. In this work, we find robust
evidence for a spatially isotropic correlation between the MF14 WISE W3 map
and ACT maps at both 150 and 220 GHz that we interpret as the first detection
of the cosmic PAH background in cross-correlation with the Cosmic Infrared
Background (CIB).

This paper is organized as follows. In Section 2, we summarize the data
products used in this analysis. In Section 3, we outline our methodologies for
power spectrum computation, uncertainty quantification, and parameter
estimation. We present the results of our total intensity $\times$ total
intensity (“$TT$”) analysis and our total intensity $\times$ $E$-mode
polarization (“$TE$”) analysis in Sections 4 and 5, respectively. In Section 6
we discuss the implications of our results on both the structure of the ISM
and on the cosmic PAH background, and we conclude in Section 7.

##  2 Data

Figure 1: Locations of the 107 MF14 WISE tiles analyzed in this work overlaid
on the MF14 W3 diffuse emission map. Each $\sim 11^{\circ}$ diameter circle
represents the analysis mask adopted in the present study. Areas outside the
ACT footprint are rendered in grayscale. White regions, such as the prominent
Moon contamination features, are pixels masked by MF14.

###  2.1 ACT

ACT has measured the total intensity and linear polarization of the millimeter
sky over 18,000 deg2. ACT observed in five passbands: f030 (22–32 GHz), f040
(29–48 GHz), f090 (79–112 GHz), f150 (124–171 GHz) and f220 (171–276 GHz). In
this paper we make use of I, Q, and U Stokes maps made from f150 and f220
observations, with angular resolutions of $1.4^{\prime}$ and $1.0^{\prime}$,
respectively (Henderson et al., 2016; Thornton et al., 2016; Naess et al.,
2020). All of the f150 maps employed in this analysis co-add the three
detector arrays that observed in the band; only one array observed in the f220
band.

We use maps made using the f150 and f220 data from the 2008–2019 observing
seasons111These 2008–19 maps are not a major data release for ACT, but the DR5
2008-18 maps are available, and the upcoming DR6 maps from the 2008–22 data
will be made public. We compare the measured W3 $\times$ ACT f150 and f220
$TT$ and $TE$ spectra from the 2008–2019 dataset used in this analysis to the
spectra computed with the publicly available 2008–2018 DR5 dataset (Naess et
al., 2020) in two regions. Uncertainties on the power spectra decrease by
$\sim$20% with the inclusion of the 2019 data. A model fit yields parameters
that change by less than $1\sigma$ for both $TT$ and $TE$ fits between the
datasets. following a similar process used in making the publicly available
Data Release 5 (DR5)
maps222https://lambda.gsfc.nasa.gov/product/act/actpol_prod_table.html, which
only used data taken through 2018 (Naess et al., 2020). These maps use the
same Plate Carrée (CAR) projection as the DR5 maps and have a pixel size of
approximately $0\arcmin.5$. These maps do not include any Planck data directly
but have been scaled by an overall multiplicative factor that minimizes the
residual relative to the Planck temperature power spectrum between
$1000<\ell<2000$ (Naess et al., 2020). We use maps that include data taken
both during the night and day, and from which point sources detected at
greater than $5\sigma$ have been subtracted, resulting in a subtraction
threshold of $\sim$15 mJy in the f150 band.

We use the ACT-produced ACT+Planck co-added $Q$ and $U$ maps at f150 and f220
for our polarization analysis to enhance our sensitivity at $\ell\lesssim
2000$. The maps are constructed by combining the co-added ACT f150 and f220
maps described above with the Planck PR2 and PR3 143 and 217 GHz maps (Planck
Collaboration I, 2016, 2020) using the co-adding procedure described in Naess
et al. (2020).

The maps have associated products including inverse variance maps, which
include an estimate of the inverse variance in $\mu{\rm K}^{-2}$ per pixel, as
well as beam window functions and passbands. The maps are converted from CMB
temperature units to ${\rm MJy\,sr^{-1}}$ using conversion factors 395 and 482
$\rm{MJy~{}sr^{-1}}~{}{\rm K}^{-1}$ for f150 and f220, respectively, derived
assuming the approximate ACT band centers for a Rayleigh-Jeans spectrum (150
and 220 GHz, respectively). We consider scales $\ell>1000$ for the f150 total
intensity analysis and all polarization analyses, but only $\ell>1400$ for the
f220 total intensity analysis to avoid scales that are expected to suffer from
a loss of power due to modeling errors in the mapmaking process (Naess &
Louis, 2022). We use the HEALPix/COSMO polarization
convention333https://lambda.gsfc.nasa.gov/product/about/pol_convention.html
throughout. The overall systematic uncertainties are at the level of 10%
(Naess et al., 2020), which includes both map calibration uncertainties and
unit conversion errors from use of the ACT band centers rather than full
bandpass integration.

###  2.2 Planck

We use the Planck 353 GHz maps444https://pla.esac.esa.int/#maps produced with
the NPIPE data processing pipeline (Planck Collaboration Int. LVII, 2020) to
extend our analysis to higher frequencies and provide a point of comparison
with the ACT measurements. The Planck 353 GHz channel is the highest Planck
frequency designed to measure linear polarization and is the channel most
sensitive to polarized thermal dust emission. The maps are provided in units
of ${\rm K}$, which we convert to $\rm{MJy~{}sr^{-1}}$ via a conversion factor
of 287.5 $\rm{MJy~{}sr^{-1}}~{}{\rm K}^{-1}$ (Planck Collaboration III, 2020).

###  2.3 WISE

The Wide-field Infrared Survey Explorer (WISE) satellite observed the full sky
in four bands across the mid-infrared (MIR; Wright et al., 2010). The WISE W3
band is centered at $\sim$12$\,\micron$ and spans the $8.6$, $11.3$ and
$12.7\,\micron$ vibrational emission features of PAHs that dominate Galactic
cirrus emission at these frequencies (e.g., Mattila et al., 1996; Ingalls et
al., 2011).

MF14 produced a map of diffuse Galactic emission from the full-sky W3 data by
modeling and subtracting emission from point sources, solar system objects,
diffraction spikes, compact sources, and Moon and zodiacal light
contamination. MF14 combined the W3 data with Planck 857 GHz (350 $\micron$)
data to recover extended emission at scales greater than $2^{\circ}$. The MF14
data are provided as a set of 430 $12.5^{\circ}~{}\times~{}12.5^{\circ}$ tiles
smoothed to an angular resolution of $15\arcsec$ FWHM. Approximately 200 of
the 430 tiles overlap significantly ($>50\%$) with the ACT footprint. We
smooth these data to a final resolution of $45\arcsec$ before reprojecting
onto the ACT pixellization and CAR Projection using the pixell555Available
online at www.github.com/simonsobs/pixell package (Naess et al., 2021).

The WISE tiles are provided as fluxes $F_{\it W3}$ in units of digital number
(DN), which we convert to MJy sr-1 following Cutri et al. (2012):

$\displaystyle I_{\nu}$ $\displaystyle=\frac{F_{\nu 0}}{\theta_{\rm
pix}^{2}}10^{-\left(M_{0,\rm inst}+\Delta m-8.926\right)/2.5}\ \frac{F_{\it
W3}}{\rm DN}$ $\displaystyle=0.0135\frac{F_{\it W3}}{\rm DN}\ {\rm
MJy\,sr}^{-1}~{}~{}~{},$ (1)

where $F_{\nu 0}=31.674$ Jy is the W3 zero-magnitude flux density of a source
with $F_{\nu}\propto\nu^{-2}$ (Jarrett et al., 2011), $\theta_{\rm
pix}=2.75\arcsec$ is the W3 pixel size (Mainzer et al., 2005), $M_{0,\rm
inst}=17.800$ is the instrumental zero point (Cutri et al., 2012), $\Delta
m=5.174$ is the conversion from the WISE Vega-system magnitudes to AB
magnitudes (Jarrett et al., 2011), and 8.926 is the factor relating flux
density in Janskies to AB magnitudes (Tokunaga & Vacca, 2005). In detail, the
conversion between DN and MJy sr-1 depends on the spectrum of the
source—Equation (1) is strictly accurate only for spectra with
$F_{\nu}\propto\nu^{-2}$. However, the conversion factor differs by $\lesssim
13$% for power law indices between $-3$ and 3 (Wright et al., 2010).

##  3 Methodology

###  3.1 Power spectrum estimation

We compute the angular power spectra in each tile using standard partial sky
pseudo-$C_{\ell}$ methods (e.g., Hivon et al., 2002). We compute the $TT$,
$TE$, and $TB$ cross-spectra between the WISE I total intensity map, and the
ACT I, Q, and U maps, using the nawrapper666Introduced in Li et al. (2020) and
available through GitHub at https://github.com/xzackli/nawrapper interface to
the NaMaster software (Alonso et al., 2019). The mask is described in Section
3.2. We bin the measurements with $\ell(\ell+1)$ weighting, with equally-
spaced logarithmic bins from $\ell=1000$ to $10,000$ for f150 and $\ell=1400$
to $10,000$ for f220 (see Section 2.1).

To compute the uncertainty on a binned cross-spectrum $C_{b}^{xy}$, where $x$
and $y$ are any of $T$, $E$, or $B$, we start from the analytic expression
(Knox, 1995; Hivon et al., 2002)

$\sigma^{2}\left(C_{b}^{xy}\right)=\frac{(C_{b}^{xy})^{2}+(C_{b}^{xx}C_{b}^{yy})}{f_{\rm
sky}\frac{w_{2}^{2}}{w_{4}}(2\ell+1)\Delta\ell}~{}~{}~{},$ (2)

where $f_{\rm sky}$ is the fractional area of the sky, $\Delta\ell$ the bin-
width, and $\ell$ the bin midpoint. The correction factor
${w_{2}^{2}}/{w_{4}}$ accounts for the apodization of the mask. The $w_{i}$
factors are defined as $w_{i}\equiv\sum\Omega_{j}W_{j}^{i}$, where
$\Omega_{j}$ is the pixel area and $W_{j}$ is the value of the apodized mask
in pixel $j$ (Hivon et al., 2002). The auto power spectrum term is the sum of
a signal and a noise term (i.e., $C_{b}^{xx}=S_{b}^{xx}+N_{b}^{xx}$), while
the cross power spectrum is signal-only (i.e., $C_{b}^{xy}=S_{b}^{xy}$).

However, Equation (2) neglects the effects of mode-coupling induced by the
mask. To account for this, we employ the NaMaster implementation of analytic
methods to compute the full covariance matrix (García-García et al., 2019).
Since calculation of the full covariance matrix for all the cross-spectra used
in our analysis is computationally expensive, we compute it only for the WISE
W3 $\times$ f150 spectrum in each region in $TT$. We then approximate the per-
band uncertainties of each of the other spectra ($TT$ at f220 and both $TE$
and $TB$ at f150, f220, and 353 GHz) using the simpler expression in Equation
(2) rescaled by the ratio between the uncertainty computed from the diagonal
of the full covariance matrix and from Equation (2) at f150. This increases
the uncertainties by $\sim$10%, with some scale dependence. Note that only the
diagonal entries of the covariance matrix are used for plotting error bars and
in parameter estimation.

We apply the same correction derived from the $TT$ spectrum to both the $TE$
and $TB$ uncertainties. Because the same sky mask is used for all frequencies
and for all of $TT$, $TE$, and $TB$, this is a good approximation.

To further ensure the robustness of this approach, we exclude regions from our
analysis for which the off-diagonal terms of the covariance matrix are large.
Specifically, if including the off-diagonal elements of the covariance matrix
changes the $\chi^{2}$ of our best-fit model (see Section 3.3) by more than 2
for 20 degrees of freedom, then that region is discarded. In practice, this
removes regions with particularly complicated masks, e.g., many disconnected
regions.

Equation (2) accounts for the contribution of sample variance to the total
uncertainty. It would be appropriate to neglect sample variance in reporting
the cross-power spectrum of emission in a given region, since the error bars
are reflecting only how well the particular spectrum of that region is being
measured (see, e.g., the discussion in Planck Collaboration Int. XXX, 2016).
In contrast, sample variance should be included when fitting a model for the
underlying spectrum from which the observed spectrum is drawn. As such model
fitting is a principal focus of our analysis, we include sample variance in
all error bars throughout this work. Sample variance is typically $\lesssim
10\%$ of the total uncertainty at all $\ell$.

###  3.2 Masking

Our primary mask for each region is a circle in R.A. and Decl. with a diameter
of approximately 11∘. We apodize this mask with the `C1` cosine taper
implemented in NaMaster using an apodization scale of $1^{\circ}$.

The MF14 W3 data include a bitmask of data quality flags. We mask all pixels
affected by saturated point sources, first and second latent points of source
ghosts, line-like effects, and Moon and solar system object contamination
(corresponding to Flags 0, 3, 8, 14, 15, 18, and 20), as well as all pixels
without a flux measurement. Point sources brighter than 15 mJy at f150,
predominantly active galactic nuclei (AGN Marsden et al., 2014), are masked
with circular holes of radius $5\arcmin$. We apply an additional mask of 14
extended sources identified in ACT maps to mitigate contamination from objects
like planetary nebulae and resolved galaxies. This mask was created for the
upcoming ACT DR6 power spectrum analysis by visual inspection of the maps
after the initial source mask was applied.

We do not mask pixels flagged by MF14 as “compact resolved sources” (Flag 7).
These are mostly galaxies, most of which would not be detected in the ACT
maps. We do not wish to remove extragalactic signal, particularly in an
inhomogeneous way across the sky, as this would complicate our analysis of
extragalactic cross-correlations. On the other hand, it is likely that Flag 7
also identifies some compact resolved Galactic sources that, if retained,
would affect our power spectra particularly at high-$\ell$. On balance, we
prefer to leave all pixels flagged by Flag 7 in our analysis but note that a
careful separation of Galactic versus extragalactic compact resolved sources
could improve the analysis presented in this work.

We apodize the source and artifact mask at an apodization scale of $18\arcmin$
using the `C1` cosine taper implemented in NaMaster. The $18\arcmin$
apodization scale follows the nominal choices for the upcoming ACT CMB power
spectrum analysis. Our final mask combines the source and artifact mask with
the primary circular mask. Each region has a single mask that is used for all
analyses at all frequencies.

We select WISE tiles that overlap entirely with the ACT footprint, that retain
at least 33 deg2 of sky area after the mask is applied, and that pass the
covariance matrix criterion described in Section 3.1. This results in 107
regions encompassing 6,190 square degrees, or 15% of the sky, after masking.
These regions are illustrated in Figure 1 and span roughly two orders of
magnitude in dust column density.

###  3.3 Power Law Fitting

We fit simple power-law models to the $TT$ and $TE$ spectra. Following similar
analyses (e.g., Planck Collaboration XI, 2020), the Galactic signal at each
frequency is modeled as a power law of the form
$A_{d}(\ell/\ell_{0})^{\alpha_{d}}$. We expect $\alpha_{d}\approx-3$ and for
$A_{d}$ to scale with the dust intensity.

For the $TT$ spectra we also include an extragalactic signal that can arise
from the cross-correlation of the CIB seen by both WISE and ACT. We use a
three-component model containing the Galactic component, a clustered
extragalactic component (CIB-C), and a Poisson-like extragalactic component
(CIB-P), with total power spectrum given by

$C_{\ell}^{TT}=A_{d}\,\left(\frac{\ell}{1000}\right)^{\alpha_{d}}+A_{C}\,\left(\frac{\ell}{7000}\right)^{\alpha_{C}}+A_{P}~{}~{}~{}.$
(3)

Here $A_{C}$ is the amplitude of the clustered component, and $\alpha_{C}$ is
its power law index. $A_{P}$ is the scale-independent amplitude of the Poisson
component. We define $A_{\rm CIB}=A_{C}+A_{P}$ as a measure of the total
extragalactic power at $\ell=7000$.

We use the Markov Chain Monte Carlo (MCMC) methods implemented in emcee
(Foreman-Mackey et al., 2013) to estimate parameters, using a Gaussian
likelihood

$-2\ln{\cal
L}=(C_{b}-w_{b\ell}C_{\ell})^{T}Q^{-1}(C_{b}-w_{b\ell}C_{\ell})~{}~{}~{},$ (4)

up to an additive constant. Here $C_{\ell}$ is the model vector, $w_{b\ell}$
are the bandpower window functions that weight the model given the effects of
the mode coupling matrix, $C_{b}$ is the binned data vector, and $Q$ is the
diagonal covariance matrix.

We do not expect a significant extragalactic contribution to the $TE$ spectra,
so we use a simpler two-parameter Galactic model of the form:

$C_{\ell}^{TE}=A_{d}\,\left(\frac{\ell}{1000}\right)^{\alpha_{d}}~{}~{}~{}.$
(5)

Since the signal-to-noise ratio of the $TE$ spectra is lower than for $TT$, we
fix $\alpha_{d}=-2.5$ in the baseline analysis following fit values at lower
multipoles (Planck Collaboration XI, 2020).

##  4 Total Intensity Analysis

Figure 2: WISE W3$\times$ ACT f150 $TT$ spectra of selected regions of varying
column density. The slopes of the best fit power laws are consistent with
previous measurements of Galactic dust power spectra ($\alpha_{d}\sim-3$). In
lower column density regions we see evidence of an extragalactic component.
The best fit model for each region is shown with a dotted line. Figure 3:
Posteriors for the clustered ($A_{C}$) and Poisson ($A_{P}$) components of the
fit to the isotropic high-$\ell$ signal in f150 over six low column density
regions. We marginalize over the Galactic dust emission in each region by
fitting a Galactic dust amplitude ($A_{d}$) in each region assuming
$\alpha_{d}=-3$. A non-zero signal ($A_{C}+A_{P}$) is detected at $30\sigma$
significance. Table 1: Isotropic Background Parameters

|  W3$\times$f150 |  W3$\times$f220
---|---|---
| [$10^{-12}$ (MJy sr-1)2] | [$10^{-12}$ (MJy sr-1)2]
A_C | 0.24±0.04 | 1.2±0.2
A_P | 0.10±0.04 | 0.2±0.2
A_CIB | 0.34±0.01 | 1.44±0.05

Note. — Extragalactic background parameters estimated from WISE $\times$ ACT
$TT$ spectra in six regions of low dust column density. A fixed
$\alpha_{d}=-3$ and $\alpha_{C}=-1$ were assumed (see Equation (3)).

Figure 4: Example model fit of the WISE W3 $\times$ ACT f150 $TT$ spectrum in
a moderately high-latitude region ($l=134.5^{\circ},b=-52.2^{\circ}$; MF14
Tile 236). The top left panel shows the measured $TT$ spectrum (blue circles
with error bars) where $2\sigma$ upper limits are quoted for bandpowers
consistent with zero. Also shown is the total fit model (black) with its
region-specific best-fit Galactic dust component (green) along with the global
best-fit CIB-C and CIB-P components in red and purple, respectively. The
residuals of the fit are presented in the lower left panel. The upper and
lower right panels show the WISE and ACT maps of the region, respectively,
including the applied mask. The $TT$ spectrum transitions from Galactic
emission at low-$\ell$ to extragalactic emission at high-$\ell$ and is well-
fit by the model. Figure 5: Normalized histograms of $\chi^{2}$ values for the
model fits to all 107 regions at f150 (top) and f220 (bottom). The number of
degrees of freedom is indicated with the black dashed line, 20 for f150 and 17
for f220 due to the higher $\ell_{\rm min}$, while the red solid line is the
$\chi^{2}$ probability density function (PDF) for the indicated number of
degrees of freedom. Figure 6: The measured $C_{\ell}^{TT}$ at $\ell=9000$ at
f150 (teal circles) and f220 (orange squares) in each of the 107 regions is
plotted against the fit $A_{d}$. The red dotted line is
$A_{d}\left(9000/1000\right)^{-3}$, corresponding to the Galactic dust term of
Equation (3) for $\alpha_{d}=-3$. The black lines correspond to Equation (3)
with $\alpha_{d}=-3$ and extragalactic parameters at f150 (dashed) and f220
(dash-dotted) from Table 1. The model provides a good description of the data
at both frequencies. Figure 7: The top panel presents a WISE W3 $\times$ ACT
f150 $TT$ spectrum of a region near the Galactic plane ($\ell=193^{\circ}$,
$b=-15^{\circ}$) where the model fit (black dashed) is poor (PTE = 0.06%).
Down arrows indicate 95% upper limits on bandpowers consistent with zero. The
fit residuals are in the bottom panel, demonstrating breakdown of the power
law parameterization at $\ell<3000$. Figure 8: WISE W3 $\times$ ACT f150 $TT$
spectra in two regions (MF14 tiles 170 and 239, located at Galactic latitudes
$b=-24.8^{\circ}$ and $-38.2^{\circ}$, respectively). The best-fit CIB model
has been subtracted from each and the best-fit Galactic dust models are
presented as dashed lines. Notably, the measured power law index of the
Galactic dust $TT$ spectrum ($\alpha_{d}$) differs between these tiles at
$\sim$3$\sigma$ significance ($-2.87\pm 0.07$ versus $-3.12\pm 0.03$). Figure
9: The $\alpha_{d}^{TT}$ slopes and uncertainties estimated from WISE W3
$\times$ ACT $[150,\,220]$ GHz. Regions with $>5\sigma$ detections of Galactic
dust $A_{d}^{TT}$ in both frequencies are opaque. The dashed line shows the
one-to-one line. Histograms show the distribution of $\alpha_{d}$ for the two
frequencies, where $>5\sigma$ measurements of $A_{d}^{TT}$ are opaque. The
dust spectral indices are correlated between f220 and f150.

###  4.1 Power spectra

We detect robust correlation between the WISE and ACT maps over nearly the
entire region of sky analyzed. Specifically, for 106 of 107 regions in our
analysis, we reject the model $C_{\ell}^{TT}=0$ at $>3\sigma$ for both
W3$\times$f150 and W3$\times$f220. A selection of five W3$\times$f150 $TT$
spectra spanning a range of column densities is presented in Figure 2. The
highest column density tiles shown have roughly power-law spectra with
$C_{\ell}\propto\ell^{-3}$ as has been seen for Galactic dust emission (e.g.,
Gautier et al., 1992; Bracco et al., 2011; Hajian et al., 2012; Planck
Collaboration XXX, 2014). While the lower column density regions of Figure 2
are consistent with this behavior at larger scales ($\ell\lesssim 3000$), all
have a rising spectrum in $\mathcal{D}_{\ell}\equiv\ell(\ell+1)C_{\ell}/2\pi$
at higher multipoles. Similar behavior is observed with f220 (not shown).

As an initial validation check of the $TT$ spectra, we perform a null test by
computing the cross-spectra between WISE data in one region with ACT data in
different regions of sky. Specifically, for each of a set of four WISE maps,
we compute the W3$\times$f150 spectrum with ten different ACT regions at the
same declination. As expected, all $TT$ spectra are consistent with zero. For
the 40 $TT$ spectra consisting of 22 bins each, we find $\chi^{2}=838$ for
$C_{\ell}^{TT}=0$ versus 880 degrees of freedom (PTE = 84%).

We therefore seek a physical explanation of the rising $TT$ spectrum at high
multipoles. In the following section, we demonstrate that this signal is
compatible with extragalactic background fluctuations correlated between WISE
and ACT frequencies.

###  4.2 Fitting the Extragalactic Background

Given that a single power law in $\ell$ is an inadequate description of the
$TT$ spectra presented in Figure 2, we consider the model described in
Equation (3) that includes an extragalactic component. Extragalactic emission
should have the same amplitude across the sky while the Galactic dust emission
should vary from region to region. To avoid assuming a constant Galactic dust
$\alpha_{d}$ across all tiles, we proceed in two steps. First, we infer the
extragalactic component parameters $A_{C}$ and $A_{P}$ from a simultaneous fit
to six regions (Tiles 70, 133, 159, 166, 246, 255) of low column density
(median $N_{\rm HI}\lesssim 2.5\times 10^{20}\,{\rm cm}^{-2}$). We then hold
$A_{C}$ and $A_{P}$ fixed to their best fit values to fit $A_{d}$ and
$\alpha_{d}$ in each region separately.

To fit Equation (3) to the six selected regions simultaneously, we first fix
$\alpha_{d}$ and $\alpha_{C}$ to representative values of $-3$ and $-1$,
respectively (see, e.g., Addison et al., 2012). We then fit the six $A_{d}$
parameters, one $A_{C}$ parameter, and one $A_{P}$ parameter, all with
positive definite priors using the methods described in Section 3.3. The
$A_{C}$ and $A_{P}$ posteriors for the f150 fit are presented in Figure 3 with
best fit values from both the f150 and f220 fits listed in Table 1. We have
verified that a simple joint Gaussian likelihood fit to all 107 tiles
simultaneously with $\alpha_{d}=-3$ and $\alpha_{C}=-1$ yields consistent
best-fit values of $A_{C}$ and $A_{P}$.

We find that $A_{\rm CIB}=A_{C}+A_{P}>0$ at 30$\sigma$ significance at both
f150 ($A_{\rm CIB}=3.4\pm 0.1\times 10^{-13}\,{\rm(MJy~{}sr^{-1})^{2}}$) and
f220 ($A_{\rm CIB}=1.44\pm 0.05\times 10^{-12}\,{\rm(MJy~{}sr^{-1})^{2}}$).
The fits yield a ratio $\mathcal{F}_{iso}\equiv A_{\rm CIB}^{\rm f220}/A_{\rm
CIB}^{\rm f150}=4.2\pm 0.2$. The systematic uncertainty on these numbers is of
order 10% (see Section 2.1). Dunkley et al. (2013) performed multi-frequency
fits to the CIB signal in ACT f150 and f220 using a modified blackbody model
$I_{\nu}\propto\nu^{\beta}B_{\nu}(T_{d})$ for the frequency dependence, where
$B_{\nu}(T)$ is the Planck function and $T_{d}$ is the dust temperature. They
found $\beta=2.2\pm 0.1$ for a fixed dust temperature $T_{d}=9.7\,{\rm K}$,
equivalent to $\mathcal{F}_{iso}=4.1\pm 0.2$. Thus, our derived value is
consistent with previous ACT measurements of the CIB at millimeter
wavelengths. We have likewise verified that the fitted Galactic dust $A_{d}$
parameters at f150 and f220 are consistent with a typical frequency scaling
for Galactic dust emission, although, by design, the Galactic dust signal is
weak in these regions and the constraints are not stringent.

The amplitudes $A_{C}$ and $A_{P}$ are strongly anti-correlated, and the
Poisson component is not measured at high significance. To investigate this
further, we perform a fit with just one extragalactic component where we fit
for both amplitude and slope. For the latter, we impose a uniform prior
$[-2,1]$. This fit excludes $\alpha=0$ at $>5\sigma$, i.e., a pure Poisson
component is strongly disfavored. On the other hand, if we fix $\alpha_{C}=-1$
and $A_{P}=0$, we find little degradation in the goodness of fit (PTE = 0.42
versus 0.51 for the fiducial model). In another variation, we fit $A_{P}$,
$A_{C}$, and $\alpha_{C}$, imposing a uniform prior on $\alpha_{C}$ of
$[-1.5,-0.5]$. We find that $\alpha_{C}\leq-0.76$ at 95% confidence. The data
therefore require a component that resembles the clustered component of the
CIB but do not require a Poisson component. However, given the extent of the
degeneracy between $A_{C}$ and $A_{P}$ (see Figure 3), we cannot place strong
constraints on their relative amplitudes.

We have found a high-$\ell$ $TT$ correlation between the WISE and ACT maps
that is well-fit with a single amplitude across six regions and with
amplitudes at f150 and f220 consistent with the frequency scaling of the CIB.
We therefore interpret this signal as correlation between galaxies observed by
both WISE and ACT. Possible origins of the extragalactic component are
discussed further in Section 6.4, but we will first verify that this component
is indeed of constant amplitude across the remaining 101 regions.

For all subsequent $TT$ fits, we fix $A_{C}$ and $A_{P}$ parameters to their
fit values at a given frequency (see Table 1) and fix $\alpha_{C}=-1$.

###  4.3 Fitting the Galactic Dust TT Spectrum

####  4.3.1 Goodness of Fit

In the previous section, we derived the best fit values of the extragalactic
parameters in our $TT$ model ($A_{C}$ and $A_{P}$ in Equation (3)) based on a
set of six regions. In this section, we fix $A_{C}$ and $A_{P}$ to these
values (see Table 1) and perform another MCMC fit separately in each of the
107 regions to derive $A_{d}$ and $\alpha_{d}$. In these fits, $A_{d}$ is
required to be positive and a conservative uniform prior of $[-4,-2]$ is
imposed on $\alpha_{d}$ based on measurements in the literature (Miville-
Deschênes et al., 2016; Planck Collaboration XI, 2020).

An example fit to a WISE W3 $\times$ ACT f150 $TT$ spectrum is presented in
Figure 4. This region, Tile 236, is centered on Galactic longitude
$l=134.5^{\circ}$ and Galactic latitude $b=-52.2^{\circ}$ and transitions from
being dominated by a Galactic dust spectrum at low-$\ell$ to an extragalactic
spectrum at high-$\ell$. The Galactic component is best fit with
$\alpha_{d}=-2.8\pm 0.2$, and overall the parametric model provides an
excellent fit to the data ($\chi^{2}=9.6$ for 20 d.o.f). Overall, the model
provides a good description of the data in all regions. As illustrated in
Figure 5, the distribution of $\chi^{2}$ values across all 107 regions at both
f150 and f220 is broadly consistent with expectations, though some outliers
have high $\chi^{2}$ values.

Another visualization of the model fit to all 107 regions is presented in
Figure 6, which plots the best-fit Galactic dust amplitude $A_{d}$ (see
Equation (3)) against the measured $\ell=9000$ bandpower. In high column
density regions (large $A_{d}$), the high-$\ell$ spectrum is dominated by
Galactic dust emission and there is a strong linear correlation between
$C_{\ell=9000}^{TT}$ and $A_{d}$. As $A_{d}$ decreases to lower column
densities, however, $C_{\ell=9000}^{TT}$ asymptotes to a roughly constant
value in both f150 and f220. This is the extragalactic signal common to all
regions. The sum of the Galactic and extragalactic model components, plotted
in black lines assuming $\alpha_{d}=-3$, provides a good description of the
measurements.

Although the model fits are generally good, Figure 5 demonstrates that the
distribution of $\chi^{2}$ values is biased toward higher values than
expected. Twelve regions in f150 and nine in f220 have fits with PTE $<1$%. We
identify two possible explanations for the model failures: (1) Galactic dust
$TT$ spectra that differ from a pure power law and (2) unmasked compact
Galactic sources. Of the 38 (31) regions with PTE $<\,10\%$ ($\chi^{2}>25$) in
f150 (f220), 18 (14) are regions where departures from the model are mostly at
$\ell<3000$, while the remaining 20 (17) are mostly at $\ell>3000$.

We illustrate an example of (1) in Figure 7, which presents the W3 $\times$
f150 spectrum of a region near the Galactic plane ($b=-15^{\circ}$). While the
$\ell>3000$ spectrum is well-fit by a power law ($\alpha_{d}=-2.93\pm 0.04$),
there are clear departures at lower multipoles. Indeed, the PTE of the fit is
only 0.06%. The relatively high column density of the region permits high
signal-to-noise measurements even at $\ell=10^{4}$, and thus deviations from
our simple parametric model are easier to detect. Further, this region
includes a range of dense, complex Galactic structure likely at different
distances, and so it is not unexpected that the spatial statistics are
complicated. Thus, at least in some regions, we appear to be seeing the
inability of the model to capture the complexity of the dust emission. We
discuss the implications of this further in Section 6.2.

While (1) represents a limitation of the model, (2)—the presence of Galactic
sources in the maps—is a limitation of the analysis. Sources such as stars,
planetary nebulae, and supernova remnants have been identified in ACT maps
(Naess et al., 2020) and many of these have counterparts in the WISE maps.
While we have used maps that are as cleaned of these objects as possible (see
Section 2), residual correlations from objects below our flux cuts could
contribute power particularly at high-$\ell$. Given the robustness of
extragalactic background model fit across a wide range of column densities and
Galactic latitudes (see Figure 6), it is unlikely that such sources constitute
much of the signal we have identified as extragalactic. However, they could
account for regions with high-$\ell$ power in excess of our model. Higher-
fidelity modeling of diffuse dust emission will almost certainly require a
dedicated effort to identify and mask Galactic sources at lower flux
thresholds than employed here.

####  4.3.2 Variation in Galactic dust TT spectrum

We have demonstrated that the model of Equation (3) provides a good
description of the $TT$ spectrum of most of the 107 regions analyzed here.
Although we fit for $\alpha_{d}$ in each region individually, Figure 6
illustrates that a constant $\alpha_{d}=-3$ yields a reasonable fit to the
data. In this section, we demonstrate that there is true $\alpha_{d}$
variability in our sample, justifying our choice of fitting $\alpha_{d}$ as
free parameter and having implications for modeling the dust $TT$ spectrum
more broadly.

Figure 8 illustrates an example of variation in the slope of the W3 $\times$
f150 dust $TT$ spectrum between two regions. MF14 tiles 170 and 239, both at
moderate Galactic latitudes ($b=-24.8^{\circ}$ and $-38.2^{\circ}$,
respectively), have comparable $TT$ power at $\ell\gtrsim 6000$. However, at
lower multipoles they diverge. The best fit $\alpha_{d}$ values for the two
regions are $-3.12\pm 0.03$ and $-2.87\pm 0.07$, respectively. The W3 $\times$
f220 $TT$ spectra of these two regions are best fit by $\alpha_{d}$ values of
$-3.13\pm 0.04$ and $-2.86\pm 0.08$, respectively.

The distribution of best fit $\alpha_{d}$ values is presented in Figure 9,
which shows the best fit $\alpha_{d}$ in f150 against the best fit value at
f220 for 29 regions with $A_{d}/\sigma(A_{d})>5$ at both frequencies. The best
fit values for $\alpha_{d}$ range from $-3.4$ to $-2.7$. In these 29 regions,
we find a median $\alpha_{d}=-2.95$ and $-2.96$ at f150 and f220,
respectively, while the $\alpha_{d}$ values at the two frequencies are
correlated at a level of Pearson $r=0.5$. To assess the statistical
significance of region-to-region variation in $\alpha_{d}$, we consider two
models: (1) $\alpha_{d}$ is constant over all regions, and is estimated as the
inverse variance weighted mean of the 58 $\alpha_{d}$ values in the 29 regions
at the two frequencies; and (2) $\alpha_{d}$ differs from region to region,
and is estimated in each region as the weighted mean of the two $\alpha_{d}$
fits at f150 and f220. Model (1) has $\chi^{2}=151$ for $57$ degrees of
freedom while Model (2) has $\chi^{2}=37$ for $29$ degrees of freedom. The
likelihood ratio test assuming the $\alpha_{d}$ posteriors are Gaussian yields
a $7\sigma$ preference for the model with region-to-region $\alpha_{d}$
variations. The model with constant $\alpha_{d}$ has a PTE of $2\times
10^{-10}$, corresponding to exclusion at $6\sigma$ significance.

The correlation of the fit $\alpha_{d}$ values between f150 and f220 is
consistent with true astrophysical variations in the value of $\alpha_{d}$.
However, some of the variation could be driven by fitting degeneracy between
the amplitudes of the Galactic component and the extragalactic signal. We find
that $\alpha_{d}$ is negatively correlated with column density, i.e., with
shallower slopes at low column densities. This behavior is as predicted by
fitting degeneracy but is not inconsistent with astrophysical variation. The
$5\sigma$ cut employed above mitigates the effect of this fitting degeneracy
on our analysis of $\alpha_{d}$ variations. If these are indeed physical
variations, they could be confirmed with other tracers.

##  5 Polarization Analysis

Figure 10: A selection of the highest signal-to-noise W3 $\times$ f150 $TE$
spectra. The best fit power law $C_{\ell}^{TE}\propto\ell^{-2.5}$ is shown
(red dashed) along with its associated $\chi^{2}$ value. Figure 11: Left: The
best-fit $A_{d}^{TE}$ (see Equation (5)) in each of the 107 regions for each
of f150 and f220 is shown as a function of Galactic latitude. The $A_{d}^{TE}$
have been normalized by the uncertainty of the fit. Vertical red lines connect
selected pairs of f150 and f220 values for the same region. Right: Histogram
of $A_{d}^{TE}/\sigma^{TE}$ with a unit Gaussian overlaid in dark brown,
illustrating bias toward positive $TE$. Figure 12: Distribution of $\chi^{2}$
values across all regions for the fits employing Equation (5). While a simple
power law model $C_{\ell}^{TE}\propto\ell^{-2.5}$ is broadly consistent with
the data, there is evidence for departures. Figure 13: A comparison of the W3
$\times$ f150, f220, and Planck 353 GHz $TE$ spectra of Tile 239 ((l, b)
=(172∘.2, -38∘.2)). The f220 and 353 GHz spectra have been scaled to 150 GHz
assuming a modified blackbody spectral energy distribution with $\beta=1.5$
and $T_{d}=20$ K, corresponding to multiplicative factors of 0.286 and 0.065,
respectively. The Planck 353 GHz spectrum is truncated at $\ell=2000$ due to
lack of Planck sensitivity at higher multipoles. The three spectra are broadly
consistent in amplitude (after scaling) and shape over the full $\ell$ range
shown. Figure 14: Best-fit $A_{d}$ values (see Equation (5)) for $TE$ (x-axis)
and $TB$ (y-axis) in each of the 107 regions at both f150 and f220. Both
quantities have been normalized by the $1\sigma$ uncertainty of the fit with
the shaded regions corresponding to $\pm 3\sigma$. Also shown are the
corresponding histograms with a unit Gaussian overlaid in dark brown. The
distribution of $A^{TE}$ (top) is biased positive (as seen in Figure 11), but
the distribution of $A^{TB}$ (right) shows no significant trend.

In this Section, we analyze the cross power spectra between the W3 total
intensity map and millimeter polarization maps at 150, 220, and 353 GHz,
focusing on the $TE$ spectra. Unlike the $TT$ analysis in Section 4, for the
$TE$ analysis we employ ACT maps that have been co-added with Planck data
(Naess et al., 2020) to enhance signal-to-noise at $\ell\lesssim 2000$ (see
Section 2.1 for details).

We present a sample of nine of the highest signal-to-noise W3$\times{\rm
f}150$ $TE$ spectra in Figure 10. In all nine regions,
$\mathcal{D}_{\ell}^{TE}>0$ over most or all of the $\ell$ range considered
($10^{3}<\ell<10^{4}$).

We first check if systematic temperature-to-polarization leakage in the ACT
maps could bias our measurement of $TE$. Following Lungu et al. (2022), the
expected $T\rightarrow E$ leakage is modeled as a leakage beam
$B^{T\rightarrow E}_{\ell}$ that is determined from ACT $Q$ and $U$ maps of
Uranus. To first order, the expected bias to the $TE$ spectrum is given by
$C_{\ell}^{TT}B^{T\rightarrow E}_{\ell}/B_{\ell}$, where $C_{\ell}^{TT}$ is
the measured WISE $\times$ ACT cross spectrum and $B_{\ell}$ is the
instrumental beam. The leakage signal is found to be largest in f220 and peaks
at the highest multipoles included our analysis ($\ell\sim$4000), but is still
$\lesssim 1\%$ of the measured $TE$ signal. We therefore ignore the systematic
$T\rightarrow E$ leakage in the analysis.

We next apply the data model in Equation (5) to characterize the spectra of
all 107 regions. In our fiducial analysis, we use a uniform prior on $A_{d}$
and fix $\alpha_{d}^{TE}=-2.5$, representative of measurements over large sky
areas at $40<\ell<600$ (Planck Collaboration XI, 2020). We fit over the range
$10^{3}<\ell<10^{4}$ for both f150 and f220. We repeat the analysis with
Planck 353 GHz data, but restrict the fits to $1000<\ell<2000$ given the lack
of constraining power of the Planck polarization data at higher multipoles.

The fits to the ${\it W3}\times{\rm f}150$ $TE$ spectra of the nine selected
regions presented in Figure 10 demonstrate broad but imperfect agreement with
a power law model. In all cases, $A_{d}^{TE}>0$ with
$A_{d}^{TE}/\sigma\left(A_{d}^{TE}\right)>3$, indicating robust detections of
a positive $TE$ signal in all nine regions. Figure 11 shows the
$A_{d}^{TE}/\sigma\left(A_{d}^{TE}\right)$ values of all regions as a function
of Galactic latitude. In total, there are 26 regions with
$A_{d}^{TE}/\sigma\left(A_{d}^{TE}\right)>3$ for f150, 17 regions for f220,
and 13 regions with $A_{d}^{TE}/\sigma\left(A_{d}^{TE}\right)>3$ in both f150
and f220. Most significant detections of non-zero $TE$ signal are in regions
at $|b|<40^{\circ}$ where the dust emission is brightest. There are no regions
in which the $A_{d}^{TE}<0$ at $3\sigma$ significance at either frequency. The
distribution of fit $A_{d}^{TE}$ values is clearly biased toward
$A_{d}^{TE}>0$.

In detail, however, the $TE$ spectra are not all well-described by a power law
in $\ell$. The $\chi^{2}$ distributions of the ensemble of f150 and f220 fits
are presented in Figure 12. Agreement with the expected $\chi^{2}$ probability
density function (PDF) is generally good, but there are more regions with
large $\chi^{2}$ values than expected from chance. Allowing $\alpha_{d}^{TE}$
to vary yields only marginal improvements in most regions, suggesting that the
data have little constraining power on the $TE$ spectral index.

The shape of the $TE$ power spectrum may vary because of the structure of the
diffuse dust emission itself, which is not necessarily well-described by a
power law in $\ell$. On the other hand, the measured $TE$ power spectrum may
be affected by any compact sources in the map that were not identified by the
flags. Further, the WISE data still contain a number of unmitigated data
artifacts that could also affect the spectra particularly at high $\ell$.
Given the limitations imposed by these potential systematics and lack of
evidence for variability in $\alpha_{d}^{TE}$, we do not draw strong
conclusions on the shape of the $TE$ spectra.

Applying our fitting framework to the Planck 353 GHz data yields 19 regions
with $A_{d}^{TE}>0$ at $3\sigma$ confidence. Five of these non-zero $TE$
detections are unique to the 353 GHz analysis while the remaining regions were
identified with f150 (11), f220 (11), or both (8). To assess the impact of the
Planck data in the ACT+Planck co-adds, we repeat the f150 and f220 fits over
the multipole range $2000<\ell<10^{4}$, where the Planck data have little
constraining power. We find that the number of $3\sigma$ positive-$TE$
detections falls from 26 to 17 for f150 and from 17 to 5 for f220.

As further illustration of the relative constraining power of the three
frequencies on the $TE$ spectrum, Figure 13 presents the f150, f220, and
Planck 353 GHz $TE$ spectra of Tile 239, centered on (l, b) =(172∘.2, -38∘.2).
To facilitate direct comparison, we scale the f220 and 353 GHz spectra to 150
GHz assuming a modified blackbody emission law with $\beta=1.5$ and $T_{d}=20$
K. The Planck 353 GHz spectrum is truncated at $\ell=2000$ due to lack of
sensitivity at higher multipoles. There is broad agreement in the amplitude
(after scaling) and shape of the $TE$ spectrum across the three frequencies.
The ACT data are consistent with positive $C_{\ell}^{TE}$ out to scales
$\ell\simeq$6000\. To our knowledge, these are the smallest-scale measurements
of a Galactic dust $TE$ spectrum to date. Using the Capitanio et al. (2017) 3D
reddening map, we estimate that the dust emission toward the center
coordinates of Tile 239 primarily originates from a distance of about 150 pc
from the Sun. At 150 pc, our measurements constrain the dust $TE$ down to
physical scales of $\sim$0.05$~{}\mathrm{pc}$.

The fit $A_{d}^{TE}$ values are highly correlated with the fit $A_{d}^{TT}$
values, as expected. However, we find that the relation between these
quantities is sub-linear—$A_{d}^{TE}$ scales roughly as
$\left(A_{d}^{TT}\right)^{0.8}$. Two physical effects may contribute to this
relation. First, as the dust column density increases, so too does the number
of distinct structures that may be superimposed along the line of sight. While
the $TT$ correlation is unaffected by such superposition, the $TE$ correlation
may be weakened by depolarization and by loss of apparent filamentarity in the
integrated map. Second, at higher column densities a loss of alignment between
the local magnetic field and dust filaments is observed (Planck Collaboration
Int. XXXII, 2016). This should weaken the $TE$ correlation, which is positive
in diffuse regions due to a preferred alignment between dust structures and
the magnetic field (Planck Collaboration Int. XXXVIII, 2016; Clark et al.,
2021).

We repeat this analysis on the $TB$ spectra, finding no 3$\sigma$ detections
of non-zero $TB$ in any region. The fitted $A^{TB}$ amplitudes at both f150
and f220 are compared to the $A^{TE}$ amplitudes in Figure 14. Unlike the fit
$A^{TE}$, the ensemble of $A^{TB}$ amplitudes show no bias toward positive or
negative values at either frequency.

##  6 Discussion

###  6.1 Mid-Infrared PAH Emission as a Spatial Template for Microwave Dust
Emission

Measurements of the CMB at small ($\sim$arcminute) angular scales is a
principal focus of current and upcoming ground-based experiments. Measurements
of lensing of the CMB constrains the growth of structure in the Universe and
the neutrino masses, while removal of the $B$-mode signature generated by
lensing will be required for constraints on primordial $B$-mode signatures at
the levels pursued by next-generation experiments. Critical to all of these
analyses is robustness to Galactic emission at small angular scales. As the
combination of sensitivity and angular resolution do not yet exist to
characterize millimeter-wavelength dust emission at these scales directly,
indirect probes of these dust properties are needed.

We have demonstrated that the WISE W3 measurements of Galactic dust emission
correlate with both total and polarized intensity millimeter-wave observations
at scales $10^{3}<\ell<10^{4}$. Further, the slope of the measured $TT$
spectrum is compatible with $TT$ spectra measured from millimeter data only at
lower multipoles. The WISE data therefore offer a means of characterizing the
spatial structure of Galactic dust emission, including non-Gaussianity, that
will be informative for millimeter-wavelength analyses.

While this is a promising direction for future analysis, we highlight a few
caveats. W3 primarily traces emission from PAHs, which are known to differ in
both emission physics and spatial distribution from the submicron grains
responsible for the bulk of the millimeter emission. The fraction of the dust
mass in PAHs is variable throughout the Galaxy, with fewer PAHs per unit dust
mass found in very dense regions, in H II regions, and in the Warm Neutral
Medium (WNM). Because PAHs undergo single-photon heating, their emission
scales with the first power of the energy density of the interstellar
radiation field $U$. In contrast, millimeter-wavelength emission from
submicron grains scales as $U^{1/4+\beta}\simeq U^{1/5.5}$. Thus, to the
extent that the radiation field heating the dust varies across the sky, the
ratio of PAH emission to millimeter dust emission likewise varies. On the
other hand, the availability of two independent probes of $U$ could help
constrain spatial variations in the dust temperature and thus frequency
decorrelation, a key concern for $B$-mode analyses.

###  6.2 Variability of the Dust Power Spectrum

Simulated maps of Galactic emission used in cosmic microwave background
analyses frequently assume that the Galactic dust $B$-mode spectrum is a power
law in $\ell$ (e.g., Hervías-Caimapo et al., 2016; Thorne et al., 2017).
Measurements with the Planck satellite demonstrate this to be a good
approximation over large sky areas for $40<\ell<600$, with
$C_{\ell}^{BB}\propto\ell^{-2.5}$ (Planck Collaboration XI, 2020). In detail,
however, the dust power spectrum is expected to vary in slope across the sky.
For instance, Marchal & Miville-Deschênes (2021) found that the steepness of
the dust $TT$ spectrum is influenced by the filling factor of the WNM on the
line of sight.

We have shown evidence of spatial variability in the steepness of the dust
$TT$ spectrum at arcminute scales. Further, we find that some of the $TT$
spectra are not well-described by power laws at all. The spatial distribution
of Galactic dust is complex, and it is unsurprising that a power law in $\ell$
is an inadequate description of the power spectrum. If this variability indeed
arises from spatial inhomogeneity of interstellar density structures and the
magnetohydrodynamic turbulence that shapes them, then corresponding
variability is expected in polarization, including the dust $TE$ and $BB$
spectra. Detailed constraints on the shape of the dust $TE$ spectrum and its
variations within the Galaxy will be possible with more sensitive polarization
data.

The observed spatial variations in the dust power spectrum underscore the need
for scrutiny of simple power law models for the scale dependence of Galactic
dust emission. Moment-based methods that can account for deviations from power
law behavior (e.g., Chluba et al., 2017; Azzoni et al., 2023; Vacher et al.,
2023) are of particular interest and could be tested on the regions identified
here.

A challenge for both the analysis presented in this work and for dust modeling
in a CMB foregrounds context is the presence of compact sources. In addition
to extragalactic sources, we have identified a number of Galactic sources such
as planetary nebulae having strong emission at both mid-infrared and
millimeter wavelengths. If left unmasked, these sources can strongly affect
the measured power spectra at high $\ell$. While we have employed aggressive
masking (as described in Section 3.2), sources below the ACT flux density cut
may still contribute non-negligible power. Dedicated identification and
characterization of such sources will be the topic of future work.

###  6.3 Dust TE and TB Correlations at Arcminute Scales

Over large regions of sky and on large angular scales, the dust total
intensity is positively correlated with the dust $E$-mode polarization (Planck
Collaboration Int. XXX, 2016; Planck Collaboration XI, 2020). This positive
$TE$ correlation is consistent with a preferential alignment between elongated
dust intensity structures and the plane-of-sky projected magnetic field
orientation traced by polarized dust emission (Planck Collaboration Int.
XXXVIII, 2016). This magnetically-aligned density anisotropy is also seen in H
I, and additionally provides a natural explanation for the observation that
the dust polarization $EE/BB>1$ (Clark et al., 2015). Filament-based models of
Galactic dust polarization that invoke this alignment also show $TE>0$ (Clark
& Hensley, 2019; Hervías-Caimapo & Huffenberger, 2022).

If magnetically aligned ISM filaments source the observed positive $TE$
correlation, this raises several observationally measurable questions: in
particular, whether and how $TE$ correlation changes as a function of scale
and/or Galactic environment. There could be an environmental $TE$ dependence
set by the relative orientation of filaments and magnetic fields in regions
dominated by different physics. Filamentary structures are strongly aligned
with the magnetic field orientation throughout the diffuse ISM (Clark et al.,
2014), but higher-density filaments are closer to being orthogonal to the
projected magnetic field orientation (Planck Collaboration Int. XXXII, 2016;
Planck Collaboration Int. XXXV, 2016; Fissel et al., 2019b). This empirical
result may be related to the mass-to-flux ratio of molecular cloud filaments
(e.g., Seifried et al., 2020).

A dust filament with a polarization structure that corresponds to a
perpendicular plane-of-sky magnetic field orientation would produce a negative
$TE$ correlation (Zaldarriaga, 2001; Huffenberger et al., 2020). We find that
the measured $TE$ correlation is generally stronger at higher column densities
and lower Galactic latitudes, where the dust is brighter and measured with
higher signal-to-noise, but even at low Galactic latitudes we find no robust
detections of negative $TE$ spectra (Figure 11). The data are thus consistent
with a general alignment between the ACT-measured magnetic field and the
density structures seen in PAH emission. The scale dependence of the $TE$
correlation is plausibly related to the physics that couples the dust density
structure to the magnetic field on a particular scale. This work measures dust
$TE$ that is generally biased toward positive values down to sub-parsec
scales.

Planck data also exhibit a non-zero $TB$ correlation over large sky areas and
large angular scales (Planck Collaboration Int. XXX, 2016; Planck
Collaboration XI, 2020). In the filament-based model, non-zero $TB$ is caused
by imperfect alignment between the long axis of a dust filament and the
magnetic field, such that non-zero $TB$ over large sky regions implies that
this misalignment has a preferred handedness (Huffenberger et al., 2020; Clark
et al., 2021; Cukierman et al., 2023). In cross-correlation with the WISE
data, we find no regions with robustly nonzero $TB$, and no preference for one
sign of $TB$ over the distribution of sky regions considered.

###  6.4 The Cosmic PAH Background

We have detected at $30\sigma$ significance a high-$\ell$ correlation between
the W3 map at 12 $\mu$m and the f150 and f220 ACT maps inconsistent with
extrapolation of the Galactic dust power spectrum. The correlation appears
spatially isotropic and well-described as a power law
$C_{\ell}\propto\ell^{-1}$. We conclude that this signal is extragalactic in
origin and most likely to arise from the correlation between PAH emission in
dusty, star-forming galaxies as seen by WISE and the CIB as seen by ACT.

The restframe MIR emission of a dusty star-forming galaxy is dominated by PAH
features that can account for up to $\sim$20% of its total infrared emission
(Smith et al., 2007). The strongest of these is the 7.7 $\micron$ feature
(Tielens, 2008). Even at $z=0$ this feature makes a non-negligible
contribution to the W3 band, and it remains within the W3 band up to
$z\simeq$1\. Using large optical galaxy catalogues, Chiang & Ménard (2019)
demonstrated that the MF14 map is correlated with galaxies in redshift bins up
to $z\sim 2$, consistent with redshifted PAH emission. Detailed modeling of
the extragalactic background light suggests that the 12 $\micron$
extragalactic sky is dominated by PAH emission from star-forming galaxies
rather than by emission associated with active galactic nuclei (AGN, e.g.,
Andrews et al., 2018). Likewise, in recent multi-wavelength fits to galaxies
detected by ACT, Kilerci et al. (2023) found that even galaxies dominated by
AGN emission at ACT frequencies could be dominated by PAH emission at 12
$\micron$. Thus, the “Cosmic PAH Background” appears the most natural
explanation for the observed correlation, though we cannot rule out
significant contributions from a $12\,\micron$ “Cosmic AGN Background” on the
basis of these data alone.

The cross-power spectrum encodes the relationship between the galaxies
producing the MIR and the millimeter wavelength emission. We find the
correlation is inconsistent with a pure Poisson spectrum, but, as illustrated
in Figure 3, the relative contributions of a clustered versus Poisson
component is not well constrained. Interpretation of shape of the power
spectrum would benefit from forward models of the Cosmic PAH Background based
on galaxy simulations.

The frequency spectrum of the emission, both in the MIR and at millimeter
wavelengths, is a window into galaxy properties. Similarly, quantification of
the level of correlation of maps of diffuse extragalactic emission at two
frequencies constrains the diversity of emission spectra and their variability
with galaxy properties and with cosmic time. For instance, with a greater
number of MIR bands, it will be possible to assess whether different PAH
features preferentially arise from galaxies with different properties.
Likewise, the level of correlation between MIR and millimeter wavelength maps
constrains the extent to which it is the same population of galaxies
responsible for the observed emission in both frequency ranges. Ultimately,
the implementation of PAH emission spectra in tools such as SIDES (Béthermin
et al., 2017) and Websky (Stein et al., 2020) could allow these data to place
constraints on the relationship between PAH-bright galaxies and those
responsible for the CIB, including how the PAH luminosity function evolves
with cosmic time.

In addition to extragalactic emission, the observed high-$\ell$ correlation
could also include a contribution from Galactic point sources. Large numbers
of dusty, compact Galactic point sources have been identified in Planck data
(Planck Collaboration XXVIII, 2016), and ACT has observed objects such as
planetary nebulae that are also bright in the W3 band (Naess et al., 2020).
While we see no evident correlation between the high-$\ell$ component we model
as extragalactic emission and Galactic latitude or dust column density,
careful treatment of Galactic contamination will be required to make
quantitative comparisons between the signal observed here and models of
extragalactic emission.

##  7 Conclusions

We have presented a correlation analysis between 12 $\mu$m emission observed
by WISE and both 150 and 220 GHz emission observed by ACT at multipoles
$10^{3}<\ell<10^{4}$. Our principal conclusions are as follows:

  * •

We report a $30\sigma$ detection of a spatially isotropic, high-$\ell$ $TT$
signal that we interpret as a correlation between the CIB at ACT frequencies
and the “Cosmic PAH Background” seen by WISE. The spectrum is well-fit by a
power law $C_{\ell}\propto\ell^{-1}$, consistent with a clustered component.
The fits do not require, but do not exclude, the presence of a Poisson
component.

  * •

The $TT$ spectrum of Galactic dust at $10^{3}<\ell<10^{4}$ is generally well-
fit by a power law in $\ell$ with $C_{\ell}^{TT}\propto\ell^{-2.9}$,
consistent with Galactic dust $TT$ spectra that have been measured at lower
multipoles (e.g., Gautier et al., 1992; Bracco et al., 2011; Hajian et al.,
2012; Planck Collaboration XXX, 2014). However, we find evidence for spatial
variability in the power law index and identify several regions where a power
law is an inadequate description of the $TT$ spectrum at the sensitivity of
the measurements. The strength of the observed correlation suggests that WISE
maps of dust emission can be used to understand the spatial statistics of
millimeter-wavelength dust emission at small angular scales.

  * •

We identify 35 regions with $>3\sigma$ detections of positive $TE$ correlation
and none with $>3\sigma$ detections of negative $TE$ correlation. We further
find that the distribution of all fit $TE$ amplitudes is biased positive. To
our knowledge, these are the highest-$\ell$ measurements of the dust $TE$
correlation to date.

This work showcases the power of high angular resolution observations of dust
emission at MIR wavelengths to understand the astrophysics of dust emission at
millimeter wavelengths. The small-scale dust morphology in the WISE maps may
be representative of what will be observed by next-generation millimeter
experiments, and so characterization of its non-Gaussianity and other
properties is a promising direction for future work.

The Cosmic PAH background provides another window into the evolution of
galaxies with cosmic time, especially the buildup of PAHs. Implementation of
PAH spectra into existing tools to model the CIB and cross-correlating the
Cosmic PAH Background with other tracers of galaxy properties (e.g., H I
emission) will be important for understanding the properties of the galaxies
giving rise this emission. Given the recent JWST detection of the 2175 Å
feature, associated with PAHs, in a $z=6.71$ galaxy (Witstok et al., 2023),
understanding how the Universe becomes enriched with PAHs is all the more
pressing.

The data underlying the analyses in this work are set to improve dramatically
in the near future. The Simons Observatory will soon begin operations in Chile
and will provide maps of millimeter dust emission and polarization with
greater sensitivity and over a wider frequency range than ACT (Ade et al.,
2019; Hensley et al., 2022). The SPHEREx satellite will soon measure the full
sky at $\sim$6$\arcsec$ angular resolution in 102 channels spanning 0.75–5
$\mu$m, each with comparable sensitivity to WISE (Crill et al., 2020). From
the Cosmic PAH Background to the morphology of dust throughout the Milky Way,
these new datasets can be used to extend the investigations presented here.

## Acknowledgments

We thank B. Draine, D. Finkbeiner, J. Greene, A. Goulding, A. Meisner, M.A.
Miville-Deschênes, and D. Spergel for helpful conversations and guidance
throughout the course of this work. Support for ACT was through the U.S.
National Science Foundation through awards AST-0408698, AST- 0965625, and
AST-1440226 for the ACT project, as well as awards PHY-0355328, PHY-0855887
and PHY-1214379. Funding was also provided by Princeton University, the
University of Pennsylvania, and a Canada Foundation for Innovation (CFI) award
to UBC. ACT operated in the Parque Astronómico Atacama in northern Chile under
the auspices of the Agencia Nacional de Investigacion y Desarrollo (ANID). The
development of multichroic detectors and lenses was supported by NASA grants
NNX13AE56G and NNX14AB58G. Detector research at NIST was supported by the NIST
Innovations in Measurement Science program.

We thank the Republic of Chile for hosting ACT in the northern Atacama, and
the local indigenous Licanantay communities whom we follow in observing and
learning from the night sky.

Computing was performed using the Princeton Research Computing resources at
Princeton University. R.C.R. acknowledges support from the Ford Foundation
Predoctoral Fellowship from the National Academy of Sciences, Engineering, and
Medicine. S.E.C. acknowledges support from the National Science Foundation
grant No. AST-2106607. S.K.C. acknowledges support from NSF award AST-2001866.
C.S. acknowledges support from the Agencia Nacional de Investigación y
Desarrollo (ANID) through FONDECYT grant no. 11191125 and BASAL project
FB210003.

## References

  * Addison et al. (2012) Addison, G. E., Dunkley, J., Hajian, A., et al. 2012, ApJ, 752, 120, doi: 10.1088/0004-637X/752/2/120
  * Ade et al. (2019) Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmology Astropart. Phys, 2019, 056, doi: 10.1088/1475-7516/2019/02/056
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127, doi: 10.1093/mnras/stz093
  * Andrews et al. (2018) Andrews, S. K., Driver, S. P., Davies, L. J. M., Lagos, C. d. P., & Robotham, A. S. G. 2018, MNRAS, 474, 898, doi: 10.1093/mnras/stx2843
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Azzoni et al. (2023) Azzoni, S., Alonso, D., Abitbol, M. H., Errard, J., & Krachmalnicoff, N. 2023, J. Cosmology Astropart. Phys, 2023, 035, doi: 10.1088/1475-7516/2023/03/035
  * Béthermin et al. (2017) Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89, doi: 10.1051/0004-6361/201730866
  * BICEP/Keck Collaboration et al. (2023) BICEP/Keck Collaboration, Ade, P. A. R., Ahmed, Z., et al. 2023, ApJ, 945, 72, doi: 10.3847/1538-4357/acb64c
  * Blagrave et al. (2017) Blagrave, K., Martin, P. G., Joncas, G., et al. 2017, ApJ, 834, 126, doi: 10.3847/1538-4357/834/2/126
  * Bracco et al. (2019) Bracco, A., Ghosh, T., Boulanger, F., & Aumont, J. 2019, A&A, 632, A17, doi: 10.1051/0004-6361/201935951
  * Bracco et al. (2011) Bracco, A., Cooray, A., Veneziani, M., et al. 2011, MNRAS, 412, 1151, doi: 10.1111/j.1365-2966.2010.17971.x
  * Capitanio et al. (2017) Capitanio, L., Lallement, R., Vergely, J. L., Elyajouri, M., & Monreal-Ibero, A. 2017, A&A, 606, A65, doi: 10.1051/0004-6361/201730831
  * Chiang & Ménard (2019) Chiang, Y.-K., & Ménard, B. 2019, ApJ, 870, 120, doi: 10.3847/1538-4357/aaf4f6
  * Chluba et al. (2017) Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195, doi: 10.1093/mnras/stx1982
  * Clark & Hensley (2019) Clark, S. E., & Hensley, B. S. 2019, ApJ, 887, 136, doi: 10.3847/1538-4357/ab5803
  * Clark et al. (2015) Clark, S. E., Hill, J. C., Peek, J. E. G., Putman, M. E., & Babler, B. L. 2015, Phys. Rev. Lett., 115, 241302, doi: 10.1103/PhysRevLett.115.241302
  * Clark et al. (2021) Clark, S. E., Kim, C.-G., Hill, J. C., & Hensley, B. S. 2021, ApJ, 919, 53, doi: 10.3847/1538-4357/ac0e35
  * Clark et al. (2019) Clark, S. E., Peek, J. E. G., & Miville-Deschênes, M. A. 2019, ApJ, 874, 171, doi: 10.3847/1538-4357/ab0b3b
  * Clark et al. (2014) Clark, S. E., Peek, J. E. G., & Putman, M. E. 2014, ApJ, 789, 82, doi: 10.1088/0004-637X/789/1/82
  * Crill et al. (2020) Crill, B. P., Werner, M., Akeson, R., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11443, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114430I, doi: 10.1117/12.2567224
  * Crovisier & Dickey (1983) Crovisier, J., & Dickey, J. M. 1983, A&A, 122, 282
  * Cukierman et al. (2023) Cukierman, A. J., Clark, S. E., & Halal, G. 2023, ApJ, 946, 106, doi: 10.3847/1538-4357/acb0c4
  * Cutri et al. (2012) Cutri, R. M., Wright, E. L., Conrow, T., & et al. 2012, VizieR Online Data Catalog, II/311.  https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/
  * Draine & Li (2007) Draine, B. T., & Li, A. 2007, ApJ, 657, 810, doi: 10.1086/511055
  * Draine et al. (2021) Draine, B. T., Li, A., Hensley, B. S., et al. 2021, ApJ, 917, 3, doi: 10.3847/1538-4357/abff51
  * Dunkley et al. (2013) Dunkley, J., Calabrese, E., Sievers, J., et al. 2013, J. Cosmology Astropart. Phys, 2013, 025, doi: 10.1088/1475-7516/2013/07/025
  * Elmegreen & Scalo (2004) Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211, doi: 10.1146/annurev.astro.41.011802.094859
  * Ferrière (2001) Ferrière, K. M. 2001, Reviews of Modern Physics, 73, 1031, doi: 10.1103/RevModPhys.73.1031
  * Fissel et al. (2019a) Fissel, L., Hull, C. L. H., Clark, S. E., et al. 2019a, Astro2020: Decadal Survey on Astronomy and Astrophysics, 2020, 193, doi: 10.48550/arXiv.1903.08757
  * Fissel et al. (2019b) Fissel, L. M., Ade, P. A. R., Angilè, F. E., et al. 2019b, ApJ, 878, 110, doi: 10.3847/1538-4357/ab1eb0
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foreman-Mackey et al. (2019) Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, The Journal of Open Source Software, 4, 1864, doi: 10.21105/joss.01864
  * García-García et al. (2019) García-García, C., Alonso, D., & Bellini, E. 2019, J. Cosmology Astropart. Phys, 2019, 043, doi: 10.1088/1475-7516/2019/11/043
  * Gautier et al. (1992) Gautier, T. N., I., Boulanger, F., Perault, M., & Puget, J. L. 1992, AJ, 103, 1313, doi: 10.1086/116144
  * Hajian et al. (2012) Hajian, A., Viero, M. P., Addison, G., et al. 2012, ApJ, 744, 40, doi: 10.1088/0004-637X/744/1/40
  * Halal et al. (2023) Halal, G., Clark, S. E., Cukierman, A., Beck, D., & Kuo, C.-L. 2023, arXiv e-prints, arXiv:2306.10107.  https://arxiv.org/abs/2306.10107
  * Henderson et al. (2016) Henderson, S. W., Allison, R., Austermann, J., et al. 2016, Journal of Low Temperature Physics, 184, 772, doi: 10.1007/s10909-016-1575-z
  * Hensley et al. (2016) Hensley, B. S., Draine, B. T., & Meisner, A. M. 2016, ApJ, 827, 45, doi: 10.3847/0004-637X/827/1/45
  * Hensley et al. (2022) Hensley, B. S., Murray, C. E., & Dodici, M. 2022, ApJ, 929, 23, doi: 10.3847/1538-4357/ac5cbd
  * Hervías-Caimapo et al. (2016) Hervías-Caimapo, C., Bonaldi, A., & Brown, M. L. 2016, MNRAS, 462, 2063, doi: 10.1093/mnras/stw1787
  * Hervías-Caimapo & Huffenberger (2022) Hervías-Caimapo, C., & Huffenberger, K. M. 2022, ApJ, 928, 65, doi: 10.3847/1538-4357/ac54b2
  * Hivon et al. (2002) Hivon, E., Gorski, K. M., Netterfield, C. B., et al. 2002, The Astrophysical Journal, 567, 2, doi: 10.1086/338126
  * Huffenberger et al. (2020) Huffenberger, K. M., Rotti, A., & Collins, D. C. 2020, ApJ, 899, 31, doi: 10.3847/1538-4357/ab9df9
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ingalls et al. (2011) Ingalls, J. G., Bania, T. M., Boulanger, F., et al. 2011, ApJ, 743, 174, doi: 10.1088/0004-637X/743/2/174
  * Jarrett et al. (2011) Jarrett, T. H., Cohen, M., Masci, F., et al. 2011, ApJ, 735, 112, doi: 10.1088/0004-637X/735/2/112
  * Kalberla & Haud (2019) Kalberla, P. M. W., & Haud, U. 2019, A&A, 627, A112, doi: 10.1051/0004-6361/201834533
  * Kilerci et al. (2023) Kilerci, E., Hashimoto, T., Goto, T., et al. 2023, ApJS, 265, 45, doi: 10.3847/1538-4365/acb8ba
  * Knox (1995) Knox, L. 1995, Phys. Rev. D, 52, 4307, doi: 10.1103/PhysRevD.52.4307
  * Li et al. (2020) Li, Z., Naess, S., Aiola, S., et al. 2020, J. Cosmology Astropart. Phys, 2020, 010, doi: 10.1088/1475-7516/2020/09/010
  * Lungu et al. (2022) Lungu, M., Storer, E. R., Hasselfield, M., et al. 2022, J. Cosmology Astropart. Phys, 2022, 044, doi: 10.1088/1475-7516/2022/05/044
  * Mainzer et al. (2005) Mainzer, A. K., Eisenhardt, P., Wright, E. L., et al. 2005, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5899, UV/Optical/IR Space Telescopes: Innovative Technologies and Concepts II, ed. H. A. MacEwen, 262–273, doi: 10.1117/12.611774
  * Marchal & Miville-Deschênes (2021) Marchal, A., & Miville-Deschênes, M.-A. 2021, ApJ, 908, 186, doi: 10.3847/1538-4357/abd108
  * Marsden et al. (2014) Marsden, D., Gralla, M., Marriage, T. A., et al. 2014, MNRAS, 439, 1556, doi: 10.1093/mnras/stu001
  * Martin et al. (2015) Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153, doi: 10.1088/0004-637X/809/2/153
  * Martin et al. (2010) Martin, P. G., Miville-Deschênes, M. A., Roy, A., et al. 2010, A&A, 518, L105, doi: 10.1051/0004-6361/201014684
  * Mattila et al. (1996) Mattila, K., Lemke, D., Haikala, L. K., et al. 1996, A&A, 315, L353
  * Meisner & Finkbeiner (2014) Meisner, A. M., & Finkbeiner, D. P. 2014, ApJ, 781, 5, doi: 10.1088/0004-637X/781/1/5
  * Miville-Deschênes et al. (2016) Miville-Deschênes, M. A., Duc, P. A., Marleau, F., et al. 2016, A&A, 593, A4, doi: 10.1051/0004-6361/201628503
  * Miville-Deschênes et al. (2003a) Miville-Deschênes, M. A., Joncas, G., Falgarone, E., & Boulanger, F. 2003a, A&A, 411, 109, doi: 10.1051/0004-6361:20031297
  * Miville-Deschênes & Lagache (2005) Miville-Deschênes, M.-A., & Lagache, G. 2005, ApJS, 157, 302, doi: 10.1086/427938
  * Miville-Deschênes et al. (2007) Miville-Deschênes, M. A., Lagache, G., Boulanger, F., & Puget, J. L. 2007, A&A, 469, 595, doi: 10.1051/0004-6361:20066962
  * Miville-Deschênes et al. (2003b) Miville-Deschênes, M. A., Levrier, F., & Falgarone, E. 2003b, ApJ, 593, 831, doi: 10.1086/376603
  * Naess & Louis (2022) Naess, S., & Louis, T. 2022, arXiv e-prints, arXiv:2210.02243, doi: 10.48550/arXiv.2210.02243
  * Naess et al. (2021) Naess, S., Madhavacheril, M., & Hasselfield, M. 2021, Pixell: Rectangular pixel map manipulation and harmonic analysis library.  http://ascl.net/2102.003
  * Naess et al. (2020) Naess, S., Aiola, S., Austermann, J. E., et al. 2020, J. Cosmology Astropart. Phys, 2020, 046, doi: 10.1088/1475-7516/2020/12/046
  * Onaka et al. (1996) Onaka, T., Yamamura, I., Tanabe, T., Roellig, T. L., & Yuen, L. 1996, PASJ, 48, L59, doi: 10.1093/pasj/48.5.L59
  * Pingel et al. (2022) Pingel, N. M., Dempsey, J., McClure-Griffiths, N. M., et al. 2022, PASA, 39, e005, doi: 10.1017/pasa.2021.59
  * Planck Collaboration XXX (2014) Planck Collaboration XXX. 2014, A&A, 571, A30, doi: 10.1051/0004-6361/201322093
  * Planck Collaboration I (2016) Planck Collaboration I. 2016, A&A, 594, A1, doi: 10.1051/0004-6361/201527101
  * Planck Collaboration XXVIII (2016) Planck Collaboration XXVIII. 2016, A&A, 594, A28, doi: 10.1051/0004-6361/201525819
  * Planck Collaboration I (2020) Planck Collaboration I. 2020, A&A, 641, A1, doi: 10.1051/0004-6361/201833880
  * Planck Collaboration III (2020) Planck Collaboration III. 2020, A&A, 641, A3, doi: 10.1051/0004-6361/201832909
  * Planck Collaboration XI (2020) Planck Collaboration XI. 2020, A&A, 641, A11, doi: 10.1051/0004-6361/201832618
  * Planck Collaboration Int. XXIX (2016) Planck Collaboration Int. XXIX. 2016, A&A, 586, A132, doi: 10.1051/0004-6361/201424945
  * Planck Collaboration Int. XXX (2016) Planck Collaboration Int. XXX. 2016, A&A, 586, A133, doi: 10.1051/0004-6361/201425034
  * Planck Collaboration Int. XXXII (2016) Planck Collaboration Int. XXXII. 2016, A&A, 586, A135, doi: 10.1051/0004-6361/201425044
  * Planck Collaboration Int. XXXV (2016) Planck Collaboration Int. XXXV. 2016, A&A, 586, A138, doi: 10.1051/0004-6361/201525896
  * Planck Collaboration Int. XXXVIII (2016) Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141, doi: 10.1051/0004-6361/201526506
  * Planck Collaboration Int. XLVIII (2016) Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109, doi: 10.1051/0004-6361/201629022
  * Planck Collaboration Int. LVII (2020) Planck Collaboration Int. LVII. 2020, A&A, 643, 42, doi: 10.1051/0004-6361/202038073
  * Seifried et al. (2020) Seifried, D., Walch, S., Weis, M., et al. 2020, MNRAS, 497, 4196, doi: 10.1093/mnras/staa2231
  * Smith et al. (2007) Smith, K. M., Zahn, O., & Doré, O. 2007, Phys. Rev. D, 76, 043510, doi: 10.1103/PhysRevD.76.043510
  * Stein et al. (2020) Stein, G., Alvarez, M. A., Bond, J. R., van Engelen, A., & Battaglia, N. 2020, J. Cosmology Astropart. Phys, 2020, 012, doi: 10.1088/1475-7516/2020/10/012
  * Stinebring et al. (2019) Stinebring, D. R., Chatterjee, S., Clark, S. E., et al. 2019, BAAS, 51, 492, doi: 10.48550/arXiv.1903.07370
  * Thorne et al. (2017) Thorne, B., Dunkley, J., Alonso, D., & Næss, S. 2017, MNRAS, 469, 2821, doi: 10.1093/mnras/stx949
  * Thornton et al. (2016) Thornton, R. J., Ade, P. A. R., Aiola, S., et al. 2016, ApJS, 227, 21, doi: 10.3847/1538-4365/227/2/21
  * Tielens (2008) Tielens, A. G. G. M. 2008, ARA&A, 46, 289, doi: 10.1146/annurev.astro.46.060407.145211
  * Tokunaga & Vacca (2005) Tokunaga, A. T., & Vacca, W. D. 2005, PASP, 117, 421, doi: 10.1086/429382
  * Vacher et al. (2023) Vacher, L., Chluba, J., Aumont, J., Rotti, A., & Montier, L. 2023, A&A, 669, A5, doi: 10.1051/0004-6361/202243913
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Witstok et al. (2023) Witstok, J., Shivaei, I., Smit, R., et al. 2023, arXiv e-prints, arXiv:2302.05468, doi: 10.48550/arXiv.2302.05468
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Zaldarriaga (2001) Zaldarriaga, M. 2001, Phys. Rev. D, 64, 103001, doi: 10.1103/PhysRevD.64.103001