The Quest for the Missing Dust. I. Restoring Large-scale Emission in Herschel Maps of Local Group Galaxies

Studies of dust emission in Local Group galaxies often suffer from the specific complexities of far-infrared (FIR) and submm observing. Instruments in this regime generally operate by scanning the sky, with the resulting detector timelines compared and combined in the data reduction process to produce maps. However, the change in flux density recorded by the detector as it scans will not only be due to emission from the target source(s), but can also be caused by drift in instrumental temperature, varying atmospheric foreground emission (for ground-based telescopes), 1/f noise, and other effects. This will introduce considerable artifacts on larger angular scales in the resulting maps (see Roussel 2012; Piazzo et al. 2015, and references therein).

A standard way of minimizing such artifacts is by observing a wide area of background devoid of emission from the target source. Source emission is then constrained relative to this background by comparing overlapping scans that cover both (Meixner et al. 2013). However, observing large amounts of sky is not possible for all observing programs. Plus, background scanning will struggle to save observations from noise that manifests at the same angular scales as the emission from the target object, or where the "empty" background also contains large-scale structure. It is practically impossible to recover emission on scales larger than the size of the map being scanned.

Another technique for removing large-scale instrumental artifacts consists of high-pass filtering the detector timelines (Griffin et al. 2010; Roussel 2012). However, this runs the risk of removing genuine large-scale astronomical emission in the observations (Pascale et al. 2011; Valiante et al. 2016). Additionally, some degree of ringing is also likely to be introduced around areas of compact bright emission (Chapin et al. 2013; Kirk et al. 2018); this ringing will manifest on the same scale as the applied filter.

Ultimately, most reduction treatments for suppressing large-scale artifacts and noise from FIR–submm scan data will result in astrophysical emission on large angular scales also being filtered out of the observations. This is particularly problematic for the Local Group, as it contains the most extended galaxies on the sky—observations of which therefore stand to suffer the most from any removal of FIR–submm emission on larger angular scales. Of particular concern is that the dust giving rise to this filtered-out flux will be the most diffuse dust in these galaxies, found in their outer regions and other low-density environments. This diffuse dust is likely to have distinct properties not found in denser areas, and is important for understanding the evolution of the ISM with environment, especially with regards to G/D (Roman-Duval et al. 2017). The loss of this flux will therefore systematically skew our understanding of the ISM.

There are, however, certain FIR–submm observations that are unaffected by the large-scale noise and filtering issues. All-sky surveys with absolute photometric calibration, such as the Cosmic Background Explorer (COBE; Boggess et al. 1992) and Planck (Planck Collaboration et al. 2011a), can accurately capture emission on all angular scales, down to their resolving limit. However, these observatories have resolutions at least an order of magnitude worse than that achieved by larger telescopes such as the Herschel Space Observatory (Pilbratt et al. 2010); but, of course, larger telescopes are the ones whose data most suffer from suppression of large-scale astrophysical emission. ⁴

This situation—high-resolution data that are missing large-scale flux, and low-resolution data that preserve it—is one that radio astronomers handle frequently. When performing multi-dish interferometry, it is common practice to use low-resolution single-dish data to restore the large-scale flux to which the high-resolution interferometric data are not sensitive, by combining the two data sets in Fourier space; this process is often referred to as "feathering." Despite being a long-established technique in radio interferometry (Bajaja & van Albada 1979), it has only rarely been applied to single-dish FIR–submm observations (e.g., Csengeri et al. 2016; Abreu-Vicente et al. 2017, M. W. L. Smith et al. 2021, submitted).

In this paper, we use the "feathering" Fourier combination approach to produce corrected versions of the Herschel maps of the Local Group galaxies M31, M33, the Large Magellanic Cloud (LMC), and the Small Magellanic Cloud (SMC).

In Section 3, we detail all of the high- and low-resolution input data we employ, including our new Herschel reductions. To feather together two observations requires us to first infer how the emission detected by the low-resolution telescopes would appear if observed in the bandpass used by the high-resolution telescopes; we describe this process in Section 4. In Section 5, we cover the specifics of our Fourier combination process. In Section 6, we describe how we subtract Milky Way foreground emission from the maps we produce; then in Section 7, we present some initial results obtained with our new maps, exploring the properties of the newly restored dust emission in our target galaxies.

The Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on board Herschel performed photometric imaging in 3 FIR bands, at 70, 100, and 160 μm (although the design of the filter wheel meant that only two bands could be observed simultaneously; 160 μm plus one of the others); the typical ⁵ FWHM resolution in these bands is 9'', 10'', and 13'', respectively.

For M31, we used the PACS data from Herschel program GT1_okrause_4. The detector timeline reduction was performed using the Herschel Interactive Processing Environment (HIPE; Ott 2010) v12. ⁶ Map-making was then performed using the Scanamorphos pipeline (Roussel 2012, 2013) v24.0, a pipeline designed to be particularly effective at preserving larger-scale emission and minimizing negative bowl features.

For M33, we retrieved the latest reductions from the Herschel Science Archive (HSA ⁷ ). Specifically, we used the Standard Product Generation (SPG; the standardized automated reductions provided by the HSA) maps, for which the detector timeline reduction was performed by the HSA using HIPE v14, and for which the map-making was carried out using the JScanam pipeline (a modified version of Scanamorphos, now included with HIPE, designed to be robust when run in an automated way; Graciá-Carpio et al. 2017). At 70 and 160 μm, we used the PACS data from Herschel program OT2_mboquien_4; at 100 μm, where data from that program were not available, we instead used the shallower (although wider area) data from program KPOT_ckrame01_1.

3.1.1. PACS Re-reduction for the Magellanic Clouds

For the Magellanic Clouds, we performed our own re-reductions. The Magellanic Clouds were observed by Herschel as part of the Herschel Inventory of the Agents of Galaxy Evolution (HERITAGE; Meixner et al. 2013) key program (KPOT1_mmeixner1_1). The HERITAGE PACS data include both the LMC and SMC at 100 and 160 μm. HERITAGE employed an unusual "basket-weave" observing strategy, where the full width of each Magellanic Cloud was observed with long continuous scans legs in alternating directions; then six months later, a matching set of orthogonal cross-scans was observed. This approach was taken for several reasons: it ensured that all scans sampled the "empty" sky beyond the target galaxies; it maximized the number of cross-scans that each scan overlapped, to allow instrumental variation to be disentangled more effectively; and it increased the fraction of the observing time that was spent integrating along scan legs, as opposed to time spent in turnaround or slewing.

However, this observing strategy had the unintended consequence of leaving the PACS data severely affected by 1/f noise, arising from the extremely long scan legs (8° for the LMC and 6° for the SMC), which are the longest of any Herschel observations. In the PACS detector timelines, instrumental baseline drift entirely dominates over astrophysical signal. While the redundancy from the multiple orthogonal cross-scans means this drift can be well accounted for, doing so requires holding all of the scans in memory, which was not plausible at the time of the original HERITAGE data release. Instead, they adopted a number of other strategies to remove the baseline drift from the scans; see Section 3 of Meixner et al. (2013) for details. The baseline drift removal strategies employed by Meixner et al. (2013) led to significant artifacts around bright sources (illustrated in the left panel of Figure 2), and relied upon tying the surface brightness at the end of each scan leg to a linear interpolation from IRAS and COBE data, leading to conspicuous discontinuities.

Zoom In Zoom Out Reset image size

We therefore needed to perform our own new reductions of the HERITAGE observations of the Magellanic Clouds. We obtained the detector timelines for each individual observation from the HSA, where they had already been reduced using the SPG pipeline in HIPE v14, with the most recent calibration products. To create combined maps from these processed timelines, we used the UNIMAP pipeline (Piazzo et al. 2015), v7.1.0. UNIMAP is a generalized least-squares map-maker, specifically designed to handle data for which the detector timelines suffer from 1/f noise. While most other PACS map-makers (such as JScanam) require specific matched pairs of scan and cross-scan observations in order to function, UNIMAP can operate with arbitrary sets of overlapping scans—a necessary feature given the HERITAGE observing strategy. Recent releases of HIPE provide the option to call UNIMAP for the map-making stage of PACS data reduction, and this was how we used UNIMAP to produce maps for the LMC and SMC, with all processing options set to the default settings. This makes our LMC and SMC reductions substantially similar to those produced by the Herschel Science Centre as part of its PACS "Highly Processed Data Products" data release (Calzoletti 2017), which also used UNIMAP. However, we use a slightly newer version of HIPE than Calzoletti (2017); also, their reductions for the SMC break the data into separate maps for the bar and bridge of the SMC, whereas we produce one contiguous map. A comparison of the old HERITAGE reduction and our new reduction for a portion of the SMC is shown in Figure 2. As an example, our PACS 100 μm reduction for the LMC is shown in the right panel of Figure 3.

