Impact of Dust on Spectral Distortion Measurements of the Cosmic Microwave Background

Many physical processes can lead to departures from the blackbody spectrum of cosmic radiation (Chluba 2014; Chluba & Jeong 2014; Tashiro 2014; Hill et al. 2015)—for example, reionization and structure formation, decaying or annihilating particles, dissipation of primordial density fluctuations, cosmic strings, primordial black holes, small-scale magnetic fields, adiabatic cooling of matter, and recombination.

The distortions created between double Compton scattering decoupling at $z\approx {10}^{6}$ and the decoupling of thermalization by Compton scattering at $z\approx {10}^{5}$ are usually characterized by a chemical potential and are thus called μ distortions (Zeldovich & Sunyaev 1969; Sunyaev & Zeldovich 1970; Illarionov & Siuniaev 1975; Burigana et al. 1991; Hu & Silk 1993). After thermalization decoupling, inverse Compton scattering or other similar processes produce what we call y-type distortions (Zeldovich & Sunyaev 1969; Sunyaev & Zeldovich 1972). Energy injection at intermediate redshifts ($1.5\times {10}^{4}\,\lt z\lt 2\times {10}^{5}$) produces a distortion spectrum that is intermediate between y and μ distortions (Chluba & Sunyaev 2012; Chluba & Jeong 2014; Khatri & Sunyaev 2012).

The FIRAS experiment on NASA's Cosmic Background Explorer (Boggess et al. 1992) provided upper limits on the μ and y parameters at values of $| y| \lt 1.5\times {10}^{-5},| \mu | \,\lt 9\,\times {10}^{-5}$. It put constraints on the deviation of the observed spectrum from the blackbody spectrum of $\tfrac{{\rm{\Delta }}{I}_{\nu }}{{I}_{\nu }}\leqslant {10}^{-5}-{10}^{-4}$ (Mather et al. 1994; Fixsen et al. 1996; Fixsen 2009). Although these constraints were obtained with hardware designed nearly 40 years ago, they have never been superseded. With modern technology, the sensitivity to spectral distortions could be improved by a few orders of magnitude.

The Primordial Inflation Explorer (PIXIE) (Kogut et al. 2011, 2016; Næss et al. 2019; Kogut & Fixsen 2020) is a proposed satellite mission that aims to map the absolute intensity and linear polarization (Stokes I, Q, and U parameters) of the CMB. It would use 416 spectral channels from 14.4 GHz to 6 THz to map the whole sky with an angular resolution of 26. Like its predecessor, FIRAS (Mather et al. 1994; Fixsen et al. 1996; Fixsen 2009), PIXIE would use a polarizing Michelson interferometer with a Fourier transform spectrometer, but with 76 times greater sensitivity.

This increased sensitivity would both improve spectral distortion constraints and constrain primordial gravity waves from the inflationary epoch to $r\lt {10}^{-3}$ at more than 5σ based on the CMB anisotropy.

PIXIE is very sensitive to thermal radiation from dust, and will provide important information about the interstellar medium (ISM). However, the dust foreground poses a challenge for the detectability of spectral distortions. The sensitivity of PIXIE to spectral distortions depends on how well the dust emission can be modeled, which in turn depends on both the correctness of the dust model and its degeneracy with y and μ. In the FIRAS era, the dust spectral energy distribution (SED) was represented by a modified blackbody (MBB) (Reach et al. 1995) or a combination of MBBs (Finkbeiner et al. 1999). This was used and improved by the Wilkinson Microwave Anisotropy Probe (WMAP) team (Bennett et al. 2003) and by Planck (Planck Collaboration 2014a, 2016a). Abitbol et al. (2017) studied the overall impact of foregrounds while using this dust model, and mentioned the need for an analysis using a more comprehensive dust model. Chluba et al. (2017), Remazeilles & Chluba (2020), and Rotti & Chluba (2021) considered how using a moment expansion approach together with internal linear combination (ILC) may provide a way to make use of the spatial information for foreground subtraction for extracting average-sky signals and CMB anisotropies, as an alternative to exactly modeling foregrounds. At this time, the various models describing temperature distribution, the abundance and composition of the dust, and the magnetic alignment of aspherical grains are unconstrained. This requires a broader look at the possible size distribution of dust grains beyond the MBB dust model. Our goal is to provide a better estimation of the impact dust can have on the PIXIE spectral distortion measurement by studying the impact that varying the dust composition, size distribution, and interstellar radiation field (ISRF) can have on the model.

Of primary concern is thermal emission from dust, both in the Milky Way and in distant galaxies. The size and composition of dust grains are variable from place to place in our Galaxy, and possibly across cosmic time. We have ideas about the chief constituents of dust, but detailed knowledge is elusive. Nevertheless, the thermal emission from a wide range of possible dust grains is somewhat similar, and it is plausible that the variability of the emission spectrum can be expressed by a few parameters.

The primary goal of this work is to quantify the variability of dust emission in the PIXIE frequency range, and estimate the sensitivity of PIXIE to y and μ distortions after marginalizing over dust emission. In order for this estimate to be meaningful, it is not required that the dust model be correct, merely that it reflect the range of possible variations in the size distribution and composition for interstellar dust. A model with reasonable constituents and a plausible range of grain size distributions, constrained by well-motivated physical priors, provides at least a lower bound on the variation in the emission spectrum.

This work presents three main results. First, we start from the parameters describing the size distribution of the grains of dust, which have been constrained to match the existing extinction law variability by Zelko & Finkbeiner (2020). The goal is to determine the feasibility of dust foreground subtraction for PIXIE. By exploring the space that the size distributions can span while still maintaining what we know about the dust reddening curve and its spatial variation, we have the opportunity to assess the impact of a broad class of dust models on PIXIE sensitivity. We can also explore the effect of having dust grains exposed to different ISRFs.

For each sample from the size distribution parameter space, we calculate the dust emission at 416 frequencies to be observed by PIXIE. Each of these spectra can be thought of as a point in a 416-dimensional vector space. This ensemble of points occupies a low-dimensional subspace, which we can explore by computing the principal components (PCs) in this space. We then use the PCs to perform a Markov Chain Monte Carlo (MCMC) and Fisher information matrix analysis to determine the impact this improved dust modeling will have on the detectability of y and μ distortions. Our foreground model includes synchrotron radiation, cosmic infrared background (CIB), and free–free emission. We also include distortions coming from the measurement of the temperature of the blackbody to the required precision. The effect of using the incorrect model to describe ISM dust is also characterized.

