<figure><img src="image_1.png" /></figure>

# Assessing the Cosmic Infrared Background Monopole from Far-infrared to Millimeter Wavelengths

Jean-Marc Casandjian<fnref n="1" /> ID, Jean Ballet<fnref n="1" /> ID, and Isabelle Grenier<fnref n="2" /> ID

Received 2024 January 29; revised 2024 April 19; accepted 2024 April 24; published 2024 July 4

## Abstract

We measured the cosmic infrared background (CIB) monopole for the COBE Diffuse Infrared Background Experiment (DIRBE) and Planck High-Frequency Instrument (HFI) bands with an updated model for the Galactic dust emission. This model includes a dust excess recently observed in 25% of the sky mainly at high latitude compared to the prediction from $N_{HI}$. We correlated observations from COBE/DIRBE and Planck-HFI with this model to extract the zero levels of the sky maps. We corrected for the isotropic interplanetary dust (IPD) emission and calibration gains and obtained CIB values of $1.4 \pm 8.0, 24.5 \pm 3.9, 15.4 \pm 4.9, 6.8 \pm 2.0, 3.2 \pm 0.3, 1.5 \pm 0.1, 0.40 \pm 0.05, 0.11 \pm 0.04, 0.014 \pm 0.027$, and $0.008 \pm 0.012 \text{nW m}^{-2} \text{sr}^{-1}$ at 60, 100, 140, and 240 µm, and 857, 545, 353, 217, 143, and 100 GHz. We compared those numbers with previous direct CIB measurements and extragalactic source counts. We obtain CIB values lower than previous measurements for wavelengths above 140 µm. Below this value, the large uncertainty related to the IPD emission prevents a clear interpretation.

*Unified Astronomy Thesaurus concepts: Interstellar thermal emission (857); Interstellar medium (847); Dust continuum emission (412)*

## 1. Introduction

The cosmic infrared background (CIB) monopole is the infrared radiation originating from faraway unresolved sources. While the near-infrared part is dominated by the redshifted radiation of stars from galaxies, the far-infrared (FIR) component comes from the emission of galaxies where the dust is heated by ultraviolet photons and reemits in the FIR domain. To a lesser extent, it also comes from the emission of dusty tori around the central black holes of active galaxies. Its intensity is about half of all starlight ever emitted; it is a direct link with contributors to the star formation rate density of the Universe. In the text, we refer to the CIB monopole as CIB.

Direct detection of the CIB is challenging since the strong contribution from foreground emission must first be removed. That encompasses the emission arising within the instrument, the Earth and solar system, and the Galaxy.

The CIB was first firmly detected with the COBE Far Infrared Absolute Spectrophotometer (FIRAS) instrument (Puget et al. 1996; Fixsen et al. 1998) at wavelengths from ~200 to ~2000 µm. In 1998, Schlegel et al. (1998) reported a possible CIB detection at 140 and 240 µm followed by the Diffuse Infrared Background Experiment (DIRBE) team (Arendt et al. 1998; Dwek et al. 1998; Hauser et al. 1998; Kelsall et al. 1998), which detected a CIB of (14 ± 3) nW m⁻<fnref n="2" /> sr⁻<fnref n="1" /> and (25 ± 7) nW m⁻<fnref n="2" /> sr⁻<fnref n="1" /> at 140 and 240 µm. In that work they used a composite map of the DIRBE and Infrared Astronomical Satellite (IRAS) maps at 100 µm as a template for the interstellar medium and extracted the CIB at this wavelength from a correlation with $N_{HI}$. They also

accounted for dust temperature at 240 µm with a linear combination of the 100 and 140 µm maps.

In order to lower the uncertainty inherent to the interplanetary dust (IPD) model in DIRBE observations, Finkbeiner et al. (2000) used two methods based on the inclination of the IPD dust related to the ecliptic and on its ecliptic latitude dependence. They studied the DIRBE weekly averaged sky maps and studied the annual variations of two dimensionless variables almost independent of the choice of IPD dust model. They observed a CIB of $28.1 \pm 7.2 \text{nW m}^{-2} \text{sr}^{-1}$ and $24.6 \pm 8.4 \text{nW m}^{-2} \text{sr}^{-1}$ at 60 and 100 µm.

In Lagache et al. (2000) the DIRBE data were fitted against a linear combination of $N_{HI}$ and the $H_\alpha$ intensity obtained by the Wisconsin H-Alpha Mapper. They derived CIB values of $23.4 \pm 6.3, 24.2 \pm 11.6$, and $11 \pm 6.9 \text{nW m}^{-2} \text{sr}^{-1}$ at 100, 140, and 240 µm.

One can also measure the CIB indirectly from the apparent attenuation of TeV distant sources due to $\gamma_{TeV} + \gamma_{CIB} \rightarrow e^+ + e^-$ electron-positron pair production. This opacity from the CIB can be used to derive its intensity, provided the intrinsic spectrum and redshift of the distant source are known. The energy corresponding to the peak of the $\gamma\gamma$ cross section for FIR photons is ~70 TeV.

In Casandjian et al. (2022, hereafter Paper I) we recalculated the zero levels of the Planck High-Frequency Instrument (HFI) detectors outside regions of the lowest $N_{HI}$. That led to lower zero levels compared to previous works based on correlation at the lowest $N_{HI}$ and therefore to a higher Galactic dust column density. The emission of this dust is a foreground to the CIB and must be removed for an accurate determination of this background. In this work we compared DIRBE and Planck-HFI sky maps to the Galactic dust emission calculated from the modified blackbody (MBB) parameters derived in Paper I and extracted the CIB that we compared to previous measurements.

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France</footnote>

<footnote n="2">² Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France; casandjian@cea.fr</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Diagram summarizing the methodology employed in this study to derive the CIB values. Additionally, we outline the procedure used to derive the dust model in Paper I from the Generalized Needlet Internal Linear Combination (GNILC) maps.</figcaption></figure>

Figure 1 provides a summary of the methodology applied in this study.

## 2. Data

### 2.1. DIRBE