Zoom In Zoom Out Reset image size

3.1.2. Herschel-PACS PSF

The Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) on board Herschel provided simultaneous photometric imaging in three submm bands, centered at 250, 350, and 500 μm; these bands achieved resolutions of 18'', 25'', and 36'' FWHM, respectively.

For M31, we used the SPIRE data from Herschel program GT1_jfritz_1, as re-reduced by M. W. L Smith (2021, private communication). The detector timeline reduction was performed using the Bright Galaxy Adaptive Element (BriGAdE; Smith et al. 2012; Smith 2013) pipeline, run with HIPE v12 (see footnote 6). In the case of M31, these data have fewer thermal drift artifacts than the SPG reductions from the HSA.

For M33, we took the latest SPG reductions from the HSA, which use the up-to-date HIPE v14 photometric calibrations. These maps use the data from Herschel program KPOT_ckrame01_1.

In all instances, the data were reduced with the relative gains of the SPIRE bolometers optimized for detecting extended emission, using the values for the beam area values provided in HIPE v15; specifically, 469.7, 831.7, and 1793.5 arcsec² at 250, 350, and 500 μm, respectively. Additionally, all maps were produced using the SPIRE de-striper to mitigate large-scale artifacts arising from instrumental drift.

3.2.1. SPIRE Re-reduction for the Magellanic Clouds

3.2.2. Herschel-SPIRE PSF

The Planck satellite (Planck Collaboration et al. 2011a) surveyed the entire sky in nine submm–microwave bands from 350 μm to 10 mm. Its mission to map the Cosmic Microwave Background required its instruments to have highly accurate absolute calibration, with sensitivity to emission on angular scales spanning the whole range from its primary beam up to the all-sky dipole (Planck Collaboration et al. 2016a, 2016b). This makes it perfectly suited to provide the absolutely calibrated low-resolution data we need at wavelengths ≥350 μm.

We retrieved the 2018 Planck data release all-sky intensity maps from the Planck Legacy Archive. ⁹ Specifically we obtained the 350 and 550 μm maps, both from Planck's high-frequency instrument (Planck Collaboration et al. 2018). These maps have had the zodiacal light contribution removed, are already provided in our desired surface brightness units of MJy sr⁻¹, and are calibrated assuming a constant ν S_ν reference spectrum. The telescope achieves a FWHM resolution of 46 and 48 at 350 and 550 μm, respectively.

For each of our galaxies, we produced a square cutout map centered on the target, with a width equal to 20 times its R₂₅ (the isophotal semimajor axis at which the optical surface brightness falls beneath 25 $\mathrm{mag}\,{\mathrm{arcsec}}^{2};$ see Table 1). We made these cutouts by reprojecting from the HEALPix projection (Górski et al. 2005) of the all-sky maps to a standard east–north gnomic tan projection, using the Python package reproject.

3.3.1. Planck PSF

The Infrared Astronomical Satellite (IRAS; Neugebauer et al. 1984) mapped 96% of the sky in four IR bands, from 12 to 100 μm. The IRAS observations at 100 μm provide the necessary shorter-wavelength counterpart to the Planck data, with a well-matched angular resolution of 5'. Together, IRAS and Planck cover the entire wavelength range of the Herschel data we are reprocessing.

There exist two sets of all-sky IRAS maps; the IRAS Sky Survey Atlas (ISSA; Wheelock et al. 1994), and the Improved Reprocessing of the IRAS Survey (IRIS; Miville-Deschênes et al. 2005). For our purposes, both sets of maps have pros and cons. The newer IRIS maps have been de-striped, de-glitched, zodiacal light subtracted, and had their absolute flux calibration pegged to that of COBE. However, C. Bot et al. (2021, in preparation) have found that the IRIS maps nonetheless show considerable (up to ∼20%) deviation from the COBE data over large angular scales (>a few degrees) in certain regions, such as the Magellanic Clouds—precisely where we need to be able to depend on large-scale fidelity in the IRAS data. While the ISSA maps do not suffer from this specific systematic, they do however lack the improved absolute calibration of the IRIS data, and have a highly nonlinear response function that varies with both the angular scale and surface brightness of the emission being observed. Indeed, the ISSA explanatory supplement (Wheelock et al. 1994) suggests that users should assume photometric uncertainties of up to 30% and 60% at 60 μm and 100 μm, respectively.

For this reason, we opted to use IRIS for our source of IRAS maps. But this required us to first find a way to correct the extended surface brightness deviations in the IRIS data around the Magellanic Clouds. Fortunately, this problem can be corrected in exactly the same manner as the poor preservation of extended emission in the Herschel data.

First, we corrected the IRIS data by feathering it with data from an even-lower-resolution telescope, one that has accurate absolute calibration and preservation of large-scale emission. Then, we used this corrected IRIS data (in concert with the Planck data) to provide the low-resolution counterpart to correct the Herschel and Spitzer data in a second round of feathering, to produce our final high-resolution maps. The obvious choice for the even-lower-resolution data with which to correct the IRIS maps is COBE, which we describe below, in Section 3.5.

We retrieved the IRAS-IRIS data we used from the Infrared Science Archive (IRSA ¹⁰ ). For this work, we only required the IRAS 100 μm maps. For each our galaxies, we produced a square cutout map centered on the target, with a width equal to 20 times its R₂₅. We made these cutouts by mosaicking together all of the individual IRIS maps that covered the region of interest—by first projecting these individual maps to a shared pixel grid (using the Python package reproject), then taking the mean value in any pixel where multiple maps overlapped. This is the same process used in the mosaicking scripts released by the IRIS team. ¹¹ As an example, the IRIS 100 μm map of LMC is shown in the center panel of Figure 3.

3.4.1. IRAS-IRIS PSF

The Diffuse Infrared Background Experiment (DIRBE; Silverberg et al. 1993) instrument on COBE was a single-pixel absolute photometer that observed the entire sky in 10 bands from 1.25 to 240 μm, at an instrumentally limited resolution of 07. DIRBE required accurate absolute calibration for COBE's mission to map both the Cosmic Infrared Background and Cosmic Microwave Background. As a result, its accuracy has been well characterized (Fixsen et al. 1997), and DIRBE serves as the primary photometric reference standard for extended emission across most of its wavelength range (Miville-Deschênes et al. 2005; Gordon et al. 2006).

We retrieved the all-sky COBE-DIRBE maps at 100, 140, and 240 μm from the Legacy Archive for Microwave Background Data Analysis (LAMBDA ¹² ). Specifically, we used the Zodiacal-Light-Subtracted Mission Average Maps. These all-sky maps are presented in the pre-HEALPix quadrilateralized spherical cube projection (Torres et al. 1989). We used these all-sky data to create square cutouts, of width 30°, centered on each of our sample galaxies. To reproject the maps from the quadrilateralized spherical cube projection to the standard gnomonic TAN projection for our cutouts, we used the algorithm supplied by the COBE team (see footnote 12), albeit adapted from the original FORTRAN code into Python (O. Lomax 2021, private communication). As an example, the DIRBE 100 μm map of LMC is shown in the left panel of Figure 3.

We note that there is an alternate set of DIRBE maps hosted at IRSA, stored in a HEALPix projection. However, the surface brightness in these maps, compared to those hosted at LAMBDA, do not agree. When re-gridded to the same projection to allow direct comparison, they can disagree by up to ±15%, particularly surrounding bright sources. It seems that this is due to a slight difference in resolution between the two maps. In communication with personnel at IRSA, it appears that these maps were produced and delivered by the Planck team as part of their data releases (potentially just for quick-look purposes). As the maps hosted at LAMBDA are the original and official maps delivered by the COBE team, we opted to use those. People seeking DIRBE data online should be aware that the different archives provide differing versions of the data.

3.5.1. COBE-DIRBE PSF