Second, this work explores the result of performing a CIB fit to contributions coming from dust in galaxies at different temperatures. We create a mock CIB SED, from a superposition of MBBs at different temperatures. We compute the biases that result from fitting the CIB with a simple MBB, a smoothed MBB, and a PC analysis (PCA).

Finally, we study the effect of neglecting the absorption of the CMB's monopole component by interstellar dust and by dust in other galaxies on the detectability of spectral distortions. Nashimoto et al. (2020) studied the impact of the "CMB shadow" caused by galactic matter on measurements of CMB polarization and temperature anisotropy. Prompted by their work, we use an MCMC to explore the offset introduced by failing to model this effect both in the ISM dust and in the CIB, and the corresponding deviations in the spectral distortion parameters.

In Section 2 we explain the modeling of spectral distortions, other foregrounds, and the PIXIE mission configuration. Section 3 describes the modeling of ISM dust from the dust properties constrained by Zelko & Finkbeiner (2020) and the PCA method applied to characterize ISM dust emission. In Section 4 we characterize the impact of ISM dust modeling using both the Fisher information matrix method and an MCMC for foreground analysis. The results of using an incorrect dust model are also explored. Section 5 shows the analysis of a mock CIB created from superposition of multiple temperature sources. Section 6 presents the analysis and results from including the CMB shadows given by ISM dust and dust in other galaxies as part of the model. Finally, we conclude in Section 7.

y Distortion. PIXIE can shed light on the history of star formation by looking at the spectral distortions produced during the epoch of reionization. CMB photons inverse Compton scatter off of gas that has been ionized by early stars, and generally acquire a small amount of energy from each scattering. These Compton distortions are parameterized by a parameter, y, proportional to the energy gain per scattering and the probability of scattering. In practice, it is proportional to the ionized gas pressure, integrated along the line of sight.

The intensity contribution is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{y}={I}_{0}\displaystyle \frac{{x}^{4}{e}^{x}}{{\left({e}^{x}-1\right)}^{2}}\left[x\coth \left(\displaystyle \frac{x}{2}\right)-4\right]y,\end{eqnarray} \tag{ 1 }$$

where ${I}_{0}=\tfrac{2h}{{c}^{2}}{\left(\tfrac{{{kT}}_{0}}{h}\right)}^{3}=270$ MJy sr⁻¹ for ${T}_{0}=2.72548$ K, and $x=\tfrac{h\nu }{{{kT}}_{0}}$.

Hill et al. (2015) estimated that reionization, galaxy groups and clusters, and the intracluster medium will create a total signal of $y=1.77\times {10}^{-6}$. Following Abitbol et al. (2017), we take $y=1.7\times {10}^{-6}$ as our fiducial value (Figure 1). The spectrum of the y distortion is negative below some frequency (a deficit of photons that have been upscattered) and positive above it. The zero crossing of Equation (1) is the solution of $x\coth \left(x/2\right)=4$, which can also be written as $\exp (x)=\tfrac{4+x}{4-x}$. Its positive solution is $x\approx 3.83$, which corresponds to a frequency of 217.5 GHz.

μ Distortion. The amplitude of density fluctuations during inflation translates into energy injection in the CMB (Chluba et al. 2019) for redshifts ${10}^{5}\lt z\lt {10}^{7}$, which leads to μ distortions characterized by the chemical potential, with the intensity given by Equation (2). The chemical potential expresses energy that is absorbed or released when the number of particles in a system changes. Here it is given in units of kT₀ and is therefore dimensionless:

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mu }={I}_{0}\displaystyle \frac{{x}^{4}{e}^{x}}{{\left({e}^{x}-1\right)}^{2}}\left[\displaystyle \frac{1}{\beta }-\displaystyle \frac{1}{x}\right]\mu ,\end{eqnarray} \tag{ 2 }$$

with ${I}_{0}=\tfrac{2h}{{c}^{2}}{\left(\tfrac{{{kT}}_{0}}{h}\right)}^{3}\,=\,270$ MJy sr⁻¹ for ${T}_{0}=2.72548$ K, $\beta =2.1923$, and $x=\tfrac{h\nu }{{{kT}}_{0}}$.

The expected value of μ depends on various assumptions, including assumptions about possible energy injection from new physics. Because there is no consensus about what the value should be, we follow Abitbol et al. (2017) and Chluba et al. (2012) in taking $\mu =2\times {10}^{-8}$ as a fiducial value (Figure 1), in the range sometimes considered for standard ΛCDM. This value is well below the limit detected by FIRAS of $9\times {10}^{-5}$, and in the range plausibly accessible to PIXIE. Its zero crossing occurs at $x=\beta$, which corresponds to a frequency of 124.5 GHz.

Blackbody temperature distortion. The current estimate of the CMB blackbody temperature is ${T}_{0}=2.72548\pm 0.00057{\rm{K}}$ (Fixsen 2009), based on the recalibration of the FIRAS data using WMAP and other measurements in the literature. PIXIE will be able to measure the blackbody temperature to a higher precision than that. This will give rise to a deviation ${{\rm{\Delta }}}_{T}=({T}_{\mathrm{CMB}}-{T}_{0})/{T}_{0}$ from the fiducial T₀ value, creating a spectral distortion (Chluba & Jeong 2014)

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{{{\rm{\Delta }}}_{T}}={I}_{0}\displaystyle \frac{{x}^{4}{e}^{x}}{{\left({e}^{x}-1\right)}^{2}}{{\rm{\Delta }}}_{T},\end{eqnarray} \tag{ 3 }$$

where ${I}_{0}=\tfrac{2h}{{c}^{2}}{\left(\tfrac{{{kT}}_{0}}{h}\right)}^{3}\,=\,270$ MJy sr⁻¹ for ${T}_{0}=2.72548$ K, and $x=\tfrac{h\nu }{{{kT}}_{0}}$.