We used the observations of the DIRBE instrument at 60, 100, 140, and 240 µm (bands 7, 8, 9, and 10).<fnref n="3" /> DIRBE used bolometers at 140 and 240 µm and less noisy Germanium photoconductors at 60 and 100 µm. The DIRBE beam size is ~0.7° at all frequencies. The instrumental offsets due to electronic pickup or stray radiation were measured several times per orbit by closing a cold shutter located at the prime focus of the photometer.<fnref n="4" /> DIRBE products have the IPD light and point sources removed. The modeling of the thermal emission (and scattered light) from IPD dust was based on the time-dependent semiphysical parametric model of Kelsall et al. (1998). Measurements were calibrated absolutely using Voyager observations and atmospheric models of Jupiter (140 and 240 µm) and Uranus (60 and 100 µm). DIRBE flux densities are derived assuming a constant $vI_ν$ source spectrum in each band. When the spectrum differs, we must apply color corrections to derive true

flux densities. We used the color-correction coefficients given by the COBE/DIRBE collaboration.<fnref n="5" />

We transformed all DIRBE intensity maps as well as other sky maps used in this work into HEALPix standard (Górski et al. 2005) with a resolution parameter $N_{side} = 512$ corresponding to mean pixel spacing of 0.11°. The observations display fluctuations at very small angular scales due to instrumental noise and CIB anisotropy. To reduce those fluctuations we smoothed the sky maps with a Gaussian symmetric beam of $σ = 1°$ for the 60 and 100 µm bands and $σ = 2°$ at 140 and 240 µm.

We show the resulting DIRBE maps in the first column of Figure 2. The second column of the same figure displays the DIRBE maps minus the zero levels measured in this work (see Section 5).

### 2.2. Planck-HFI

We downloaded the Planck-HFI sky maps<fnref n="6" /> at 857, 545, 353, 217, 143, and 100 GHz (350, 550, 849, 1382, 2096, and

<footnote n="3">³ https://lambda.gsfc.nasa.gov/product/cobe/dirbe_zsma_data_get.cfm</footnote>

<footnote n="4">⁴ https://lambda.gsfc.nasa.gov/product/cobe/dirbe_exsup.html</footnote>

<footnote n="5">⁵ https://lambda.gsfc.nasa.gov/data/cobe/dirbe/ancil/colcorr/DIRBE颜色校正表ASC</footnote>

<footnote n="6">⁶ HFI_SkyMap_857-field-Int_2048_R2.02_full.fits, HFI_SkyMap_545-field-Int_2048_R2.02_full.fits, HFI_SkyMap_353-field-IQU_2048_R2.02_full.fits, HFI_SkyMap_217-field-IQU_2048_R2.02_full.fits, HFI_SkyMap_143-field-IQU_2048_R2.02_full.fits, HFI_SkyMap_100-field-IQU_2048_R2.02_full.fits.</footnote>
2998 µm) from the Planck Legacy Archive.<fnref n="7" /> The processing applied to construct the Planck-HFI maps is described in Planck Collaboration et al. (2016a). The Planck maps were corrected for IPD light emission by fitting the emissivity of various components of the model derived by Kelsall et al. (1998). Observations were also corrected for solar dipole resulting from the motion of the solar system with respect to the CMB. We subtracted the CMB calculated by Bobin et al. (2016). This particular choice has no bearing on our study since all CMB publications based on Planck are quasi-identical outside the Galactic plane. The map units are $K_{CMB}$ for frequencies up to 353 GHz and MJy sr⁻¹ above. We converted $K_{CMB}$ to MJy sr⁻¹ with the routine provided by the Planck collaboration.<fnref n="8" /> We selected the 2015 Planck Public Data Release 2 since this release was used to derive the dust parameters $\tau_{353}$, T, and β in Paper I.

Planck was not designed to measure the zero levels of its detectors. In Odegard et al. (2019) the authors describe the careful calibration of the instrumental offsets of HFI detectors with the COBE/FIRAS observations. The FIRAS calibrator was a temperature-controlled external blackbody (with a temperature varying from 2.2 to 20 K) that was moved about 3 days per month into the sky horn of the spectrometer (Mather et al. 1999). This provided a reference spectrum by which the instrument could be calibrated. The instrumental offsets were measured with the external calibrator kept cold, emitting essentially zero radiation. Odegard et al. (2019) subtracted the CMB monopole and dipole, and the IPD light contributions from FIRAS data and applied weights to match the HFI band frequency response. They correlated the resulting FIRAS data with the 2015 release of HFI maps for each HFI band. They observed an offset of 0.045 ± 0.014, 0.050 ± 0.015, 0.035 ± 0.016, 0.029 ± 0.023, 0.001 ± 0.019, and 0.004 ± 0.014 MJy sr⁻¹ at respectively 857, 545, 353, 217, 143, and 100 GHz (Table 5 of Odegard et al. 2019). We added those offsets to the Planck-HFI maps so that the monopole of those maps should only represent the CIB contribution. Since we used the same Planck map as Odegard et al. (2019), these constants correct for any instrumental offsets and for the presence of any monopole in the Planck maps. We did not perform synchrotron or free-free subtraction in this work, as these components are subdominant to thermal dust emission at the wavelengths considered.

Since HFI spectral responses are not uniform, we must account for the intensity variation within a detector bandpass due to the signal spectral profile. We obtained those color-correction coefficients from dedicated software routines provided by the Planck collaboration.<fnref n="9" />

We also smoothed the sky maps with a Gaussian symmetric beam of $\sigma = 1^{\circ}$. We show the resulting HFI maps in the first column of Figure 3. The second column of the same figure displays the HFI maps minus the zero levels measured in this work (see Section 5).

### 3. Modified Blackbody Model for Galactic Dust

The dust emission observed at FIR, submillimeter, and millimeter wavelengths was modeled with an MBB function: $I_{\nu}(\tau_{353}, T, \beta) = \tau_{353} B_{\nu}(T) (\nu/353 \text{ GHz})^{\beta}$ where $B_{\nu}$ is the Planck

function and where the parameters $\tau_{353}$, T, and $\beta$ represent respectively the optical depth at 353 GHz, the temperature, and the opacity spectral index. We used the parameters $\tau_{353}$, T, and $\beta$ derived in Paper I.

In Paper I we measured the zero levels (instrumental offsets plus CIB) of the HFI maps at 857, 545, and 353 GHz and of the 100 µm composite map from COBE/DIRBE and IRAS. We first correlated $I_{857GHz}$ and $N_{H1}$ in a mask corresponding to the region of longitude 0°–180° and of latitude -60° to -15° in which we removed pixels with dark neutral medium (DNM). We extracted an initial zero level of $I_{857GHz}$ from the y-intercept of an affine fit for $N_{H1}$ between 2 and $5 \times 10^{20} \text{ cm}^{-2}$.