We perform our SED fitting using the Dust Brute Force Fitter (DustBFF; Gordon et al. 2014). DustBFF is a grid-based SED fitting code. This is the most efficient solution for fitting the SED of large numbers of pixels in a set of maps.

The BEMBB model implementation in DustBFF does not explicitly parameterize β₂, instead using the 500 μm excess, e₅₀₀, which describes the relative excess in flux at 500 μm, above what would be expected from an MBB model with no break (where negative values indicate a 500 μm deficit), given by:

$$\begin{eqnarray}&&{e}_{500}={\left(\displaystyle \frac{{\lambda }_{\mathrm{break}}}{500\,\mu {\rm{m}}}\right)}^{{\beta }_{2}-{\beta }_{1}}-1.\end{eqnarray} \tag{ 5 }$$

With our BEMBB model, the free parameters are therefore Σ_d , T_d , β₁, λ_break, and e₅₀₀. The parameter grid we use is described in Table 2. When fitting DIRBE data alone (see Section 4.2), we set e₅₀₀ = 0 and fix λ_break, effectively adopting an MBB model, as DIRBE data do not provide sufficient long wavelength coverage to constrain any SED break.

The grid spacing for Σ_d and T_c is base-10 logarithmic; this is to capture the wide dynamic range of possible dust mass surface densities, and the fact that the FIR–submm bands we are modeling are more sensitive to smaller changes in cooler dust temperatures than hotter ones. Our parameter grid effectively imposes a flat prior across the range of values modeled for each parameter (or flat in logarithmic space for Σ_d and T_c ), with a likelihood of zero outside these ranges. The full grid contains over 24 million models. DustBFF computes the probability of fitting the data for each model in the grid. The mathematical formalism under which DustBFF operates is provided in full in Section 4 of Gordon et al. (2014).

The DustBFF formalism incorporates a full covariance matrix, ${ \mathcal C }$, given by:

$$\begin{eqnarray}&&{ \mathcal C }={{ \mathcal C }}_{\mathrm{calib}}+{{ \mathcal C }}_{\mathrm{instr}}\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal C }}_{\mathrm{calib}}$ is a matrix incorporating the covariances between bands due to calibration uncertainty, and ${{ \mathcal C }}_{\mathrm{instr}}$ is a matrix used to capture how the instrumental noise will effect the uncertainty for each model.

It is vital to take account of the correlated uncertainties, as they can lead to systematic errors in model results (Veneziani et al. 2013). For instance, if correlated error in the calibration of the long-wavelength points of a dust SED leads to those fluxes being recorded as greater than they are, this would be erroneously interpreted as favoring a flatter emissivity slope β; accounting for the correlated uncertainties via ${{ \mathcal C }}_{\mathrm{calib}}$ protects against this. We provide ${{ \mathcal U }}_{\mathrm{uncorr}}$ and ${{ \mathcal U }}_{\mathrm{corr}}$ for each instrument in their respective subsections below. This effect can be accounted for by explicitly parameterizing the correlated uncertainties (Galliano et al. 2011; Clark et al. 2019), or by treated them as hyperparameters in a hierarchical Bayesian model (Kelly et al. 2012), or by including them as off-diagonal terms in the covariance matrix (Smith et al. 2012; Gordon et al. 2014), as we do here.

We also incorporate the effect of the instrumental noise in the observations, via the matrix ${{ \mathcal C }}_{\mathrm{instr}}$ (which was not included in the original Gordon et al. 2014 presentation of DustBFF). This allows our fitting to take account of the fact that models that predict a higher surface brightness will be less effected by uncertainty, simply because the noise will be lower relative to the emission. All of the observations we will be fitting using DustBFF have effectively flat noise over the map areas being modeled. We therefore calculate the diagonal values of ${{ \mathcal C }}_{\mathrm{instr}}$ by taking the median of each observation's uncertainty map; the off-diagonal elements are zero. We provide ${{ \mathcal C }}_{\mathrm{instr}}$ for each instrument in their respective subsections below, with a different ${{ \mathcal C }}_{\mathrm{instr}}$ calculated for each galaxy's observations.

In previous uses of DustBFF, ${ \mathcal C }$ also incorporated the correlated uncertainty between bands arising from the uncertainty on the background subtraction (Gordon et al. 2014; Chastenet et al. 2017; Roman-Duval et al. 2017). However, as we are concerned with the entirety of the flux recorded in each pixel, we are not performing background subtraction, so background subtraction covariance is not present.

The full posterior probability distributions returned by DustBFF are extremely large, containing the probabilities of all 24 million models, calculated for every single pixel (of which there are over one million, in the case of our Planck and IRIS maps of the LMC). This is an impractical amount of data to handle, especially given that the vast majority of the models will have essentially nil probability for any given pixel. Therefore, for each pixel, we draw 500 random models from the posterior (with replacement), where the probably of any given model being drawn is proportional to its posterior probability. These 500 realizations of the posterior are what we use to constrain parameter medians, uncertainties, covariances, etc., in each pixel, and are propagated throughout our analysis. For our fiducial model in each pixel, we use the maximum a-posteriori likelihood (MAP) model; our fiducial output maps are comprised of the MAP model values for each pixel.

As described in Section 3, we required DIRBE data to correct large-scale artifacts in the IRIS maps, and feathering together these two data sets first requires fitting the DIRBE SED, pixel by pixel, to infer how the emission observed by DIRBE would appear in the IRIS 100 μm band. To do this, we modeled the pixel-by-pixel SEDs of the DIRBE 100, 140, and 240 μm (i.e., DIRBE bands 8, 9, and 10) observations of our target galaxies, using DustBFF.

All of the data are already at the same 07 instrumentally limited DIRBE resolution, so no convolution was required prior to fitting.

We fit the SED of each pixel using DustBFF. We kept e₅₀₀ = 0 and fixed λ_break at an arbitrary value, because DIRBE does not provide sufficiently long-wavelength coverage to constrain these parameters. Moreover, because we only performed this fitting to predict the IRIS flux at 100 μm, constraining the model at the longest wavelengths was less vital (unlike for the IRIS–Planck SED fitting in Section 4.3). This resulted in a grid of 79,200 models.

The uncertainty on these bands' calibration is dominated by the 10% uncertainty on the absolute gain, which is shared between the bands; the remaining 1%–3% uncertainty for each band represents the photometric repeatability error arising from fluctuations in the relative gains, and is uncorrelated between bands. Therefore, each element of the DIRBE ${{ \mathcal U }}_{\mathrm{corr}}$ matrix has a value of 0.1; while the diagonal elements of the ${{ \mathcal U }}_{\mathrm{uncorr}}$ matrix have values of 0.03, 0.01, and 0.02 (with matrix rows/columns representing the 100, 140, and 240 μm values, respectively).

We calculated the target galaxies' covariance matrices for instrumental noise, ${{ \mathcal C }}_{\mathrm{instr}}$, by taking the average absolute of each band's uncertainty map as the diagonal elements. The diagonal elements for each galaxy and band are given in Table 3; the off-diagonals are all zero. As can be seen, each band has similar noise across all maps.

We validated the SED fits by checking the residual between the model and the observed surface brightness. In Figure 4, we plot the relative residuals, defined as ${{\rm{\Delta }}}_{\mathrm{rel}}=({S}_{\mathrm{mod}}-{S}_{\mathrm{obs}})/{S}_{\mathrm{obs}}$, for the MAP model of every pixel. For all galaxies and bands, the median Δ_rel was <0.7%, with >61% of pixels having ∣Δ_rel∣ < 5%. This indicates that DustBFF suffered no systematic problems fitting the data.

Zoom In Zoom Out Reset image size

With the DustBFF SED fitting completed, we calculated how each model would appear if it were observed in the IRIS 100 μm band. This was done by convolving every model SED, in each pixel, through the IRAS 100 μm filter, color-corrected to match the IRIS reference spectrum. The result was 500 images, representing the IRIS 100 μm fluxes predicted from the realizations of the posterior of the SED fit for every pixel, along with the fiducial image, containing the MAP model prediction for every pixel—all at DIRBE resolution. An illustration of this for a particular pixel can be seen in Figure 5.

Zoom In Zoom Out Reset image size