When choosing the fiducial value for the parameter ${\rm{\Delta }}T$, any value that is below the current measurement precision of 570 μK would be reasonable. For ease of comparison with Abitbol et al. (2017), we adopt their value of ${\rm{\Delta }}T=1.2\times {10}^{-4}$ K.

Other spectral distortions . There are other spectral distortions that are not included in this analysis for simplicity, but could be relevant for PIXIE, such as relativistic temperature corrections to thermal Sunyaev–Zeldovich (y) distortions (Sazonov & Sunyaev 1998; Abitbol et al. 2017).

Synchrotron radiation. Relativistic cosmic-ray electrons are deflected by the Galactic magnetic field and emit synchrotron radiation, which is the dominant foreground at low frequency. The Planck Collaboration (2016a) modeled this radiation as a power law with a flattening at low frequencies. Abitbol et al. (2017) allowed for a more general SED without using a template by fitting a power law with logarithmic curvature:

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{{\rm{s}}}={A}_{{\rm{s}}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\alpha }_{{\rm{s}}}}\left(1+\displaystyle \frac{1}{2}{\omega }_{{\rm{s}}}{\mathrm{ln}}^{2}\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)\right),\end{eqnarray} \tag{ 4 }$$

with ${\nu }_{0}=100\,\mathrm{GHz}$. They estimated the three free parameters, ${A}_{{\rm{s}}},{\alpha }_{{\rm{s}}}$, and ${\omega }_{{\rm{s}}}$, by fitting these to the Planck synchrotron spectrum, obtaining values of 288 Jy sr⁻¹, −0.82, and 0.2, respectively. We use these values as our fiducial values as well.

CIB. Dust present in other galaxies emits thermal radiation that is redshifted on its way to us. The emission coming from all galaxies taken together is referred to as the CIB. The Planck Collaboration (2014b) modeled the CIB SED using an MBB. We take the functional form

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB}}={\alpha }_{\mathrm{CIB}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\beta }_{\mathrm{CIB}}}\displaystyle \frac{{\left(\nu /{\nu }_{0}\right)}^{3}}{{e}^{h\nu /k/{T}_{\mathrm{CIB}}}-1}\end{eqnarray} \tag{ 5 }$$

with ${\nu }_{0}=$ 100 GHz; ${\alpha }_{\mathrm{CIB}}$ is related to the parameter ${A}_{\mathrm{CIB}}$ used by Abitbol et al. (2017) via

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CIB}}={A}_{\mathrm{CIB}}{\left(\displaystyle \frac{h{\nu }_{0}}{{{kT}}_{\mathrm{CIB}}}\right)}^{{\beta }_{\mathrm{CIB}}+3}.\end{eqnarray} \tag{ 6 }$$

Using their values of ${\beta }_{\mathrm{CIB}}=0.6$, ${T}_{\mathrm{CIB}}=18.8$ K, and ${A}_{\mathrm{CIB}}=3.46\times {10}^{-1}$ MJy sr⁻¹ (which correspond to Planck's values as well), we obtain a fiducial value for ${\alpha }_{\mathrm{CIB}}$ of $1.779\times {10}^{-3}$ MJy sr⁻¹.

Abitbol et al. (2017) cautioned against extending the model to high frequencies, which for the purposes of this analysis we do up to 6 THz.

Nashimoto et al. (2020) highlighted the importance of modeling the absorption of the CMB's monopole component. In Section 6, we analyze the impact of including this effect in the foreground modeling, and conclude it is necessary. Fortunately, the emissivity function (that multiplies the Planck function) is the same for both the CMB absorption and the dust emission, so the required modification is straightforward:

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB},{\rm{s}}}={\alpha }_{\mathrm{CIB}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\beta }_{\mathrm{CIB}}}\left(\displaystyle \frac{{\left(\nu /{\nu }_{0}\right)}^{3}}{{e}^{\tfrac{h\nu }{{{kT}}_{\mathrm{CIB}}}}-1}-\displaystyle \frac{{\left(\nu /{\nu }_{0}\right)}^{3}}{{e}^{\tfrac{h\nu }{{{kT}}_{0}}}-1}\right).\end{eqnarray} \tag{ 7 }$$

Free–free emission. Collisions between electrons and ions in ionized gas produce bremsstrahlung (thermal free–free) emission. Following Abitbol et al. (2017), we use an approximation of the spectrum derived from Draine (2011). Defining ${\nu }_{\mathrm{ff}}={\nu }_{\mathrm{FF}}{\left({T}_{{\rm{e}}}/{10}^{3}{\rm{K}}\right)}^{3/2}$, with ${T}_{{\rm{e}}}=7000\,{\rm{K}}$ and ${\nu }_{\mathrm{FF}}\,=255.33\,\mathrm{GHz}$, we have

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mathrm{FF}}={A}_{\mathrm{FF}}\left(1+\mathrm{ln}\left[1+{\left(\displaystyle \frac{{\nu }_{\mathrm{ff}}}{\nu }\right)}^{\sqrt{3}/\pi }\right]\right).\end{eqnarray} \tag{ 8 }$$

In our analysis, we let only the parameter ${A}_{\mathrm{FF}}$ vary. We use a fiducial value of 300 Jy sr⁻¹, which Abitbol et al. (2017) obtained by fitting the model to the free–free spectrum from the Planck Collaboration (2016a).

Other foregrounds not included. In this paper we choose to focus on dust, both Galactic and extragalactic, because this already poses a significant challenge for PIXIE. In future work, we could include other foregrounds like CO emission from distant galaxies, emission from spinning dust grains, and intergalactic dust (Imara & Loeb 2016).

See Table 1 for a summary of the fiducial values for the spectral distortions and foreground parameters other than dust.

PIXIE is expected to observe in two different configurations, aimed at polarization measurements and at intensity measurements (Kogut et al. 2011). For our analysis, we focus on intensity measurement, and assume the mission spends 24 months in this configuration.