We derived the initial zero levels for the 545 and 353 GHz maps from the one at 857 GHz. We correlated $I_{857GHz}$ and $I_{545GHz}$, removed the zero level of $I_{545GHz}$, and then correlated $I_{545GHz}$ and $I_{353GHz}$. To reduce the effect of T and $\beta$ variation, we performed those correlations for pixels for which the MBB parameters T and $\beta$ did not depart by more than half a standard deviation from the sky average values. In both cases we obtained the initial zero levels from the y-intercept of an affine line fit.

Unfortunately, the variations of T and $\beta$ across the sky prevented obtaining a good correlation between $I_{100\mu}$ and $N_{H1}$ or $I_{857GHz}$, and we were not able to derive a precise value for the initial zero level of $I_{100\mu}$. We used as initial zero level for $I_{100\mu}$ a number randomly sampled between -0.1 and 0.2 MJy sr⁻¹.

Once the initial zero levels of the FIR maps were set, we fitted for each pixel of the sky the parameterized MBB function where the parameters $\tau_{353}$, T, and $\beta$ were free to vary. We adjusted the observations at 100 µm, 857 GHz, 545 GHz, and 353 GHz with their initial zero levels subtracted. We corrected the zero levels by measuring the monopoles of the residual maps where the model deduced from the best-fit MBB spectra was subtracted from the observations.

We repeated this procedure 500 times using, instead of the nominal values of all the parameters involved in the procedure, values sampled from normal distributions whose widths were set to the uncertainty on each value. We also sampled the initial zero level for $I_{100\mu}$ from a uniform distribution ranging from -0.1 to 0.2 MJy sr⁻¹.

With these final zero levels established, we performed the MBB fit pixel by pixel once more, allowing $\tau_{353}$, T, and $\beta$ to vary freely.

From the 500 realizations, we obtained a distribution of 500 values for the zero levels, and for each pixel, 500 sets of MBB parameters $\tau_{353}$, T, and $\beta$. We fitted the zero levels distributions with a Gaussian curve and used the centroids and standard deviations as their final values and uncertainties. We displayed in Figure 7 of Paper I and in Figure 4 of this work the average values of $\tau_{353}$, T, and $\beta$ over those realizations.

To compare the MBB model to observations, we accounted for the MBB model intensities for the different spectral responses of the detectors (color correction) using the DIRBE and HFI coefficients mentioned earlier.

In Figure 1 we provide a summary of the methodology employed in Paper I.

### 4. Mask

Since some $^{12}$CO and $^{13}$CO emission lines fall inside various Planck bandpasses (Planck Collaboration et al. 2014), we excluded in this work pixels associated with CO. We obtained

<footnote n="7">⁷ http://pla.esac.esa.int/pla/#home</footnote>

<footnote n="8">⁸ hfi_unit_conversion.pro</footnote>

<footnote n="9">⁹ https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/Unit_conversion_and_Color_correction</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** All-sky Mollweide display of DIRBE FIR intensity maps at 60, 100, 140, and 240 µm (from top to bottom) in MJy sr⁻¹. The first column shows the raw intensity maps. In the second we subtracted the zero levels calculated in this work. The third column represents the dust MBB model, and the fourth is the residuals between the FIR intensity and a fit of the MBB model based on the correlation shown in Figure 5. Gray pixels in those residuals are pixels excluded in the fit of the correlation of Figure 5. The color scale is linear with a maximum for the first three columns corresponding to twice the average of intensity measured outside the gray pixel mask represented in the fourth column. In the fourth column, the residual sky maps have plus or minus the zero levels for extrema of the color scale so that if they were not removed in the maps they would be filled with red color. We smoothed the sky maps with a Gaussian symmetric beam of σ = 1° for the 60 and 100 µm bands and σ = 2° at 140 and 240 µm. We superimposed to the maps a 30° spaced grid.</figcaption></figure>

the spatial distribution of the velocity-integrated CO brightness temperature $W(CO)$ from the Center for Astrophysics compo-site survey (Dame et al. 2001). To ensure that we covered all the CO emission we convolved the CO map with a Gaussian symmetric beam of $\sigma = 3^\circ$ and we masked pixels for which the resulting smoothed $W(CO) > 0.03$ K km s⁻¹. This condition also ensured that the correlation between observations and model predictions is studied away from regions where $N_H > 5 \times 10^{21}$ cm⁻² where a lever-arm effect could affect the y-intercept in case of nonlinearity.

To avoid a potential contamination by free-free emission in the lowest frequency bands we also masked pixels inside the brightest H II regions corresponding to an $H\alpha$ intensity (Haffner et al. 2010) above 20 Rayleighs. We also excluded the Magellanic Clouds.

This mask was applied to all 10 bands. Additionally, we masked pixels with an absolute ecliptic latitude smaller than 15° for the bands at 100 and 140 µm to avoid residual zodiacal light contamination from IPD.

We increased the absolute ecliptic latitude cut to 25° for the 60 µm band since the IPD dust emission that peaks at about 20 µm is stronger at this wavelength. In the same 60 µm band the emission arising from stochastic heating of very small dust grains causes a departure from a simple MBB form extrapolation (Draine & Li 2001). We do not expect a good correlation between the observation and our modeled emission and therefore restricted the correlation to pixels with lower $N_H$ requiring to mask those with $N_H > 7 \times 10^{20}$ cm⁻².

The masks are displayed as insets in Figure 5.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Same as in Figure 2 but for Planck-HFI maps at 857, 545, 353, 217, and 100 GHz. The maps were corrected for an offset calculated by Odegard et al. (2019) and were CMB subtracted.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** All-sky Mollweide display of the MBB parameters $τ_{353}$ (top, in log scale), $T$ (middle), and $β$ (bottom) extracted from a fit to Planck-HFI and IRAS/DIRBE intensity maps obtained in Paper I. We superimposed to the maps a 30° spaced grid.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

## 5. Zero Levels of Sky Maps