Our predicted IRIS 100 μm band maps at the DIRBE resolution scale, with uncertainties constrained by our posterior sampling, are now suitable to be feathered with the the unfeathered IRIS maps. This process is described in Section 5.1.

We fit the SED of IRIS and Planck data in order to predict how the emission observed by these facilities would appear in the Herschel bands; these predicted maps can then be used to feather together the high- and low-resolution data. Note that in this section we use the feathered IRIS data that have already been combined with the DIRBE data, as described later in Section 5.1. ¹³

Each of the IRIS and Planck bands are at slightly different resolutions, so prior to fitting we convolved each to the 48 poorest resolution of the Planck 550 μm band. We did this using conversion kernels creating with the Python package photutils, using the instrumental PSFs for each band.

We modeled the pixel-by-pixel SEDs of the IRIS 100 μm and Planck 350 and 500 μm observations of our target galaxies using DustBFF. We used the full grid of 24,235,200 models. The correlated uncertainty matrix used for the IRIS–Planck SED fitting was:

$$\begin{eqnarray}{{ \mathcal U }}_{\mathrm{corr}}=\left[\begin{array}{ccc}0 & 0 & 0\\ 0 & 0.05 & 0.05\\ 0 & 0.05 & 0.05\end{array}\right]\end{eqnarray} \tag{ 7 }$$

where the rows/columns represent the 100, 350, and 550 μm values, respectively. The uncorrelated uncertainty matrix, ${{ \mathcal U }}_{\mathrm{uncorr}}$, had diagonal elements of 0.13, 0.014, and 0.011 (with off-diagonals of zero).

The 5% correlated error between the Planck bands comes from uncertainty on the emission models of Uranus and Neptune, the calibrator sources for Planck's high-frequency instrument (Bertincourt et al. 2016). As the IRIS 100 μm band is calibrated independently from the Planck bands, it has no correlated uncertainty here. The 13% uncorrelated uncertainty on the IRIS 100 μm calibration is taken from Miville-Deschênes et al. (2005). This combines the 10% uncertainty on the absolute calibration of DIRBE, with which IRIS is calibrated, and an additional 3.7% from uncertainty in the color corrections used by Miville-Deschênes et al. (2005) to relate IRIS to DIRBE. The 1.4% and 1.1% uncorrelated uncertainties on the Planck 350 and 500 μm bands are the statistical repeatability noise (Bertincourt et al. 2016).

The diagonal values of the instrumental noise covariance matrices, ${{ \mathcal C }}_{\mathrm{instr}}$, are given in Table 4, having been calculated from the average of the absolute values in the uncertainty map for each galaxy and band. For the 100 μm IRIS data, the uncertainty maps are part of the output from the DIRBE–IRIS feathering (detailed in Section 5.1).

As with the DIRBE SED fitting, we checked the SED fits by inspecting the relative residuals between the MAP model and observed surface brightnesses, for every pixel. This is plotted in Figure 6. There are definite systematic offsets apparent in the relative residuals. Fortunately, however, they remain very small in absolute terms. For IRIS 100 μm, the band with the worst offset, the median Δ_rel was only −0.91%, with a standard deviation of 0.99%, with >95% of pixels having ∣Δ_rel∣ < 2.1%. For the Planck bands, median Δ_rel values are −0.40% and −0.29% at 350 and 500 μm, respectively, with >95% having ∣Δ_rel∣ < 0.91% in both cases. Given that the systematics are so small, being at the subpercent level in all cases, with symmetrical distributions, we believe the SED fits to be satisfactory. It is odd that the offsets are negative in all bands, but we suspect this may arise from the logarithmic spacing of the dust mass surface density (and/or temperature) increments in the parameter grid.

Zoom In Zoom Out Reset image size

Having completed the IRIS–Planck SED fitting, ¹⁴ we used the resulting posterior to predict how the modeled emission would appear in the Herschel 100–500 μm bands, convolving each model through the Herschel filters, color-corrected to match the Herschel reference spectrum. An example of this for the model SEDs of a particular pixel are shown in Figure 5. The output was 500 images, representing the Herschel fluxes in each band predicted from the realizations of the posterior of the SED fit for every IRIS–Planck pixel, along with that of the fiducial maximum a-posteriori (MAP) image, containing the MAP prediction for every pixel, all at DIRBE resolution.

We originally attempted this SED fitting process using the Planck 850 μm band as well, to provide additional constraints for the models. However, we found that this resulted in considerably larger residuals, positive and negative, dominated by the structure of Cosmic Microwave Background radiation, which is more prominent at 850 μm. For this reason, we opted to not use the Planck 850 μm data.

The only band for which we feathered together DIRBE and IRIS data was at 100 μm. The DIRBE data we used are the maps we produced in Section 4.2, predicting the emission that DIRBE would have observed in the IRIS 100 μm band. For our fiducial feathered maps, we combined the IRIS cutout for each galaxy with the image produced from the MAP outputs for the pixel-by-pixel DIRBE SED fitting (Section 4.2).

For each target galaxy, we reprojected these DIRBE data to the IRIS pixel grid. We also reprojected the DIRBE and IRIS PSF maps to have the same dimensions and pixel widths as the IRIS observation. As part of the reprojection, we also had to apodize the DIRBE maps, to remove pixel-edge artifacts from the re-gridded maps; we detail this in Appendix D.

With all of the data on the same pixel grid, we transformed the observations and beams to Fourier space. Next, we deconvolved the DIRBE data with the DIRBE beam, then convolved it with the IRIS beam. This accounts for the brightness differences between the data sets that were caused only by the differences in their resolution; i.e., where a given source would have had a lower peak brightness in the low-resolution data. The DIRBE data now have the same beam amplitude and effective resolution as the IRIS data, ¹⁵ meaning the two can be directly compared and combined.

The spatial scales at which systematic response errors are observed in the IRIS data (∼1°; C. Bot et al. 2021, in preparation) are only somewhat larger than the spatial scales above which DIRBE is sensitive, as dictated by its 07 beam. The power spectra of both data sets, for the example of the SMC, are shown in Figure 7. The amplitude of the IRIS power spectrum is not considerably different from that of DIRBE at the largest scales, because the errors in IRIS are localized only around certain parts of the LMC and SMC; the IRIS data are otherwise well behaved. This can also be seen in the fifth column of Figure 9.

Zoom In Zoom Out Reset image size

For feathering together the DIRBE and IRIS data in Fourier space, we chose to follow the weighting strategy used in the Common Astronomy Software Applications (CASA; McMullin et al. 2007) task feather. ¹⁶ In this approach, the Fourier transform of the low-resolution beam is used as the weighting function. Therefore, at the largest scales (i.e., the zeroth Fourier mode), the DIRBE data were assigned a weight of 1. Going to smaller spatial scales, the weighting of the DIRBE data decreased, while the weighting of the IRIS data increased; the sum of both weights remained 1 everywhere. At scales smaller than the DIRBE beam FWHM, the weights crossed over so that the IRIS data dominated and became increasingly heavily weighted. This is illustrated in Figure 8.

Zoom In Zoom Out Reset image size

The CASA algorithm minimizes ringing by making the weighting transition as smooth as possible, tapering according to the low-resolution beam, using its FWHM as the natural crossover scale. A disadvantage of this approach is that the high-resolution data will nonetheless be assigned some (albeit minimal) weighting at the larger scales, where they have poor sensitivity (and vice versa for the low-resolution data). However, if one were to instead attempt to "squeeze" the transition region into a smaller range of scales, where the two data sets' sensitivity overlaps, there will be considerable ringing in the combined image if the transition window is too narrow. This is what would happen if we attempted to combine the DIRBE and IRIS data with such a transition window: the small-scale edge of such a window would be dictated by the 07 DIRBE beam, while the large-scale edge of the window would be dictated by the ∼1° scale of the IRIS response errors around the LMC and SMC.

We therefore opted instead to take advantage of the ringing minimization provided by the CASA algorithm. We tested this beam tapering method by feathering simulated high- and low-resolution data, in Appendix C.

Comparison of the DIRBE and IRIS power spectra (Figure 7) suggests there are no significant deviations between the brightness scales. IRIS was constructed to be cross calibrated to the DIRBE brightness scales, so this is not surprising (indeed, it seems that the IRIS response errors around the LMC and SMC are localized failures of this cross-calibration). Because the largest scales of the feathered maps are dictated entirely by the DIRBE data, the zero level of the output maps will match that of DIRBE, by construction.