Kogut & Fixsen (2020) described the current expected PIXIE sensitivity in the spectral distortion measuring mode. ³ There are 416 frequency bins, each spaced 14.4 GHz apart, spanning a frequency range of 14.4 to 5994.1 GHz. The covariance matrix is diagonal. For 12 months of exposure, for a 1 deg² patch in the sky, at 14.4 GHz, the table gives a value of $8.66\,\times {10}^{-24}$ W m⁻² sr⁻¹ Hz⁻¹ deg⁻² = 8.66 × 10⁻⁴ MJy sr⁻¹ deg⁻² (using 1 MJy = ${10}^{6}\times {10}^{-26}$ W m⁻² Hz⁻¹). There are 41,253 deg² in the sky, of which perhaps 70% might be available for analysis (because the foreground levels are low enough to be tractable). They each give us an independent measurement, for $M=0.7\times 4\pi \times {\left(180/\pi \right)}^{2}=28,877$ independent measurements. Assuming spatially uncorrelated noise, we divide the noise in each bin by $\sqrt{M}$. Next, because the starting values assume a 12 month mission and we assume a 24 month mission, we divide by $\sqrt{2}$. For 14.4 GHz, we obtain $3.60\times {10}^{-6}$ MJy sr⁻¹. The obtained sensitivity-versus-frequency function is plotted in Figure 1. It is flat at low frequency and turns up around 1 THz due to the low-pass filters in the optics. The filters are used to reduce the response to zodiacal light, which would otherwise increase the photon noise.

This section shows the procedure we use to calculate the spectral emission of a collection of dust grains given a size distribution.

The emissivity (the power radiated per unit volume per unit frequency per unit solid angle) coming from a collection of grains is defined as

$$\begin{eqnarray}&&{j}_{\nu }=\displaystyle \sum _{i}\int {da}\tfrac{{{dn}}_{i}}{da}{C}_{\mathrm{abs},i}(\nu ,a){B}_{\nu }({T}_{\mathrm{eq}}(i,a))\end{eqnarray} \tag{ 9 }$$

where i is the index for carbonaceous, silicate, and polycyclic aromatic hydrocarbon grains; ${C}_{\mathrm{abs},i}(\nu ,a)$ is the effective extinction cross section as a function of grain radius a; $({{dn}}_{i})/(da)$ is the size distribution for each type of grain; and ${T}_{\mathrm{eq}}(i,a)$ is the equilibrium temperature of the grain of type i and radius a. The spectral intensity ${I}_{\nu }$ is defined as the emissivity integrated along the line of sight s:

$$\begin{eqnarray}&&{I}_{\nu }=\int {j}_{\nu }{ds}=\int \left(\displaystyle \frac{{j}_{\nu }}{{n}_{{\rm{H}}}}\right)(s){n}_{{\rm{H}}}(s){ds}.\end{eqnarray} \tag{ 10 }$$

Since $({n}_{i})/({n}_{{\rm{H}}})$ is assumed to be constant along the line of sight, $({j}_{\nu })/({n}_{{\rm{H}}})$ becomes constant along the line of sight as well. Using ${N}_{{\rm{H}}}=\int {n}_{{\rm{H}}}(s){ds}$ , we obtain

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu } & = & \displaystyle \frac{{j}_{\nu }}{{n}_{{\rm{H}}}}\int {n}_{{\rm{H}}}(s){ds}=\displaystyle \frac{{j}_{\nu }}{{n}_{{\rm{H}}}}{N}_{{\rm{H}}}\\ & = & {N}_{{\rm{H}}}\displaystyle \sum _{i}\int {da}\displaystyle \frac{1}{{n}_{{\rm{H}}}}\displaystyle \frac{{{dn}}_{i}}{da}{C}_{\mathrm{abs},i}(\nu ,a){B}_{\nu }({T}_{\mathrm{eq}}(i,a)).\end{array}\end{eqnarray} \tag{ 11 }$$

Nashimoto et al. (2020) demonstrated the importance of modeling the absorption of the CMB monopole (see Section 6). In this work we take this effect into account, and consider the net dust emission (subtracting the absorbed CMB):

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu } & = & \displaystyle \frac{{j}_{\nu }}{{n}_{{\rm{H}}}}\int {n}_{{\rm{H}}}(s){ds}=\displaystyle \frac{{j}_{\nu }}{{n}_{{\rm{H}}}}{N}_{{\rm{H}}}\\ & = & {N}_{{\rm{H}}}\displaystyle \sum _{i}\int {da}\displaystyle \frac{1}{{n}_{{\rm{H}}}}\displaystyle \frac{{{dn}}_{i}}{da}{C}_{\mathrm{abs},i}(\nu ,a)\left({B}_{\nu }({T}_{\mathrm{eq}}(a))-{B}_{\nu }({T}_{0})\right)\end{array}\end{eqnarray} \tag{ 12 }$$

where T₀ is the CMB temperature.

For each sample point from Zelko & Finkbeiner (2020) with its corresponding size distribution, using the precomputed values of the temperature T for each radius and type of grain, we integrate to calculate the specific intensity at each of the 416 PIXIE frequencies (Figure 1). Thinking of each such spectrum as a point in a 416-dimensional space, we can use PCA (Shlens 2014) to find a low-dimensional subspace that contains most of the variance. In other words, the ensemble of models from Zelko & Finkbeiner (2020) were generated from an 11-parameter size distribution and a range of ISRF heating, but they span a relatively low-dimensional "dust space" of possible emission spectra. To the extent that the μ and y spectral distortions are orthogonal to this dust space, they can be measured by PIXIE.

When computing the PCs, we first scale the dust emissivity by the PIXIE sensitivity, as described in Section 2.3. Then, we subtract the average.

To reconstruct the spectrum from the PCs, one uses the equation

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mathrm{dust}}=\left(\sum _{i}{W}_{i}{{PC}}_{i}+\mathrm{mean}\right)\times {\sigma }_{\nu }^{\mathrm{PIXIE}}.\end{eqnarray} \tag{ 13 }$$

We perform two different analyses: one in which the ISRF is fixed, and another where we multiply the ISRF by a factor ${\chi }_{\mathrm{ISRF}}$ varying between 0.5 and 2, which generates different equilibrium temperatures for the grains, and thus a different SED.

The 20 PCs we obtain can be seen in Figure 2. Table 2 shows the variance explained by each PC. We see that in both cases, the first component explains most of the variance, but the relative importance of it increases when we let the radiation field vary.

Zoom In Zoom Out Reset image size