We correlated the observations at 60, 100, 140, and 240 µm and 857, 545, 353, 217, 143, and 100 GHz with the emission predicted by the MBB model at those wavelengths outside the masks described above. We display in Figure 5 the correlation plots in all the bands. In Paper I to study similar correlations we first derived and plotted the centroid position in each X-axis bin with a Gaussian curve fit. However, the centroid of a vertical slice may be biased by the potential asymmetric distributions of adjacent slices. Here we measured the centroid positions in slices perpendicular to the correlation axis. The correlation in the scatterplots shown in Figure 5 may be biased due to the presence of pixels with very low intensity that dominate over the larger ones. So we fitted the centroids between the 10th and 99.9th percentiles of the dust MBB intensity distribution. To characterize the correlation we fitted the centroids with a quadratic function. We obtained a very good linear correlation (except at 60 µm for the reason we already mentioned), but we preferred using a quadratic fit to account for a possible small nonlinearity in the correlation. The residuals between the observations and fit of the model intensities are displayed in the fourth column of Figures 2 and 3 in which we masked all pixels not used in the polynomial fit.

For information, in Figure 6 we show the unmasked residuals. We used tighter scale extrema for the Planck-HFI sky maps. In the residuals of Figures 2 and 3 we used as scale limits the zero levels of the maps measured as the y-intercepts of the correlations. We note that in average the residuals are significantly smaller than the zero levels and therefore that the MBB model gives an accurate representation of the observations. In Figure 6 we observe the contribution of the $J=2 \rightarrow 1$ and $J=1 \rightarrow 0$ $^{12}$CO and $^{13}$CO emission lines. We also observe

some large-scale structures in the lowest frequencies. At 60 and 100 µm we observe residuals corresponding to incorrect IPD light modeling. This is also observed in the correlation plot at 60, 100, and 140 µm in Figure 5 where an excess of MBB model compared to observations is seen (positioned horizontally). The contours are spaced in log, and those excesses have no influence on the fit. As expected, the maps at 140 and 240 µm are noisier than the rest of the bands. We also performed a linear fit to the correlation in order to extract the calibration gain that we used when we converted the zero levels into CIB.

## 6. Correction for Isotropic Interplanetary Dust Component

Early measurements of the CIB at 60 and 100 µm prompted interrogations (Finkbeiner et al. 2000; Protheroe & Meyer 2000; Renault et al. 2001; Dwek & Krennrich 2005) regarding their compatibility with the observed γ-ray spectra of Mrk 501 (Aharonian et al. 1999) and Mrk 421 (Aharonian et al. 2002) above 10 TeV. At those very high energies the γ-ray opacity on FIR CIB photons induces an apparent attenuation of the intrinsic spectrum. A sharp upturn in the Mrk 501 spectrum or pileup at the highest energies was investigated and discarded (see Costamante 2013, and references therein).

Based on the observation of significant high-latitude residuals at 25 µm when using the IPD cloud parameters obtained in Kelsall et al. (1998), Gorjian et al. (2000) derived an alternative set of IPD parameters by imposing an additional condition the absence of residuals at 25 µm. Wright (2004) used this model to measure the CIB in the FIR range. They observed a strong reduction of the CIB at 100 µm to a value of $12.5 \pm 5$ nW m$^{-2}$ sr$^{-1}$ relieving the tension between direct and indirect CIB measurements. At 60 µm they obtained
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Correlations between DIRBE/HFI observations and MBB Galactic dust model emission for the 10 bands. We display 10 contours spaced in log between 1 and 10,000 pixel counts. The red points are the centroids of the distribution of slices perpendicular to the correlation axis. We reprojected the corresponding error bars into a component along the observation axis and a component along the MBB one. Except at 60 µm, we shifted the X-axis by the same amount as the y-intercept of the quadratic fit so that a correlation coefficient of 1 corresponds to the bisector of the plot. The value of the y-intercept is materialized by a horizontal dashed-dotted line. The black lines correspond to the quadratic fit of the histogram points between the 10th and the 99.9th percentile of the dust MBB model intensity distribution. The fit boundaries are represented by vertical dashed-dotted lines. We masked data associated with the white pixels of the masks shown in the insets and detailed in Section 4.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Unmasked residuals between the FIR intensity and a fit of the MBB model based on the correlation shown in Figure 5. Those residuals are similar to those displayed in the fourth column of Figures 2 and 3 with tighter scale extrema for Planck-HFI frequencies.</figcaption></figure>

a formally negative CIB value of −8 ± 14 nW m⁻² sr⁻¹. In Dole et al. (2006), similarly to what was done in Renault et al. (2001), the authors corrected the 100 µm CIB value from Lagache et al. (2000) and the 140 and 240 µm ones from Hauser et al. (1998). They redistributed and subtracted the excess observed at 12 and 25 µm at those larger wavelengths assuming that the IPD excess not accounted for in the Lagache et al. (2000) and Hauser et al. (1998) has the same spectrum as the one used in the model. Additionally, they calibrated the 140 and 240 µm data with the FIRAS photometric scale instead of

applying the DIRBE photometric calibration. As a consequence, the CIB at 100, 140, and 240 µm was reduced to 14.4 ± 6.3, 12 ± 6.9, and 12.25 ± 2.5 nW m⁻² sr⁻¹.

We note that while Finkbeiner et al. (2000) do not use the model of Kelsall et al. (1998) to derive the CIB at 60 and 100 µm, footnote 1 of de Jager & Stecker (2002) indicates that the high CIB values they found could also be due to IPD emission contamination.