We carried out the feathering process 500 times, once for each sample of the DIRBE SED fitting posterior interpolations produced in Section 4.2. We also feathered the fiducial DIRBE data with the unfeathered IRIS data to produce our fiducial IRIS+DIRBE map. For each of these 500 iterations, we permutated the surface brightness values in each IRIS pixel according to the IRIS uncertainty maps, assuming a random Gaussian noise distribution. The 500 resulting feathered maps therefore provide a bootstrapped Monte Carlo estimate of the uncertainty on the feathered output in each pixel.

5.1.1. Validation of Feathering IRIS with DIRBE

The results of the IRIS+DIRBE feathering are shown in Figure 9. The first three columns compare the IRIS, DIRBE, and feathered IRIS+DIRBE data for each galaxy. The last two columns show our two main validation tests for the feathering process.

Zoom In Zoom Out Reset image size

The fourth column of Figure 9 shows the relative residuals between the DIRBE map and the feathered map; for this comparison, we convolved the feathered map to the DIRBE resolution. Ideally this residual should be 0% everywhere, as the feathered maps should conform to the DIRBE data at scales ≥ the DIRBE beam. There are slight ringing artifacts visible around the LMC and SMC, but at worst they are ±3% (for the LMC). This is much lower than either instruments' calibration uncertainty.

The very small residuals between the feathered and IRIS data also compare very favorably to the fifth column of Figure 9. This column shows the relative residuals between the feathered data and the unfeathered IRIS data; this illustrates the magnitude of the correction applied to the IRIS data by the feathering process. As before, this comparison is plotted at the resolution of DIRBE. The feathered data have significantly more surface brightness over the LMC and SMC versus in the unfeathered IRIS data, correcting the response errors. The change in surface brightness is up to 10% for the LMC and up to 5% for the SMC. M31 and M33 are not bright and/or extended enough to have suffered from the IRIS response issue, so do not display these changes in the feathered maps. The data for all four galaxies also show a reduction in surface brightness over the background surrounding all galaxies, of ∼8% in all cases; this is due to an apparent difference in the zero level of the DIRBE and IRIS data.

The peak of 10% extra flux in bright regions, superimposed over an ∼8% reduction in background level, corresponds to an up to 18% correction to the background-subtracted surface brightness for the Magellanic Clouds. This correction is substantial, and is clearly necessary to obtain good science analysis with the data; it is also in line with the ≈20% IRIS response errors measured in other work (C. Bot et al. 2021, in preparation).

Here, we feather together our Herschel reductions with our inferred maps of large-scale emission in Herschel bands (produced using IRIS and Planck data in Section 4.3). We carried out this process for the 100, 160, 250, 350, and 500 μm data.

We first generated power spectra for the IRIS–Planck and Herschel data for each of our galaxies and bands, some examples of which are shown in Figure 10; these plots also show the ratio between the amplitudes of the two power spectra. At the largest angular scales, the Herschel power spectrum tends to fall beneath the IRIS–Planck power spectrum, by as much as 50% (e.g., see the SMC 500 μm panel of Figure 10).

Zoom In Zoom Out Reset image size

We note that in some instances no power on large scales appears to be missing, such as at 500 μm for M31 or M33; this is potentially due to the smaller angular size of these galaxies, care taken during the reduction process, and the reduced filtering needed at 500 μm thanks to its larger beam. However, it should be remembered that even if the power spectra for Herschel show no deficit in power compared to IRIS–Planck, that does not preclude the possibility of the emission being distributed differently. For instance, as described in Section 5.2.1, the LMC, SMC, and M31 data all show gradients in the foreground cirrus emission in the IRIS–Planck data that are considerably attenuated in the Herschel data; the presence (or not) of these gradients significantly effect the emission measured in the outskirts of our target galaxies after background subtraction (see Section 6). We therefore apply the full feathering process to all bands, even when no difference appears in the power spectra, to deal with these sorts of problems, and in order to treat all of our data consistently throughout.

Inspecting the power spectra in Figure 10 indicates that there is a reasonably wide range of spatial scales over which both sets of data are sensitive to emission. We therefore used a feathering technique that could take advantage of this. Within this window of overlapping spatial scales, a smooth tapering function was used to mediate the transition in weights, from using the low- to high-resolution data. At spatial scales larger than the taper window, only the low-resolution data were used; and at scales smaller, only the high-resolution data. The overlap in spatial sensitivity between both sets of data also made it possible to cross-calibrate them, to ensure that the low- and high-resolution data were both on the same brightness scale before being combined.

To find an appropriate tapering window, we inspected the power spectra to identify a specific spatial scale interval for which the ratio between the high- and low-resolution amplitudes remained relatively constant. This constant ratio indicates that the two sets of data are maintaining a constant spatial sensitivity, except for a difference in flux calibration (if the ratio is 1, then the two data sets already have the same flux calibration). It was necessary to balance the benefits of a larger tapering window (less ringing in the feathered outputs; more values with which to perform brightness scale cross-calibration), versus the benefits of a narrower tapering window (less risk of extending sampling to spatial scales where either of the data sets starts to lose sensitivity).

Even within the tapering window, the ratio between the high- and low-resolution amplitudes does not remain perfectly constant. To check if there were systematic variations in the ratio within each tapering window, we fit straight lines to the amplitudes inside them; ideally the gradient of such a line should be zero. We performed a Monte Carlo bootstrap resampling to evaluate the gradient uncertainty. For 17 of the 20 power spectra (i.e., for each band for each galaxy), the gradient of the ratio within the tapering window was compatible with being 0, to within the uncertainty. For the remaining three, the deviations were 1.03σ, 1.20σ, and 2.15σ. Given that 1 out of 20 should deviate by >2σ on average for Gaussian statistics, we are satisfied that this indicates the amplitude ratios within the tapering windows are well behaved.

After iterating through a range of options and inspecting the outputs, particularly the residual plots (see Section 5.2.1 and Figure 12), we found that a tapering window of 30'–60' worked best; this range is marked in Figure 10. Specifically, within this window of overlapping sensitivity, we used a half cosine bell function (i.e., half of a Hann window function; Harris 1978) to smoothly mediate the transition in weights applied to the data sets. This transition is illustrated in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based upon a wide suite of feathering simulations, Kurono et al. (2009) find that the high- and low-resolution observations should have a spatial scale overlap spanning at least a factor of 1.7, in order to achieve outputs with minimal errors (above this overlap factor, they find minimal, asymptotic improvements). Given that our 30'–60' tapering window spans a factor of 2, we comfortably satisfy this criterion.

When feathering the PACS data for the LMC, we found that shifting the window function to larger angular scales was necessary, due to an apparent nonlinearity in the response function of our IRIS–Planck data around the extremely high surface brightness star-forming complex of 30 Doradus. This response nonlinearity caused the IRIS–Planck data to overestimate the surface brightness at the center of 30 Dor, and underestimate it in the surrounding area (similar response issues around 30 Dor were discussed in Meixner et al. 2013). These errors appear at scales of up to ∼40'. As we describe in Appendix E, it appears that the PACS data correctly recover the emission around 30 Dor at these scales. By putting our tapering window at a larger scale, we ensure the emission from 30 Dor is correctly reproduced in the final feathered maps. We experimented with a range of tapering windows, and found that 45'–90' stops emission from being lost around 30 Dor, while still allowing the IRIS–Planck data to correct emission at larger scales. This maintains a tapering window spanning a factor of 2 in angular scale, so will satisfy the Kurono et al. (2009) factor >1.7 overlap criterion just as well as for the other galaxies.

To cross-calibrate the two data sets to the same brightness scale (i.e., ensure the gain for both agree), we divided the low-resolution amplitudes by the high-resolution amplitudes within the overlap window, then took the median of these ratios. We multiplied the high-resolution amplitudes (over all scales) by this ratio to fix them to the brightness scale of the low-resolution data (which have absolute calibration ultimately pegged by DIRBE and Planck). We estimated uncertainties on the cross-calibration factors by performing 100 bootstrap resamplings (with replacement) of all of the amplitudes in the overlap window, and re-computing the correction factor each time; the standard deviation of these values was taken as the uncertainty.