A reference intensity is created from all the spectral distortions and foregrounds:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{I}_{\nu }={\rm{\Delta }}{I}_{\nu }^{\mathrm{dust}}+{\rm{\Delta }}{I}_{\nu }^{y}+{\rm{\Delta }}{I}_{\nu }^{\mu }+{\rm{\Delta }}{I}_{\nu }^{{{\rm{\Delta }}}_{T}}+{\rm{\Delta }}{I}_{\nu }^{\mathrm{sync}}\\ \,+\ {\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB}}+{\rm{\Delta }}{I}_{\nu }^{\mathrm{FF}}.\end{array}\end{eqnarray} \tag{ 14 }$$

To generate it, for the spectral distortions and foregrounds other than dust we use the fiducial values in Table 1. For the coefficients of the PCs of the dust, see Table 3.

The MCMC likelihood is defined as $\mathrm{ln}{ \mathcal L }=-0.5\,\times {\rm{\Delta }}{\chi }^{2}$, with ${\rm{\Delta }}{\chi }^{2}={\sum }_{\nu }{\left(({\rm{\Delta }}{I}_{\nu }^{\mathrm{total}}-{\rm{\Delta }}{I}_{\nu }^{\mathrm{reference}})/{\sigma }_{\nu }^{\mathrm{PIXIE}}\right)}^{2}$, where ${\rm{\Delta }}{I}_{\nu }^{\mathrm{reference}}$ is the total intensity calculated with the fiducial parameters, as described above, and ${\rm{\Delta }}{I}_{\nu }^{\mathrm{total}}$ is the intensity obtained as we let the MCMC explore different parameter ranges.

No priors are imposed on the parameters of the MCMC, beyond just physical boundaries.

The following information applies to all the MCMC runs in this paper, including the ones in Section 6: the MCMC uses the ptemcee ⁴ Vousden et al. (2016) package, which uses parallel tempering. This allows for a more efficient exploration of the parameter space than something like the Metropolis–Hastings algorithm. We use five temperatures and 300 walkers, running for 50,000 steps. The runs are thinned by a factor of 50 to reduce autocorrelation. To discard the burn-in phase of the chain, only the last 50% of the steps are retained.

The posteriors appear to be reasonably Gaussian for all parameters except for the synchrotron radiation (Figure 3, Table 4). This is in agreement with Abitbol et al. (2017), who obtained non-Gaussian distributions for the case where they did not impose priors on the synchrotron background, and only obtained Gaussian posteriors once the synchrotron priors were in place. Still, at the 1σ level, the synchrotron posteriors are also Gaussian, so this validates our use of the Fisher information matrix.

Zoom In Zoom Out Reset image size

Table 4. The Results of the MCMC for the Configuration Checking that the Posterior Is Gaussian for All Foregrounds and the Four-PC ISM Dust

Parameter Name	y	μ	${{\rm{\Delta }}}_{T}$	WC0	WC1	WC2	WC3	${A}_{{\rm{s}}}$	${\alpha }_{{\rm{s}}}$	${\omega }_{{\rm{s}}}$	${A}_{\mathrm{FF}}$	${\alpha }_{\mathrm{CIB}}$	${\beta }_{\mathrm{CIB}}$	${T}_{\mathrm{CIB}}$
			(K)					(MJy sr−1)			(MJy sr−1)	(MJy sr−1)		(K)
Fiducial Value	1.70e−06	2.00e−08	1.20e−04	8.30e+05	−2.36e+05	−3.63e+03	−6.74e+03	2.88e−04	−8.20e−01	2.00e−01	3.00e−04	1.78e−03	8.60e−01	1.88e+01
Posterior Median	1.70e−06	3.95e−08	1.20e−04	8.29e+05	−2.36e+05	−3.64e+03	−6.74e+03	2.77e−04	−8.03e−01	2.11e−01	3.07e−04	1.78e−03	8.60e−01	1.88e+01
${\sigma }_{+}$	8.41e−09	1.66e−07	6.37e−08	8.41e+02	1.22e+02	2.41e+01	6.48e+00	7.01e−05	1.34e−01	1.45e−01	2.97e−05	8.77e−07	3.74e−04	2.82e−03
${\sigma }_{-}$	8.33e−09	1.69e−07	6.37e−08	8.43e+02	1.14e+02	2.62e+01	6.29e+00	6.77e−05	2.48e−01	1.54e−01	3.35e−05	9.57e−07	3.46e−04	2.91e−03

Note. ${\sigma }_{+}$ represents the difference between the 84th and 50th quantiles of the posterior, and $\sigma -$ the one between the 50th and 16th.

Download table as:  ASCIITypeset image

Assuming a Gaussian likelihood, one can use the Fisher matrix formalism to calculate the parameter uncertainties. Our likelihood takes the form

$$\begin{eqnarray}\begin{array}{rcl}L & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{n/2}| \det C{| }^{1/2}}\\ & & \times \,\exp \left[-\displaystyle \frac{1}{2}\sum _{{ij}}({p}_{i}-{p}_{0,i}){C}_{{ij}}^{-1}({p}_{j}-{p}_{0,j})\right],\end{array}\end{eqnarray} \tag{ 15 }$$

where the 0 subscript indicates the true parameter values. Near the peak of the likelihood, we can make a Taylor series expansion:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}L=\mathrm{log}L({{\boldsymbol{p}}}_{0})\\ \,+\,\displaystyle \frac{1}{2}\displaystyle \sum _{{ij}}({p}_{i}-{p}_{0,i})\displaystyle \frac{{\partial }^{2}\mathrm{log}L}{\partial {p}_{i}\partial {p}_{j}}{\left.\right|}_{{{\boldsymbol{p}}}_{0}}({p}_{j}-{p}_{0,j})+...\,.\end{array}\end{eqnarray} \tag{ 16 }$$

Taking the average over the random independent variable, the second factor becomes

$$\begin{eqnarray}&&{F}_{{ij}}=-\left\langle \displaystyle \frac{{\partial }^{2}\mathrm{log}L}{\partial {p}_{i}\partial {p}_{j}}\right\rangle .\end{eqnarray} \tag{ 17 }$$

The Fisher matrix is related to the PIXIE noise covariance, C_ab , via