We used the work of Sano et al. (2020) to correct for IPD mismodeling in our study. In that paper, the authors
<table><caption>**Table 1**
Intermediate Parameters and Calculated CIB Values</caption><thead><tr><th>Band</th><th>Sky Map Zero Levels<br>(MJy sr<sup>-1</sup>)</th><th>Sky Map Zero Levels<br>(nW m<sup>-2</sup> sr<sup>-1</sup>)</th><th>IPD Correction<br>(nW m<sup>-2</sup> sr<sup>-1</sup>)</th><th>Gain</th><th>Color Correction</th><th>CIB<br>(nW m<sup>-2</sup> sr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>60 µm</td><td>0.45 ± 0.04</td><td>22.3 ± 2.2</td><td>21.1 ± 6.1</td><td>...</td><td>1.23</td><td>1.4 ± 8.0</td></tr><tr><td>100 µm</td><td>0.93 ± 0.12</td><td>27.9 ± 3.6</td><td>4.4 ± 1.3</td><td>0.993 ± 0.001</td><td>1.03</td><td>24.5 ± 3.9</td></tr><tr><td>140 µm</td><td>0.83 ± 0.24</td><td>17.9 ± 5.2</td><td>1.5 ± 0.4</td><td>1.125 ± 0.001</td><td>1.06</td><td>15.4 ± 4.9</td></tr><tr><td>240 µm</td><td>0.59 ± 0.16</td><td>7.4 ± 2.1</td><td>0.27 ± 0.08</td><td>1.0485 ± 0.0005</td><td>1.00</td><td>6.8 ± 2.0</td></tr><tr><td>857 GHz</td><td>0.38 ± 0.03</td><td>3.3 ± 0.3</td><td>...</td><td>1.0170 ± 0.0001</td><td>1.00</td><td>3.2 ± 0.3</td></tr><tr><td>545 GHz</td><td>0.28 ± 0.02</td><td>1.5 ± 0.1</td><td>...</td><td>0.9717 ± 0.0002</td><td>0.96</td><td>1.5 ± 0.1</td></tr><tr><td>353 GHz</td><td>0.12 ± 0.02</td><td>0.42 ± 0.05</td><td>...</td><td>0.9963 ± 0.0003</td><td>0.95</td><td>0.40 ± 0.05</td></tr><tr><td>217 GHz</td><td>0.05 ± 0.02</td><td>0.12 ± 0.04</td><td>...</td><td>0.983 ± 0.002</td><td>0.95</td><td>0.11 ± 0.04</td></tr><tr><td>143 GHz</td><td>0.009 ± 0.018</td><td>0.013 ± 0.025</td><td>...</td><td>0.912 ± 0.003</td><td>1.01</td><td>0.014 ± 0.027</td></tr><tr><td>100 GHz</td><td>0.008 ± 0.022</td><td>0.008 ± 0.012</td><td>...</td><td>0.97 ± 0.02</td><td>0.97</td><td>0.008 ± 0.012</td></tr></tbody></table>

investigated a possible contribution from an isotropic IPD component originating from the dust grains supplied by Oort Cloud comets. The Opik–Oort cloud is a reservoir of icy planetesimals lying beyond the heliosphere (Oort 1950). The observation of nearly isotropically distributed long-period comets with very large orbits suggests that this reservoir must have a quasi-spherical symmetry (Morbidelli 2005, and references therein). This contribution of Oort Cloud comet dust grains has not been accounted for by Kelsall et al. (1998) because they tuned their model parameters using fits to seasonal variation of the sky brightness. In Sano et al. (2020) the authors observed the unaccounted IPD emission component in the residuals between the DIRBE observations and the Galactic emission, including the Kelsall et al. (1998) model. They observed for each DIRBE band below 60 µm a systematic variation of the residuals with the solar elongation angle for DIRBE elongation from 64° to 124°. The shape of this variation was not the one expected from a purely isotropic cloud but can be explained by a spheroidal cloud of higher IPD density farther away from the Sun. They fitted the solar elongation angle dependence of the residuals at 12 µm with a polynomial function and derived an average excess of 206.12 nW m⁻² sr⁻¹ (for Model B corresponding to the raw model of Kelsall et al. 1998). We redistributed this value to FIR wavelengths using the color ratios from (Figure 9 bottom of Kelsall et al. 1998). In the fourth column of Table 1 we show the resulting correction attributed to the IPD emission not accounted for in the Kelsall et al. (1998) model.

## 7. Systematic Uncertainties

We estimated the statistical and systematic uncertainties by applying the same method as in Paper I. We repeated 500 times the whole analysis using, instead of nominal parameter values, numbers sampled assuming a normal distribution with the parameter uncertainties as standard deviation. The uncertainties of the final CIB values derive from the ones associated with FIR data, the MBB dust model, the determination of the y-intercept of the quadratic fit, and the IPD isotropic cloud correction.

The uncertainties associated with zero levels of the HFI detectors were evaluated by Odegard et al. (2019) as 0.014, 0.015, 0.016, 0.023, 0.019, and 0.014 MJy sr⁻¹ at respectively 857, 545, 353, 217, 143, and 100 GHz. The systematic uncertainties in the measurements of the DIRBE instrumental offsets are 1.34, 0.81, 5, and 2 nW m⁻² sr⁻¹ at 60, 100, 140, and 240 µm, respectively (Hauser et al. 1998). Since we were unable to measure the calibration gain of DIRBE at 60 µm we

included an uncertainty of 10.4% (Hauser et al. 1998) in the 60 µm calibration gain.

To assess the uncertainty associated with the dust model we could not sample $\tau_{353}$, T, and $\beta$ from a Gaussian distribution, as T and $\beta$ exhibit complex correlations. As already described, in Paper I we replicated the analysis 500 times, randomly sampling the parameters used to derive the MBB model. That resulted in 500 sets of all-sky maps of $\tau_{353}$, T, and $\beta$. We preserved these 500 fits files and used them in this work to assess the uncertainties related to the MBB model of the Galactic dust.

Additionally, seven parameters were used to extract the y-intercept of the correlation shown in Figures 2 and 3. Most parameters are related to the derivation of the masks used in this analysis and are displayed as insets in Figure 5. We used an uncertainty of 10% for the W(CO) cut at 0.03 K km s⁻¹ and for the $H_\alpha$ intensity cut of 20 Rayleighs. We chose a standard deviation of 1° for the absolute ecliptic latitude cuts of 15° and 25° and an uncertainty of $0.5 \times 10^{20} \text{ cm}^{-2}$ for the $N_H$ cut at $7 \times 10^{20} \text{ cm}^{-2}$. We smoothed the sky maps with a Gaussian symmetric beam of $\sigma = 1°$ or $\sigma = 2°$. In both cases we used an associated standard deviation of 0.2°. As mentioned earlier, we fitted the correlation with a quadratic and linear function between the 10th and 99.9th percentile of the dust MBB intensity distribution. We selected an uncertainty of 0.2 percentile for those two values.

Finally we accounted for the uncertainty in the modeling of the IPD emission. DIRBE modeling of the IPD dust emission was based on Kelsall et al. (1998) with associated uncertainties of 27, 6, 2.3, and 0.5 nW m⁻² sr⁻¹ at 60, 100, 140, and 240 µm (Arendt et al. 1998). Those numbers were deduced from sets of different shapes for the cloud density distribution that gave comparably good fits to the seasonal variation of the sky brightness. Since we corrected for misprediction of this model, we did not use those uncertainties but used those of the correction we applied based on Sano et al. (2020) as follows. We redistributed to FIR wavelengths the excess of 206.12 nW m⁻² sr⁻¹ that they found at 12 µm. As pointed out by Sano et al. (2020), an estimation of the uncertainty related to this number can be found by looking at its variation with the solar elongation angle. The variation of their residuals obtained after subtracting from the DIRBE 60 µm intensity their Galactic plus IDP dust emission model (including the isotropic IPD emission) can be seen in Figure 20(g) of Sano et al. (2020). We measured a half variation of 5 nW m⁻² sr⁻¹ between the maximum and minimum values of their residuals. Additionally, Sano et al. (2020) estimated at 1.08 nW m⁻² sr⁻¹ the
<table><caption>**Table 2**
Contribution of Systematics per Band</caption><thead><tr><th rowspan="2">Band</th><th colspan="4">Contribution of Systematics</th></tr><tr><th>Detector Zero Levels (%)</th><th>MBB Model (%)</th><th>Analysis Cuts (%)</th><th>IPD Emission (%)</th></tr></thead><tbody><tr><td>60 µm</td><td>14</td><td>19</td><td>&lt;1</td><td>67</td></tr><tr><td>100 µm</td><td>16</td><td>58</td><td>3</td><td>24</td></tr><tr><td>140 µm</td><td>66</td><td>26</td><td>2</td><td>5</td></tr><tr><td>240 µm</td><td>69</td><td>27</td><td>1</td><td>3</td></tr><tr><td>857 GHz</td><td>31</td><td>68</td><td>1</td><td>0</td></tr><tr><td>545 GHz</td><td>54</td><td>44</td><td>2</td><td>0</td></tr><tr><td>353 GHz</td><td>82</td><td>17</td><td>1</td><td>0</td></tr><tr><td>217 GHz</td><td>95</td><td>4</td><td>1</td><td>0</td></tr><tr><td>143 GHz</td><td>98</td><td>2</td><td>&lt;1</td><td>0</td></tr><tr><td>100 GHz</td><td>98</td><td>1</td><td>1</td><td>0</td></tr></tbody></table>

uncertainty at 60 µm related to the fact that the spectral energy distribution of the isotropic IDP emission may not be the same as the one of the Kelsall et al. (1998) model used for the distribution at higher wavelengths. We therefore obtained an uncertainty of 6.08 nW m⁻² sr⁻¹ at 60 µm that we extrapolated at 100, 140, and 240 µm using the spectral energy distribution of the total IPD model (Figure 9 bottom of Kelsall et al. 1998). The uncertainties at those wavelengths are provided in the fourth column of Table 1.

In Table 2, we have included the distribution per band of the various categories of systematic uncertainties expressed as percentages of the total systematic uncertainty.

## 8. Results and Interpretation

For each band we set the zero level of the map as the y-intercept of the quadratic fit of the correlation between the data and the MBB model. But we did not use directly the y-intercepts shown in Figure 5 and given in the titles of the column 2 plots of Figures 2 and 3. We used instead the 500 realizations derived to estimate the uncertainties. We obtained a distribution of 500 values of the y-intercepts of the quadratic fit. We fitted those distributions with a Gaussian curve and used the centroids and standard deviations as the final zero levels and uncertainties. Those values are provided in the first column of Table 1. We converted them from surface brightness in units of MJy sr⁻¹ to νI<sub>ν</sub> in units of Wm⁻²sr⁻¹ with νI<sub>ν</sub> (nW m⁻² sr⁻¹) = 3000 µm/λ I<sub>ν</sub> (MJy sr⁻¹) and displayed them in the second column of Table 1.

We proceeded similarly for the calibration gains and fitted with a Gaussian curve the 500 slopes of the same correlation. The gains for all wavelengths but the 60 µm one are shown in the fifth column of Table 1. They are slightly different from those of Paper I because of the different methods used to calculate the correlation centroids, as explained earlier. We note that, except at 60 µm, it could not be measured. The gains are only a few percent off from the nominal value for all wavelengths except at 140 µm and 143 GHz. At 140 µm the gain of 1.125 ± 0.001 is slightly outside the quoted uncertainty of 10.6% (Hauser et al. 1998). At 143 GHz the quoted calibration uncertainty is 7% (Planck Collaboration et al. 2016a) when we find a gain of 0.912 ± 0.003.

We then subtracted the IPD isotropic correction from the zero levels (correlation y-intercepts) and divided the result by the calibration gain assuming 1 at 60 µm. We then fitted the

resulting preliminary CIB values with an MBB function and calculated the color-correction coefficients corresponding to this spectral shape for both DIRBE and HFI. We provide this color correction in column 6 of Table 1. We multiplied the preliminary CIB by the color correction to obtain the final CIB values given in column 7 of Table 1. Our CIB can be correctly modeled by an MBB with the following parameters: $τ_{353} = 1.04 \times 10^{-6}$, $T = 40$ K, $β = -0.26$

We plotted our CIB values with their MBB fit in Figure 7 together with those from Wright (2004), Dole et al. (2006), and Odegard et al. (2019). Those works also provided direct CIB measurements for which the IPD model was corrected or for which no correction was needed. The intensity of our Galactic dust emission model has increased compared to models constructed using very low $N_{H1}$ where $N_{H1}$ alone does not trace correctly the dust column density (Paper I). We therefore observe a decrease of the CIB for wavelengths above 200 µm compared to previous measurements. We find an average ratio of 54% at 240 µm compared to the average measurements of Wright (2004) and Dole et al. (2006), and 65%, 75%, and 76% at 857, 545, and 353 GHz compared to Odegard et al. (2019). We verified in the Appendix that this difference is due to the MBB dust model used here. We repeated this analysis using a model based on Planck (Planck Collaboration et al. 2016b) and we observed a good agreement with previous CIB measurements. At 140 µm, we obtained a CIB in agreement with Wright (2004) and Dole et al. (2006).

At 100 µm, our measured value exceeds that reported in Lagache et al. (2000) and Dole et al. (2006). Our initial zero levels for the sky map, derived from correlating the DIRBE intensity map at 100 µm with the MBB model (as depicted in Figure 5), was 24% higher than the value presented in Dole et al. (2006). Their value was determined through a correlation with $N_{H1}$ and did not account for temperature corrections (Figure 3 of Lagache et al. 2000). Furthermore, our correction for the mismodeling of the IPD, based on Sano et al. (2020), was only half of the correction applied by Dole et al. (2006). While it seems that a more substantial correction may be required to better align with CIB models at 100 µm (as discussed in this section), the extensive correction employed by Dole et al. (2006) surpasses the CIB value at 60 µm.

We also plot in Figure 7 the CIB lower limits derived from galaxy number counts. Dole et al. (2006) used a stacking method to extract the FIR flux density of the Spitzer 24 µm sources above 60 µJy at 70 and 160 µm. While the value at 70 µm is difficult to interpret due to the large uncertainty of our measurement at 60 µm, the interpolation of our CIB fit at 160 µm shows that the value of Dole et al. (2006) represents 75% of our CIB measurement.

In Béthermin et al. (2012) the authors extracted the bright Herschel Spectral and Photometric Imaging Receiver (SPIRE) sources in the COSMOS field and applied a stacking method to derive source counts, including fainter SPIRE sources down to ~2 mJy in the three SPIRE bands at 250, 350, and 500 µm. Note that we did not plot here the results when the counts are extrapolated down to zero flux density so the points should serve as lower limits to the CIB. These limits represent 104%, 124%, and 94% of our CIB at 250, 350, and 500 µm. While these values are within our uncertainties, we observe that there is no space for the contribution of sources with a flux density lower than 2 mJy.
<figure><img src="image_10.png" /></figure>

**Figure 7.** This graph shows the CIB measured in this work from far-infrared to millimeter wavelengths together with the best MBB fit (red-filled circles and red line). We also plot the direct CIB measurements reported by Finkbeiner et al. (2000), Lagache et al. (2000), Wright (2004), Dole et al. (2006), and Odegard et al. (2019) with empty symbols and source counts results from Dole et al. (2006); Béthermin et al. (2012); Oteo et al. (2016), and Duivenvoorden et al. (2020) as solid triangles. The black solid line is the model from Khaire & Srianand (2019) and the dashed line is the model from Béthermin et al. (2011). FIRAS data (gray line) were obtained from Hervé Dole (https://www.ias.u-psud.fr/irgalaxies/files/ebl_sed_hdole.tgz).

The poor spatial resolution of single-dish telescopes leads to source blending, and faint clustered sources can bias the measured flux (Béthermin et al. 2017). To reduce this effect Oteo et al. (2016) used ALMA observations in 69 calibrator fields reaching a flux density limit of 0.2 mJy at 1.2 mm. Our CIB measurement interpolated at this wavelength gives a contribution of only 44% for those resolved sources.

Hsu et al. (2016) used the SCUBA-2 camera on the James Clerk Maxwell Telescope in six lensing cluster fields and three blank fields to measure the counts at fluxes fainter than 1 mJy at 450 µm. Since the integrated source brightness is not provided in the original paper of Hsu et al. (2016) we used the number given in Odegard et al. (2019). We find that those sources contribute to 81% of our CIB interpolated at 450 µm.

Duivenvoorden et al. (2020) applied an improved analysis of the Herschel SPIRE maps at 250, 350, and 500 µm in the COSMOS field by stacking and fitting galaxies detected in eight different prior catalogs. The contribution of those sources is higher than our CIB measures by a ratio of 147%, 184%, and 190% at 250, 350, and 500 µm.

So we observe that some source counts are very close to, or even above, our CIB measurements, which is not expected since the source counts should represent a lower limit for the CIB. However, as pointed out among others by Odegard et al. (2019) large uncertainties exist in those source counts due to source confusion for low-resolution experiments or due to cosmic variance for high-resolution small sample size measurements.

We added to Figure 7 the model of Béthermin et al. (2011). It is a parametric model of the evolution of galaxies. The number of sources per logarithm of luminosity is a classical double exponential function (including a continuous evolution with the redshift), and the spectral energy distribution of the galaxies (starburst or normal) is provided by Lagache et al. (2004).

The parameters of the Béthermin et al. (2011) model (13 free parameters and eight calibration parameters) were adjusted with the number counts of Spitzer/MIPS (at 24, 70, and 160 µm), Herschel/SPIRE (at 250, 350, and 500 µm), the AzTEC camera from the James Clerk Maxwell Telescope (at 1.1 mm) and with various monochromatic luminosity functions from IRAS and Spitzer at various redshifts (local, z = 0.6,1,2). Model parameters were also adjusted on the absolute measurements of the CIB derived by Lagache et al. (1999) from FIRAS at 250, 350, and 500 µm.

We also show in Figure 7 the model of Khaire & Srianand (2019) based on the extinction curve of the Large Magellanic Cloud supergiant shell LMC2. Using this extinction curve they derived the combination of global average star formation rate density and the FUV dust attenuation in galaxies from observed luminosity functions at different wavelengths and redshifts. From this combination they calculate the average energy absorbed by the interstellar dust in the UV to NIR wavelengths and then estimate the FIR emission from galaxies using the local galaxy FIR templates and energy conservation arguments (Khaire & Srianand 2015).

We observe none of the models is in good agreement with our CIB measurements, particularly at 100 µm.

## 9. Conclusions

In this work we correlated Planck-HFI maps, calibrated by Odegard et al. (2019) and COBE/DIRBE maps with dust emission maps from Paper I to extract the CIB, which we further corrected for the isotropic IPD emission from Sano et al. (2020).

In Paper I we derived an MBB dust emission model based on maps from COBE/DIRBE and IRAS at 100 µm and from Planck-HFI, for which we recalculated the offsets. This model is more emissive than previous ones based on calibration at
very low $N_{H\ i}$. We observed for example in Paper I a nearly uniform increase of the dust optical depth at 353 GHz of $7.1 \times 10^{-7}$ compared to the data release 2 Planck map (Planck Collaboration et al. 2016b) and an increase of dust emission in all the bands considered here. We correlated COBE/DIRBE and Planck-HFI intensity sky maps with the dust predictions outside regions of large $N_{H}$. We used the work of Odegard et al. (2019) to correct for the instrumental offsets of HFI detectors. We deduced the zero levels of the DIRBE and HFI maps from the y-intercepts of a quadratic fit of the correlations. The intensity maps used here had the zodiacal emission from IPD removed except from a likely isotropic component that we account for using the work of Sano et al. (2020). We obtained CIB values of $1.4 \pm 8.0, 24.5 \pm 3.9, 15.4 \pm 4.9, 6.8 \pm 2.0, 3.2 \pm 0.3, 1.5 \pm 0.1, 0.40 \pm 0.05, 0.11 \pm 0.04, 0.014 \pm 0.027, 0.008 \pm 0.012 \text{nW m}^{-2} \text{sr}^{-1}$ at 60, 100, 140, and 240 µm, and 857, 545, 353, 217, 143, and 100 GHz. An expected consequence of the increased Galactic dust emission is that the CIB values are lower than previous direct measurements for wavelengths above 140 µm. They are in tension with lower limits provided by source counts that were globally in agreement with previous CIB measurements.

## Acknowledgments

The authors would like to thank Hervé Dole for helpful conversations.

Software: NumPy (van der Walt et al. 2011) Matplotlib (Hunter 2007) FILLTEX (Gerosa & Vallisneri 2017) HEALPix (Górski et al. 2005) ROOT (Brun & Rademakers 1997).

# Appendix
## Cosmic Infrared Background Monopole Using Planck MBB Model

We repeated the CIB measurements using a dust model based on the MBB parameters of Planck Collaboration et al. (2016b). Since the uncertainties related to the MBB dust model are similar to those of our model, we used the same uncertainties for the CIB as those we derived in our work. To simplify the procedure we did not run the 500 realizations to extract the centroids of the parameters in this case but used directly the y-intercept of the correlation. Since we found previously a difference from -0.8% to 10% between the two values this has a minor impact on the final results. We plotted the CIB value corrected for the isotropic IPD emission, calibration gains, and color correction in Figure 8. We observe a very good agreement with the works of Wright (2004), Dole et al. (2006), and Odegard et al. (2019). As previously noted, the disparities at the lowest wavelengths are attributed to the different corrections applied to mitigate the shortcomings of the Kelsall et al. (1998) model. In our approach, we employed the correction derived by Sano et al. (2020) instead of redistributing the entire excess observed at 25 µm. We applied an MBB fitting to those CIB values and derived the following results: $\tau_{353} = 1.82 \times 10^{-6}$, $T = 33$ K, $\beta = -0.12$.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** CIB and associated best MBB fit (red filled circles and red line) obtained with a Galactic dust model based on the MBB parameters of Planck Collaboration et al. (2016b) instead of our model from Paper I. For comparison we show the results of Finkbeiner et al. (2000), Lagache et al. (2000), Wright (2004), Dole et al. (2006), and Odegard et al. (2019).</figcaption></figure>
This graph shows the CIB measured in this work from far-infrared to millimeter wavelengths together with the best MBB fit (red filled circles and red line).

## ORCID iDs

Jean-Marc Casandjian https://orcid.org/0000-0003-1404-0524

Jean Ballet https://orcid.org/0000-0002-8784-2977

Isabelle Grenier https://orcid.org/0000-0003-3274-674X

## References

Aharonian, F., Akhperjanian, A., Barrio, J., et al. 2002, A&A, 384, L23

Aharonian, F. A., Akhperjanian, A. G., Barrio, J. A., et al. 1999, A&A, 349, 11

Arendt, R. G., Odegard, N., Weiland, J. L., et al. 1998, ApJ, 508, 74

Béthermin, M., Dole, H., Lagache, G., Le Borgne, D., & Penin, A. 2011, A&A, 529, A4

Béthermin, M., Le Floc'h, E., Ilbert, O., et al. 2012, A&A, 542, A58

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Bobin, J., Sureau, F., & Starck, J. L. 2016, A&A, 591, A50

Brun, R., & Rademakers, F. 1997, NIMPA, 389, 81

Casandjian, J.-M., Ballet, J., Grenier, I., & Remy, Q. 2022, ApJ, 940, 116

Costamante, L. 2013, IJMPD, 22, 1330025

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792

de Jager, O. C., & Stecker, F. W. 2002, ApJ, 566, 738

Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417

Draine, B. T., & Li, A. 2001, ApJ, 551, 807

Duivenvoorden, S., Oliver, S., Béthermin, M., et al. 2020, MNRAS, 491, 1355

Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1998, ApJ, 508, 106

Dwek, E., & Krennrich, F. 2005, ApJ, 618, 657

Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 2000, ApJ, 544, 81

Fixsen, D. J., Dwek, E., Mather, J. C., Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123

Gerosa, D., & Vallisneri, M. 2017, JOSS, 2, 222

Gorjian, V., Wright, E. L., & Chary, R. R. 2000, ApJ, 536, 550

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Haffner, L. M., Reynolds, R. J., Madsen, G. J., et al. 2010, in ASP Conf. Ser. 438, The Dynamic Interstellar Medium: A Celebration of the Canadian Galactic Plane Survey, ed. R. Kothes, T. L. Landecker, & A. G. Willis (San Francisco, CA: ASP), 388

Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25

Hsu, L.-Y., Cowie, L. L., Chen, C.-C., Barger, A. J., & Wang, W.-H. 2016, ApJ, 829, 25

Hunter, J. D. 2007, CSE, 9, 90

Kelsall, T., Weiland, J. L., Franz, B. A., et al. 1998, ApJ, 508, 44

Khaire, V., & Srianand, R. 2015, ApJ, 805, 33

Khaire, V., & Srianand, R. 2019, MNRAS, 484, 4174

Lagache, G., Abergel, A., Boulanger, F., Désert, F. X., & Puget, J. L. 1999, A&A, 344, 322

Lagache, G., Dole, H., Puget, J. L., et al. 2004, ApJS, 154, 112

Lagache, G., Haffner, L. M., Reynolds, R. J., & Tufte, S. L. 2000, A&A, 354, 247

Mather, J. C., Fixsen, D. J., Shafer, R. A., Mosier, C., & Wilkinson, D. T. 1999, ApJ, 512, 511

Morbidelli, A. 2005, arXiv:astro-ph/0512256

Odegard, N., Weiland, J. L., Fixsen, D. J., et al. 2019, ApJ, 877, 40

Oort, J. H. 1950, BANS, 11, 91

Oteo, I., Zwaan, M. A., Ivison, R. J., Smail, I., & Biggs, A. D. 2016, ApJ, 822, 36

Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016a, A&A, 594, A8

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A13

Planck Collaboration, Aghanim, N., Ashdown, M., et al. 2016b, A&A, 596, A109

Protheroe, R. J., & Meyer, H. 2000, PhLB, 493, 1

Puget, J. L., Abergel, A., Bernard, J. P., et al. 1996, A&A, 308, L5

Renault, C., Barrau, A., Lagache, G., & Puget, J. L. 2001, A&A, 371, 771

Sano, K., Matsuura, S., Yomo, K., & Takahashi, A. 2020, ApJ, 901, 112

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Wright, E. L. 2004, NewAR, 48, 465