The cross-calibration factors for each band and galaxy are listed in Table 5. At 500 μm, relatively little correction is necessary, with all factors being >0.928. However, going to shorter wavelengths, especially into the PACS bands, larger corrections are necessary; the average correction at 100 μm is 0.706. The fact that PACS measures emission as being 20%–30% brighter than IRIS (and Planck) is all the more striking given that we have shown that PACS can miss a considerable amount of emission at large scales. This means that PACS must be overestimating the brightness it does detect by even more than 20%–30%. When comparing PACS and IRIS–Planck at a common resolution without applying the cross-calibration corrections, we found that regions of compact sources were indeed much brighter in PACS, over 40% so in some cases.

PACS is flux calibrated by reference to a set of five standard stars (Balog et al. 2014); the fact this calibration is based off point sources may explain why the calibration appears to struggle at the highly extended scales we are working with. Unlike SPIRE, PACS makes no corrections to the relative gains of the bolometers to account for the different illumination by the beam when making maps of extended emission. This may contribute to the need for a large gain correction. We also note that the two largest correction factors are for the LMC and SMC 100 μm maps; it is therefore possible that the issue is partly due to the UNIMAP pipeline, perhaps in how it handles the unusual scan strategy for the LMC and SMC observations.

It is conceivable that this calibration discrepancy is not present at scales smaller than the 6' low-resolution beams, beneath which we cannot compare the Herschel data to the IRIS–Planck data. It is therefore possible that point sources and compact sources are accurately represented in the unfeathered data. However, this would require a very abrupt shift from well-calibrated to poorly-calibrated as angular scale increased. If this is the case, then our cross-calibrated maps would make point sources too faint, while accurately representing extended emission (compared to the original maps making diffuse emission too bright, but accurately representing point sources).

Previous authors have compared Herschel data to all-sky survey data, such as from IRAS and Planck, to check for differences in gain. For instance, Molinari et al. (2016) compared Herschel and IRIS data for observations of the Galactic plane (which therefore contain highly extended emission), and did not find evidence for gain discrepancies >10%. However, they made this comparison pixel by pixel, not in Fourier space, which could give undue weight to compact features. In contrast, Abreu-Vicente et al. (2017) do make such a comparison in Fourier space, also for Galactic plane fields—including data from Molinari et al. (2016)—and find average cross-calibration factors of 0.84, 0.88, 0.91, and 0.97 necessary at 160, 250, 350, and 500 μm, respectively (measured between 7'–100'). These factors all lie within the range of values we find for each of these bands.

Regardless, the feathering process requires us to cross-calibrate the two data sets being combined to the same brightness scale, otherwise we find that the differences lead to very pronounced artifacts—and over the scales we are able to probe, the Herschel response (especially for PACS) is too bright compared to IRIS–Planck (and therefore DIRBE also). Given that the science in this paper, and the subsequent papers in this series (C. J. R. Clark et al. 2021, in preparation), is concerned primarily with the diffuse ISM, our priority is that this is correctly calibrated—which is what our cross-calibration ensures.

We also note that the cross-calibration factors determined are quite robust against the specific choice of tapering window. As can be seen from inspecting Figure 10 (especially the plots for PACS bands), the offsets between the uncorrected Herschel amplitudes and the IRIS–Planck amplitudes stay remarkably constant over all scales large enough for IRIS–Planck to be sensitive.

Before feathering, we reprojected the low-resolution data to the pixel grid of the Herschel data in each band, apodizing it to remove pixel-edge artifacts (as per Appendix D). For each band, both sets of data, and their corresponding beams, were Fourier transformed. The low-resolution data were deconvolved with their beam, then re-convolved with the band's corresponding Herschel beam. Having done this, we combined the two sets of data in Fourier space, using our tapering function to govern the transition within the tapering window.

As with IRIS+DIRBE, we performed the feathering 500 times for each galaxy and band; once for each sample of the SED fitting posterior in each pixel. For each Monte Carlo iteration, we drew surface brightness values for each Herschel pixel according to the uncertainty map in each band, assuming a Gaussian uncertainty. Recall that the uncertainties on the IRIS–Planck SED fitting factored in the uncertainties on the IRIS+DIRBE feathering, which in turn factored in the uncertainties on the DIRBE SED fitting (and the IRIS photometric uncertainties). So the 500 bootstrap samples produced for each feathered map encompass the propagated uncertainty from all previous stages of our process.

The abrupt edge of the Herschel maps results in considerable edge effects at the borders of the feathered maps. We therefore constructed a mask to describe the portion of each feathered map far enough from the edge to not be affected. Specifically, we found that excluding a border of three times the low-resolution beam (i.e., 145) provided a sufficient buffer. We excluded this region from our various diagnostic plots assessing the feathered data.

5.2.1. Validation of Feathering Herschel with IRIS–Planck

The results of feathering together the Herschel and IRIS–Planck data are shown in Figure 12, for the example case of the 250 μm data for each; corresponding figures for 100, 160, 350, and 500 μm are shown in Appendix F. A zoomed-in view of the SMC wing, contrasting the unfeathered and feathered maps, is shown in Figure 13.

Zoom In Zoom Out Reset image size

The fourth column of these figures show the relative residuals between the IRIS–Planck map and the feathered map (with the feathered map convolved to the low-resolution beam for comparison); Ideally this residual should be 0% everywhere. The fifth columns show the relative residuals between the feathered data and the unfeathered Herschel data, illustrating the degree to which the unfeathered was corrected by feathering; note that the zero level of these plots is somewhat arbitrary. ¹⁷ If the correction in the fifth column is larger than the residual in the fourth column, then the feathering process has done "more good than harm."

In Figures F1 and F2 in Appendix F, we can see that the average residual in the fourth column is zero for every galaxy, as we would wish. There is some ringing apparent surrounding bright features, such as a star-forming complex in the SMC wing, and the rings of M31; these residuals peak in the ±5%–8% range (the background regions where there is no ringing have residuals ≤1%). By contrast, the flux corrections in the fifth column show that the surface brightness was increased by 10%–13% over most of the outskirts of the target galaxies. The distribution of emission restored by feathering is exactly as we would expect—regions of bright compact emission show little-to-no correction, whereas in the diffuse outer regions there has been significant addition. We can also see that feathering has corrected the flat backgrounds that the Herschel reduction tends to impose, restoring the gradients due to foreground cirrus; this is especially clear to the northwest and southeast of the SMC.

These diagnostic plots show some noteworthy differences in the feathering results in different bands. The ringing in the residuals in the fourth column are less pronounced in the other SPIRE bands, at 350 and 500 μm, (Figure F2), falling to ±3.5% at 500 μm. The flux corrections in the fifth column also become smaller toward longer wavelengths. By 500 μm, little-to-no correction is being made to M31 and M33, with only slight changes to the foreground gradient apparent. This is borne out by the power spectra, which show no significant power missing in the Herschel data, compared to the IRIS–Planck data at 500 μm for either galaxy (e.g., see the lower left panel of Figure 10).

The feathering diagnostic plots for the PACS data at 100 and 160 μm (Figure F1) show even more dramatic corrections to the surface brightness, as can be clearly seen in the fifth columns. The restored surface brightness represents a >30% increase over large areas around all of the galaxies, ¹⁸ for both PACS bands, once again in line with expectations from the power spectra in these bands (see Figure 10). Unfortunately, the noise in the residuals between the feathered data and IRIS–Planck data in the fourth columns can be considerable. The ringing around M31 reaches 17% level for 160 μm and 20% for 100 μm. Nonetheless, this is still much smaller than the correction applied by feathering, so the new maps are still more correct than the unfeathered maps, even with this residual noise. The only exception to this is the region of sky to the southeast of M33, where there is a large ≈40% positive residual in the fourth column, while the correction in the fifth column only reaches ≈30%. In general, large residuals in the fourth column are limited to the regions of empty background; within the galaxies (as outlined by the contours), the mean absolute residuals of the fourth column are ≈4.5% at 100 μm, and ≈4.0% at 160 μm. This compares favorably to the 7% photometric calibration uncertainty of the PACS instrument.

As noted in Section 5.1.1 and Appendix E, our IRIS–Planck maps do not correctly recover the emission around 30 Dor in the LMC. So the large residuals around 30 Dor in the fourth column of Figure F1 for PACS are expected (and indeed desirable), indicating where there is a disagreement between the feathered data and the IRIS–Planck data, due to the PACS data correcting IRIS–Planck surface brightness artifacts at smaller scales.