$$\begin{eqnarray}&&{F}_{{ij}}=\sum _{a,b}\displaystyle \frac{\partial {\left({\rm{\Delta }}{I}_{\nu }\right)}_{a}}{\partial {p}_{i}}{C}_{{ab}}^{-1}\displaystyle \frac{\partial {\rm{\Delta }}{\left({I}_{\nu }\right)}_{b}}{\partial {p}_{j}}\end{eqnarray} \tag{ 18 }$$

which is called the Fisher information matrix (see for example Verde 2010). $a,b$ are indexing over frequency, and ${p}_{i},{p}_{j}$ over the different parameters of the model, which in our case are the free parameters for the spectral distortions and the foregrounds. Inverting the F_ij matrix and taking the square root of the diagonal gives us the standard deviations expected for each parameter.

To calculate F_ij , we calculate all the derivatives of the SEDs with respect to the free parameters. In our case, they take a simple analytic form, as described in Appendix A.

As the MCMC approach showed us that we can safely use the Fisher information matrix approach, we perform it using different numbers of PCs, to see what effect complicating the dust model has on the detection sensitivity. The results are displayed in Table 5. One of the principal questions that this work is aiming to answer is how much the sensitivity to μ and y degrades as a function of the complexity of the dust model. We observe that most of the PIXIE noise-weighted variance in spectral space is contained in the first two PCs, and adding up to 20 PCs increases the standard deviation for the y and μ distortion parameters by less than one order of magnitude.

Table 5. Fisher Information Matrix Analysis Results

		Number of PCs Used
	MBB	2PC	3PC	4PC	5PC	6PC	7PC	8PC	9PC	10PC
Fixed radiation field
$\sigma y$	4.81e−09	5.85e−09	6.02e−09	6.02e−09	9.31e−09	9.33e−09	9.35e−09	1.02e−08	1.10e−08	1.23e−08
$\sigma \mu$	8.28e−08	1.11e−07	1.15e−07	1.16e−07	1.79e−07	1.81e−07	1.82e−07	1.96e−07	2.42e−07	2.72e−07
Variable radiation field
$\sigma y$	4.81e−09	4.41e−09	5.58e−09	8.89e−09	9.62e−09	1.18e−08	1.20e−08	1.30e−08	1.67e−08	1.78e−08
$\sigma \mu$	8.28e−08	8.36e−08	1.14e−07	1.77e−07	2.05e−07	2.63e−07	2.63e−07	2.83e−07	3.96e−07	4.09e−07
		Number of PCs Used
	11PC	12PC	13PC	14PC	15PC	16PC	17PC	18PC	19PC	20PC
Fixed radiation field
$\sigma y$	1.41e−08	1.46e−08	1.49e−08	2.15e−08	2.28e−08	2.44e−08	2.68e−08	2.77e−08	2.82e−08	2.86e−08
$\sigma \mu$	3.15e−07	3.32e−07	3.38e−07	5.36e−07	5.87e−07	6.14e−07	6.77e−07	7.20e−07	7.29e−07	7.44e−07
Variable radiation field
$\sigma y$	1.82e−08	1.86e−08	2.02e−08	2.48e−08	2.98e−08	3.54e−08	3.57e−08	3.57e−08	3.82e−08	3.91e−08
$\sigma \mu$	4.14e−07	4.30e−07	4.74e−07	5.94e−07	8.08e−07	9.40e−07	9.44e−07	9.45e−07	9.95e−07	1.02e−06

Note. Standard deviations for the y and μ parameters, given the inclusion of various numbers of PCs. The cases for constant and variable radiation fields are shown. All other foregrounds are included in the analysis.

Download table as:  ASCIITypeset image

We explore all combinations of foregrounds and see what impact they have on the standard deviation of y and μ. The results are summarized in Table 6. Adding all foregrounds one by one produces an increase in the expected standard deviation of the y and μ parameters, with the highest σ expected when all foreground parameters need to be constrained at once. When the foregrounds are modeled only one at a time (leaving the others to be constrained by dedicated experiments), it can be seen that dust produces the lowest impact on detectability, followed by the CIB, the free–free emission, and finally the synchrotron radiation, which produces the most impact. In particular, the μ signal could only be detected if it is stronger than the fiducial value of 2 × 10⁻⁸. However, if the synchrotron background is first constrained by a dedicated experiment, we can have a good possibility of detecting it in the range sometimes considered for the standard ΛCDM of 2 × 10⁻⁸, even in the presence of our dust model.

Table 6. Fisher Information Matrix Estimates of Uncertainty in y and μ (Full-sky Average), Allowing Various Combinations of Foreground Parameters to Float (Section 4.4)

Parameter	$y\times {10}^{9}$	$\mu \times {10}^{9}$	${{\rm{\Delta }}}_{T}$	WC0	WC1	WC2	WC3	${A}_{{\rm{s}}}$	${\alpha }_{{\rm{s}}}$	${\omega }_{{\rm{s}}}$	${\alpha }_{\mathrm{CIB}}$	${\beta }_{\mathrm{CIB}}$	${T}_{\mathrm{CIB}}$	${A}_{\mathrm{FF}}$
			(nK)					(MJy sr−1)			(MJy sr−1)		(K)	(MJy sr−1)