Each of our feathered maps has an associated uncertainty map, which captures the uncertainties propagated through the entire process. In most bands, the increase in the uncertainty maps is negligible (especially in the PACS bands, where the instrumental uncertainty levels were already quite high). However, in the SPIRE bands the levels in the uncertainty maps can increase by a factor of a few, mainly because the instrumental noise levels in the original maps were very low to start with. Even with this increased uncertainty, the cirrus emission in faint sky regions away from our target galaxies is still consistently detected with SNR > 10 in the SPIRE bands.

Because the velocity of the H i emission from the Milky Way and M31 overlap, we could not use the method described above for foreground subtraction. Fortunately, M31 sits behind a region of Galactic cirrus that exhibits a fairly smooth gradient in surface brightness, so the more detailed subtraction method is not vital for success. Instead, we first masked the Herschel observation footprint from the predicted map of emission in each Herschel band derived from the IRIS–Planck SED fitting, the same as above. However, we then performed the polynomial modeling directly on the FIR maps, providing an estimate of the foreground emission in front of M31 in each band, which we then subtracted from the feathered Herschel data. We calculated the uncertainty on this subtraction in the same manner as above.

In order to compare our new, feathered Herschel maps to the unfeathered original maps, we also needed to foreground subtract the original maps. In the case of the LMC, and of the PACS data for M33, there is very little sky surrounding the target galaxies. Without sufficient area of empty sky, it was not possible to perform the same sort of subtraction as we did in Section 6, using Galactic H i data.

Nor could we use the IRIS–Planck FIR data, as the large-scale emission that fills the sky in that data is simply not present in the unfeathered Herschel data. For instance, the sky to the far southeast of the SMC in the Herschel-PACS data has the same surface brightness as the sky to the far northwest, due to the flattening imposed by the Herschel reduction—but in the IRIS–Planck data, the surface brightness of these areas of sky differs by more than a factor of 2, due to variations in diffuse Galactic emission. Subtracting this from the unfeathered Herschel maps would lead to major negative features in large portions of the map. However, the unfeathered Herschel maps do preserve some large-scale foreground emission (especially for M31 and the SMC), so we could not simply treat the background as being flat.

The tried-and-true method of placing an annulus around the target source to measure the sky level is also not ideal here. We are interested in diffuse emission at galaxy outskirts, and any emission in or beyond a sky annulus would therefore be subtracted. And because our targets are large (and sometimes irregular) galaxies, in relatively tight maps, any sky annulus we place would run into this problem. Instead, we masked out the central regions of each unfeathered Herschel map, leaving only a border along the inside edge of the Herschel footprint, ¹⁹ with width equal to 5% the R₂₅ of the target galaxy. We then fit a first-order two-dimensional polynomial (i.e., a tilted plane) to this border data to model the foreground emission, which we then subtracted. This maximized our ability to recover diffuse emission in the outskirts of the target galaxies.

Our new maps change the photometry measured for the galaxies in our sample. The global fluxes we record are given in Table 6, which also lists the fluxes measured from the unfeathered maps, for comparison.

For the LMC and SMC, these fluxes were measured by taking the total flux within the Herschel footprint of the foreground-subtracted feathered maps. We opted not to place a tighter aperture around the sources in order to not exclude any faint emission at large radii. Our foreground selection should be equally good across the entire Herschel footprint, so no systematic bias should be suffered from measuring the fluxes this way. Nonetheless, this will drive up the "aperture noise" associated with the resulting fluxes.

For M31 and M33, the fluxes were measured within elliptical apertures, centered on the coordinates given in Table 1. For M31, the aperture had $a=126^{\prime}$, $b=46^{\prime}$, and θ = 38°; for M33, the aperture had $a=46^{\prime}$, $b=28^{\prime}$, and θ = 20° (where a, b, and θ are the semimajor axis, semiminor axis, and position angle, respectively). These apertures were designed to comfortably contain all of the potential extended emission, especially that visible in the SPIRE maps. ²⁰

For the uncertainty on the recorded fluxes, we used the fractional uncertainty on the foreground subtraction, and the cross-calibration correction, added in quadrature to the instrumental calibration uncertainty (7% for PACS bands, 5.5% for SPIRE).

For the SPIRE bands, the feathered fluxes are in close agreement with the unfeathered fluxes for the LMC and M33. For the SMC, the feathered fluxes are 9%–21% brighter; and for M31, the feathered SPIRE fluxes are 9%–14% brighter (differences decreasing toward longer wavelengths, as would be expected).

In the PACS bands, there is relatively little change in the measured fluxes. Only the 160 μm flux for the LMC and the 100 μm flux for M33 differ by larger than the uncertainty. In these cases, the feathered fluxes are actually less than the unfeathered measurements, by 18% and 43%, respectively. We suspect this change is due to the feathered maps allowing better foreground subtraction; a filament of foreground cirrus passes east–west over M33, and a larger cloud of cirrus extends over the galaxy's southwest quadrant. For the other galaxies, where the change in PACS flux is within the uncertainty, the feathered fluxes are brighter in five out of six of them, for an average change of +7%.

Note, however, that all of these comparisons are occurring after the application of the cross-calibration correction factors during the feathering process, which have therefore not been applied to the unfeathered maps. These correction factors range impose a 14%–31% reduction in flux in the PACS bands, and a 0%–16% reduction in the SPIRE bands (Table 5). Therefore, had these corrections not been applied, the feathered data would have been much brighter, in almost all instances.

In short, while the surface brightness measured at any given point in our target galaxies will often have changed dramatically between the unfeathered and feathered maps (i.e., Figures 13, F1, and F2), the total measured flux tends to be quite similar—thanks to the reduction in flux due to the cross-calibration correction being comparable to the increase in flux from restored emission.

For the LMC, our PACS fluxes agree within the uncertainties with those previously reported from IRAS (Rice et al. 1988), COBE (Israel et al. 2010), and Spitzer (Dale et al. 2009), but differ from those reported by Meixner et al. (2013) with the HERITAGE data. At 100 μm, we find a flux 18% larger than did Meixner et al. (2013), and at 160 μm we find a flux 16% less than that in Meixner et al. (2013). For SPIRE, however, our LMC fluxes agree closely with those of Meixner et al. (2013; which agrees well with COBE at 240–250 μm; Israel et al. 2010).

For the SMC, our PACS fluxes once again agree with IRAS (Rice et al. 1988), COBE (Israel et al. 2010), and Spitzer (Leroy et al. 2007; Dale et al. 2009), and once again differ from those reported by Meixner et al. (2013) with the HERITAGE data; ours being 14% brighter at 100 μm, and 49% brighter at 160 μm. Our SPIRE fluxes for the SMC are also larger than those reported by Meixner et al. (2013), by 29%, 21%, and 12%, at 250, 350, and 500 μm, respectively, whereas our 250 μm flux agrees well with the published 240 μm COBE flux (Israel et al. 2010). We note that the Meixner et al. (2013) SMC SPIRE fluxes are, however, a close match to our unfeathered fluxes.

For the Magellanic Clouds, it appears that matching the calibration of our feathered data to well-calibrated all-sky surveys has fixed the discrepancy between photometry previously measured from Herschel and photometry measured from other facilities (including all-sky surveys).

For M31, our fluxes agree with those reported by Viaene et al. (2014) to within 10% at 350 μm, and to within 5% in all other bands; this is within our uncertainties in all cases. The agreement in the PACS bands indicates that the additional flux recovered in our maps closely matches the relative decrease in flux imposed by our cross-calibration, as suggested in Section 6.3.

For M33, our SPIRE fluxes agree to within 6.3% with those reported by Hermelo et al. (2016); this is within the mutual uncertainties in each case. For the PACS bands, however, our fluxes are 12% and 18% fainter at 100 and 160 μm, respectively. This indicates that relatively less flux was restored by our feathering, as compared to the cross-calibration decrease; this is not surprising, given that M33 is by far the most compact of our sources, and therefore will have suffered the least emission lost by filtering.

We constructed Σ_H maps for each of our galaxies, using 21 cm maps to trace the H i and CO maps to trace the H₂. For the LMC, we used the H i data of Kim et al. (2003) and CO data of Wong et al. (2011). For the SMC, we used the H i data of Stanimirovic et al. (1999) and the CO data of Mizuno et al. (2001). For M31 we used the H i data of Braun et al. (2009), and the CO data of Nieten et al. (2006). For M33, we used the H i data of Koch et al. (2018) and the CO map of Gratier et al. (2010) and Druard et al. (2014). We calculated the molecular gas surface density, Σ_H2, using the standard relation:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{H}}2}={\alpha }_{\mathrm{CO}}{I}_{\mathrm{CO}(1-0)}\end{eqnarray} \tag{ 8 }$$

where I_CO(1−0) is the velocity-integrated main-beam brightness temperature of the CO(1-0) line (in K km s⁻¹), α_CO is the CO-to-H₂ conversion factor (in K⁻¹ km⁻¹ s M_⊙ pc⁻²).

All of the CO observations we used were of the CO(1−0) line, except for that of M33, which was of CO(2−1). For those data, we inferred I_CO(1−0) by applying a line ratio, r_2:1 = I_CO(2−1)/I_CO(1−0). In spiral galaxies, r_2:1 varies radially (Casoli et al. 1991; Leroy et al. 2009). We therefore used the procedure laid out in Section 3.4 of Clark et al. (2019), where r_2:1 is determined as a function of galactocentric radius as a fraction of the R₂₅, using measurements made by Leroy et al. (2009)—varying from 1.0 at the galaxy center to 0.55 at the R₂₅. We thereby computed I_CO(1−0) for each pixel in the M33 I_CO(2−1) map according to its position.

To compute Σ_H2 for M31 and M33, we used the standard Milky Way value of α_CO = 3.2 K⁻¹ km⁻¹ s pc⁻², which tends to be applicable to high-metallicity spiral galaxies in general ²¹ (see Saintonge et al. 2011, Bolatto et al. 2013, and references therein). For the LMC and SMC, we used α_CO values of 6.4 and 21 K⁻¹ km⁻¹ s pc⁻², respectively (Bolatto et al. 2013). From Σ_H2 we then calculated the H₂ surface density.

All of the H i and H₂ maps we used have resolutions in the region of 1', being quite similar both to each other and to our Herschel maps—varying from 02 for the M33 H₂ data to 26 for the SMC H₂ data. All of the H i maps incorporate both interferometric and single-dish observations, ensuring that we will not miss any diffuse atomic structure due to lack of short-spacing data.

To make our maps of Σ_H, we convolved the H i and H₂ maps for each galaxy to whichever resolution was worse of the two maps, re-gridded them to the same projection, then added them together to produce the combined map of hydrogen surface density.

Dust and gas are typically well mixed in the ISM, across a wide range of surface densities (Hildebrand 1983; Eales et al. 2012; Williams et al. 2018). Given that ISM properties are strongly driven by density (see Section 1 and references therein), dust found at different locations, but sharing a given Σ_H, should have similar properties.

We convolved our Σ_H and foreground-subtracted Herschel maps to the limiting resolution for each galaxy. For the LMC and SMC, the limiting data were the Σ_H maps, with resolutions of 1' for the LMC and 26 for the SMC. For M31 and M33, the limiting data were Herschel at 500 μm, with its 36'' resolution. We then defined bins of Σ_H, each 0.025 dex wide. For each bin, all of the Σ_H pixels within that density range were identified, and all of the corresponding Herschel pixels had their mean surface brightnesses calculated, within each band.

For comparison, we also performed this process for our unfeathered maps—plus, for the LMC and SMC, the original HERITAGE maps (Meixner et al. 2013)—to allow us to evaluate the difference from using the newer reduction calibrations, pipelines, and feathering. For the purposes of an "apples-to-apples" comparison, we applied the polynomial-based background subtraction used for our unfeathered maps in Section 6.2 to the feathered and HERITAGE maps also. By comparing the results for the feathered maps when using the polynomial versus full-foreground subtractions, we are able to determine how much additional low-surface-brightness emission is retrievable thanks to the improved foreground subtraction that is only possible with the feathered maps.

In Figure 15 (LMC and SMC) and 16 (M31, and M33), we compare the relationship between H i surface density and surface brightness for the polynomial-subtracted feathered, unfeathered, and HERITAGE maps, along with the feathered maps that had undergone our full H i-based foreground subtraction. Note that it is hard to draw firm conclusions from differences between the HERITAGE reductions and our reductions, as HERITAGE used a much older version of HIPE (v7 for HERITAGE, compared to v12+ for our reductions), and there were some differences in calibration (e.g., nonlinearity corrections were added for PACS, bolometer relative gains were altered for SPIRE, and beam sizes were updated for both). However, it is still valuable to see what H column each data set is able to trace down to.

Zoom In Zoom Out Reset image size

We only performed the binning for this comparison for pixels that are covered in all of the data sets (Herschel, H i, and CO) for each given galaxy. As such, our ability to retrieve binned Herschel emission here is a lower limit (i.e., we would be able to probe even deeper if we were able to bin over the entire map area).

For the SMC (right of Figure 15) and M31 (left of Figure 16), we see particularly impressive gains with the new data. For the SMC, averaged across all five bands, the full-foreground-subtracted feathered data can detect dust emission to an H surface density a factor of 2 lower than the HERITAGE data can, and a factor of 2.3 lower than with the unfeathered data. Similarly, for M31, the feathered data allow us to detect dust emission to a mean factor of 5 lower in density (reaching a factor of 10 lower at 250 μm).

Zoom In Zoom Out Reset image size

For the LMC (left of Figure 15), we see that there are rather minimal differences between the surface brightnesses, down to the lowest H density for which we are able to trace dust emission. In each band, both the polynomial-foreground-subtracted HERITAGE data or full-foreground-subtracted feathered data were each able to trace down to similar H densities. This may be due to the fact that the LMC has a relatively abrupt "edge" to its ISM disk, and the fact that the Herschel map is relatively tight compared to this disk, with little sky around the edges. There are also smaller gains for M33 in most bands, only probing to H densities 25% lower, on average, than the unfeathered data (although the PACS bands specifically do go to 50% lower densities). As with the LMC, the M33 maps are comparatively small, with M33 lying close to the edges.

We also note that the full-foreground-subtracted feathered data also trace dust emission out to much lower H densities than the polynomial-foreground-subtracted feathered data. This well demonstrates the value of the fact that feathering the Herschel maps with all-sky survey maps allows the combined maps to properly trace the large-scale foreground cirrus—therefore making it possible to accurately subtract it by comparison to Galactic H i data.

For all galaxies and bands, our new data allow us to trace dust emission down to ISM surface densities of <10 M_⊙ pc⁻². It is interesting that the surface brightness versus H density relations do not drop off with decreasing H density; instead, the gradients remain the same, or even flatten out. This is suggestive of the dust abundance in these low-density environments—which will be the central focus of the second paper in this series (C. J. R. Clark et al. 2021, in preparation).

The magnitude and distribution of the emission restored by the feathering process varies from band to band for each galaxy. In general, our feathering process restored more diffuse emission in the shorter-wavelength data than at longer wavelengths (where there was much less filtering). This restoration was most prominent at the galaxies' outskirts. As a result, the emission at the edges of these galaxies is much bluer than it was in the unfeathered data.

This is demonstrated in Figure 17, which shows how much the S₅₀₀/S₁₆₀ color (i.e., the 500 μm to 160 μm surface brightness ratio) has changed between the unfeathered maps and the feathered maps. Bright compact regions show no change in color, as they suffered little filtering in the unfeathered data. However, diffuse regions, especially at the galaxies' peripheries, are much bluer in the feathered data, with S₅₀₀/S₁₆₀ often falling by 20%–30% over large areas. Note that this increase in blueness is despite the fact that the cross-calibration in Section 5.2 reduced brightnesses in the PACS bands by 25%–35% (see Table 5).

Zoom In Zoom Out Reset image size

Dust temperatures derived from SED fitting FIR–submm data are strongly driven by color (Bendo et al. 2012). As such, the increased blueness of the new data can be expected to affect dust temperatures—and also affect dust masses, and β (due to the temperature–β degeneracy; Kelly et al. 2012; Galliano 2018). Exploring the dust SED properties will be a central focus of the second paper in this series (C. J. R. Clark et al. 2021, in preparation).