Fiducial	1700	20	1.2e+05	8.3e+05	−2.4e+05	−3.6e+03	−6.7e+03	2.9e−04	−8.2e−01	2.0e−01	1.8e−03	8.6e−01	1.9e+01	3.0e−04
σ	0.9	9.8	1.7	1.0e+00	1.0e+00	1.0e+00	1.0e+00	⋯	⋯	⋯	⋯	⋯	⋯	⋯
σ	1.5	23.9	5.4	⋯	⋯	⋯	⋯	3.0e−06	2.0e−03	6.5e−03	⋯	⋯	⋯	⋯
σ	1.1	10.1	1.8	⋯	⋯	⋯	⋯	⋯	⋯	⋯	1.7e−08	6.6e−06	3.9e−05	...
σ	1.0	11.5	2.3	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	2.5e−07
σ	2.0	27.7	7.2	4.7e+00	2.1e+00	2.2e+00	1.5e+00	4.0e−06	9.9e−03	9.9e−03	⋯	⋯	⋯	...
σ	1.5	10.8	2.6	7.6e+02	6.3e+01	1.1e+01	4.1e+00	⋯	⋯	⋯	2.9e−07	7.6e−05	2.4e−03	...
σ	1.7	20.6	5.3	3.0e+00	1.5e+00	1.5e+00	1.2e+00	⋯	⋯	⋯	⋯	⋯	⋯	8.3e−07
σ	2.2	33.9	11.1	⋯	⋯	⋯	⋯	7.3e−06	3.0e−02	2.7e−02	9.5e−08	2.9e−05	1.5e−04	...
σ	1.7	34.1	8.6	⋯	⋯	⋯	⋯	3.8e−05	9.2e−02	5.6e−02	⋯	⋯	⋯	1.3e−05
σ	1.5	20.7	5.3	⋯	⋯	⋯	⋯	⋯	⋯	⋯	2.5e−08	8.6e−06	4.9e−05	9.2e−07
σ	7.4	140.5	58.8	8.8e+02	1.4e+02	2.7e+01	6.2e+00	4.5e−05	2.1e−01	1.9e−01	1.0e−06	4.1e−04	2.9e−03	...
σ	3.1	62.9	17.3	5.3e+00	2.3e+00	2.4e+00	1.5e+00	8.2e−05	2.0e−01	1.2e−01	⋯	⋯	⋯	2.9e−05
σ	2.0	22.0	6.0	7.6e+02	6.3e+01	1.1e+01	4.1e+00	⋯	⋯	⋯	2.9e−07	7.6e−05	2.4e−03	9.5e−07
σ	2.3	44.6	12.9	⋯	⋯	⋯	⋯	6.7e−05	1.7e−01	1.1e−01	1.0e−07	3.1e−05	1.6e−04	2.2e−05
σ	8.9	177.0	69.4	8.8e+02	1.4e+02	2.8e+01	6.6e+00	8.6e−05	2.6e−01	2.0e−01	1.1e−06	4.3e−04	3.0e−03	3.1e−05

Note. In each row, the 1σ marginalized uncertainty is given for floating parameters, and fixed parameters are indicated by ellipses (...). ISM dust is represented by four PCs, with emission spectra drawn from the Zelko & Finkbeiner (2020) posterior (variable radiation field case). Because the spectra are scaled by the PIXIE sensitivity in the PCA, the uncertainty of the PC coefficients (WCn) is 1.0 for the dust-only case, and larger for the other cases. As more foreground parameters are allowed to float, the uncertainty in y and μ grows. When all foregrounds are included (bottom row), ${\sigma }_{\mu }$ is an order of magnitude larger than the case with fixed synchrotron parameters (third from bottom). This suggests that independent constraints on the synchrotron spectrum would substantially enhance the statistical power of PIXIE.

Download table as:  ASCIITypeset image

In order to validate our methods, we compare the results from the MCMC and the Fisher information matrix. We perform both analyses using the same conditions, modeling the foregrounds for dust, synchrotron radiation, CIB, and free–free emission. The results are listed in Tables 4 and 7, where it can be seen that the two methods are in agreement.

The results presented above using the Fisher information matrix and MCMC are useful in determining the impact ISM dust PCs have on ${\sigma }_{y}$ and ${\sigma }_{\mu }$. However, by the nature of the fits, they were testing a situation in which the correct model is reproducing the data. So the question arises, what happens if the wrong model is used to model the ISM dust? An example of that case is as follows. We take as ground truth the ISM dust SED obtained from the size distributions of the dust grains from a sample of the posteriors in Zelko & Finkbeiner (2020). As a model for the fit, an MBB is used. The goal is to see what biases are introduced to the y and μ distortions.

To fit the ISM SED, an MBB with a CMB shadow contribution is used:

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mathrm{dust}\ \mathrm{MBB},{\rm{s}}}={\tau }_{353}{\left(\displaystyle \frac{\nu }{353\mathrm{GHz}}\right)}^{\beta }\left({B}_{\nu }({T}_{{\rm{D}}})-{B}_{\nu }({T}_{0})\right)\end{eqnarray} \tag{ 19 }$$

where ${T}_{{\rm{D}}}$ is the dust temperature, and T₀ is the CMB temperature.

We show the results of performing fits for each of the 45,000 SEDs obtained from the posterior of Zelko & Finkbeiner (2020) in Figure 4. Fits with 20 PCs from the ISM dust posterior were also performed for comparison. It can be seen that the difference in fits is staggering. Of course, the PCs were by definition obtained on this SED space, so including a sufficient number of them should produce good fits. The case of the MBB introduces biases to both y and μ distortions that make detecting their real values impossible. This is due to the fact that at the PIXIE sensitivity levels, an MBB is a very poor fit to the ISM dust SED created from the dust grain size distributions, as shown in Figure 5. No other foregrounds were included in the analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, an overly simplistic dust model like an MBB will bias spectral distortion measurements at a level orders of magnitude greater than the expected sensitivity of PIXIE. This could potentially be worked around by using spatial information when considering the variation of SED parameters using ILC and moment expansion techniques (Chluba et al. 2017; Remazeilles & Chluba 2020; Rotti & Chluba 2021).

To create a reference CIB spectrum, we compute a "broadened" MBB spectrum, i.e., an MBB averaged over a Gaussian distribution of temperatures (Figure 6). This average is computed by drawing 1000 values for ${T}_{\mathrm{CIB}}$ from a Gaussian with mean ${\mu }_{{T}_{\mathrm{CIB}}}=18.8\,{\rm{K}}$ and standard deviation ${\sigma }_{{T}_{\mathrm{CIB}}}$. We then examine the biases in y and μ as a function of ${\sigma }_{{T}_{\mathrm{CIB}}}$ using three different approaches:

• 1.  
  a simple MBB function as given in Equation (7);
• 2.  
  a PCA with a varying number of CIB PCs; and
• 3.  
  a function approximating an MBB function convolved with a Gaussian of the same ${\sigma }_{{T}_{\mathrm{CIB}}}$ as that of the reference broadened CIB SED.

Zoom In Zoom Out Reset image size

The total reference function for the fit is calculated using

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{I}_{\nu }^{\mathrm{reference}} & = & {\rm{\Delta }}{I}_{\nu }^{\mathrm{sync},\ \mathrm{broadened}}+{\rm{\Delta }}{I}_{\nu }^{y}+{\rm{\Delta }}{I}_{\nu }^{\mu }+\\ & & +{\rm{\Delta }}{I}_{\nu }^{{{\rm{\Delta }}}_{T}}+{\rm{\Delta }}{I}_{\nu }^{\mathrm{dust},\mathrm{MBB}}+{\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB}}+{\rm{\Delta }}{I}_{\nu }^{\mathrm{FF}}.\end{array}\end{eqnarray} \tag{ 20 }$$

The goodness of fit is given by ${\chi }^{2}\,={\sum }_{\nu }{\left(({\rm{\Delta }}{I}_{\nu }^{\mathrm{total}}-{\rm{\Delta }}{I}_{\nu }^{\mathrm{reference}})/{\sigma }_{\nu }^{\mathrm{PIXIE}}\right)}^{2}$, where ${\rm{\Delta }}{I}_{\nu }^{\mathrm{total}}$ is the total intensity calculated by summing over all foreground components, with the CIB intensity defined according to each of the three methods used in the fit.

In analogy to the ISM dust PCA above, we can use PCA to find the dimensionality of PIXIE spectral space spanned by broadened MBBs. If a modest number of PCs can describe the CIB to sufficient precision, they could perhaps be marginalized over without losing too much information about y and μ.

To obtain a set of spectra to input into the PCA, we generate 40,000 broadened MBBs using the technique above. Each MBB is broadened by a different ${\sigma }_{{T}_{\mathrm{CIB}}}$ drawn from a log-uniform distribution of 0.1–6 K, then sampled at the 416 PIXIE frequencies.

The intensity data is normalized by the PIXIE sensitivity before the PCA is done. As expected, the PC spectra show oscillatory behavior with PCn having n extrema (Figure 7). PC1 accounts for 99.9% of the variance, and keeping the first five components leaves $\lt {10}^{-10}$ of the variance behind (Table 8).

Zoom In Zoom Out Reset image size

Our third approach to CIB modeling is to explicitly broaden an MBB. There is no analytical expression for the convolution of a Planck function with a Gaussian.

Appendix B shows that the convolution of a function with a Gaussian can be approximated by an expansion in even derivatives multiplied by powers of σ (Equation (B6)).

Applying this result to the MBB function ${\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB},{\rm{s}}}$ (Equation (7)), we obtain the second-order term: ⁵

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB},{\rm{s}},(2)}=\tfrac{{\sigma }_{\mathrm{CIB}}^{2}}{2}{\alpha }_{\mathrm{CIB}}{\tfrac{\nu }{{\nu }_{0}}}^{3+{\beta }_{\mathrm{CIB}}}\\ \quad \times \left(\frac{2\exp \left(\frac{2h\nu }{{k}_{b}{T}_{\mathrm{CIB}}}\right){h}^{2}{\nu }^{2}}{{\left(\exp \left(\frac{h\nu }{{k}_{b}{T}_{\mathrm{CIB}}}\right)-1\right)}^{3}{k}_{b}^{2}{T}_{\mathrm{CIB}}^{4}}\right.\end{array}\end{eqnarray}$$

$$\begin{eqnarray}&&\begin{array}{l}\qquad -\frac{\exp \left(\displaystyle \frac{h\nu }{{k}_{b}{T}_{\mathrm{CIB}}}\right){h}^{2}{\nu }^{2}}{{\left(\exp \left(\frac{h\nu }{{k}_{b}{T}_{\mathrm{CIB}}}\right)-1\right)}^{2}{k}_{b}^{2}{T}_{\mathrm{CIB}}^{4}}\\ \quad \left.-\frac{2\exp \left(\frac{h\nu }{{k}_{b}{T}_{\mathrm{CIB}}}\right)h\nu }{{\left(\exp \left(\frac{h\nu }{{k}_{b}{T}_{\mathrm{CIB}}}\right)-1\right)}^{2}{k}_{b}{T}_{\mathrm{CIB}}^{3}}\right).\end{array}\end{eqnarray} \tag{ 21 }$$

Thus, up to the second-order term, we can write

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB},{\rm{s}},\ \mathrm{convolution}}({\sigma }_{{T}_{\mathrm{CIB}}})={\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB},{\rm{s}}}+{\rm{\Delta }}{I}_{\nu }^{\mathrm{CIB},{\rm{s}},(2)}.\end{eqnarray} \tag{ 22 }$$

More terms could be included in the expansion, but we restrict it to second-order after checking that higher-order terms do not diverge.

Figures 8 and 9 present the results of performing the fit using the three different methods, with and without other foregrounds.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both the PCA method and the convolution method perform better at fitting the broadened CIB than the MBB function does. The ${\chi }^{2}$ and the errors on y and μ are lower. For ${\sigma }_{{T}_{\mathrm{CIB}}}\lt 0.6$ K, the convolution method produces fits that show the errors on the y and μ distortions are still manageable (i.e., well below the fiducial values). At higher ${\sigma }_{{T}_{\mathrm{CIB}}}$, the biases become larger than the fiducial values, indicating that a meaningful measurement will be difficult with PIXIE data alone.

As we add more PCs, we are able to take the goodness of fit ${\chi }^{2}$ below levels of 1 per degree of freedom (for approximately 400 degrees of freedom) even at ${\sigma }_{{T}_{\mathrm{CIB}}}$ of 6K. This is expected since as we increase the number of PCs, we match the model that generated the data better and better. In a more realistic simulation with PIXIE measurement noise added, the ${\chi }^{2}$ would have a floor of 1 per degree of freedom.

We conclude that using a model that did not generate the data will lead to large biases, especially as ${\sigma }_{{T}_{\mathrm{CIB}}}$ increases beyond 0.6 K. Using a simple MBB to represent the CIB is clearly a poor approximation, though a modest number of PCs may be able to represent the CIB adequately.

Chluba et al. (2017) also looked at the effect of averaging foregrounds over the beam, along the line of sight, and across the sky. They proposed modeling the resulting average SEDs using a moment expansion of the original SED model. It would be interesting to know how many moments one would need to include for the modeling of CIB and other foregrounds for a mission like PIXIE, while checking for the impact on the detectability of CMB spectral distortions.