MODELING DUST AND STARLIGHT IN GALAXIES OBSERVED BY SPITZER AND HERSCHEL: NGC 628 AND NGC 6946

Interstellar dust affects the appearance of galaxies, by attenuating short-wavelength radiation from stars and ionized gas, and contributing IR, submillimeter, mm, and microwave emission (for a review, see Draine 2003). Dust is also an important agent in the fluid dynamics, chemistry, heating, cooling, and even ionization balance in some interstellar regions, with a major role in the process of star formation. Despite the importance of dust, determination of the physical properties of interstellar dust grains has been a challenging task—even the overall amount of dust in other galaxies has often been very uncertain.

Many previous studies have used far-infrared photometry to estimate the dust properties of galaxies. For example, Draine et al. (2007) used global photometry of 65 galaxies in the "Spitzer Infrared Nearby Galaxies Survey" (SINGS) galaxies to estimate the total dust mass and polycyclic aromatic hydrocarbon (PAH) abundance in each galaxy, and to characterize the intensity of the starlight heating the dust. For most of these galaxies the photometry extended only to 160 μm, although ground-based global photometry at 850 μm was also available for 17 of the 65 galaxies.¹⁸ Muñoz-Mateos et al. (2009) used images of the SINGS galaxies to examine the radial distribution of the dust surface density. Sandstrom et al. (2010) studied the PAHs in the Small Magellanic Cloud (SMC) on a pixel-by-pixel basis with a very similar dust model as the present work. Very recently, Totani et al. (2011) used global photometry in six bands—9, 18, 65, 90, 140, and 160 μm—to estimate dust masses for a sample of more than 1600 galaxies in the Akari All Sky Survey. In the present work, we develop state-of-the-art image processing and dust modeling techniques aiming to reliably determine the dust properties in resolved studies.

The present study makes use of combined imaging by the Spitzer Space Telescope and Herschel Space Observatory, covering wavelengths from 3.6 μm to 500 μm, to produce well-resolved maps of the dust in two nearby galaxies. The present study is focused on two well-resolved spiral galaxies, NGC 628 and NGC 6946, as examples to illustrate the methodology. Subsequent work (G. Aniano et al. 2012, in preparation) will extend this analysis to all 61 galaxies in the KINGFISH sample (Kennicutt et al. 2011).

The paper is organized as follows. In Section 2 we discuss the data sources. Background subtraction and data processing are discussed in Section 3. The dust model is summarized in Section 4, and the fitting procedure is described in Section 5. Results for NGC 628 and NGC 6946 are given in Section 6. The sensitivity of the derived parameters to the set of cameras used to constrain the dust models is explored in Section 7, where we compare dust mass estimates obtained at high spatial resolution (without using Multiband Imaging Photometer (MIPS)160 or the longest wavelength Spectral and Photometric Imaging (SPIRE) bands) with estimates made at lower spatial resolution (using all the cameras available). We also investigate the reliability of the photometry by comparing MIPS70 and MIPS160 images with PACS70 and PACS160 images. In Section 8 we compare dust mass estimates based on spatially resolved images with dust mass estimates based on global photometry, as would apply to distant, unresolved galaxies. The "goodness of fit" of different dust models is discussed in Section 9. The principal results are discussed in Section 10 and summarized in Section 11.

Appendices A–D describe the method for image segmentation (i.e., galaxy and background recognition), background subtraction, estimation of photometric uncertainties in the images, and estimation of dust modeling uncertainties. Appendix E is a comparison of Photodetector Array Camera and Spectrometer for Herschel (PACS) and Multiband Imaging Photometer for Spitzer (MIPS) photometry.

NGC 628 and NGC 6946 are part of the SINGS galaxy sample (Kennicutt et al. 2003), and were imaged by the Spitzer Space Telescope (Werner et al. 2004). A large subset of the SINGS galaxies are also included in the Herschel Space Observatory Open Time Key Project KINGFISH (Kennicutt et al. 2011) and were observed with the Herschel Space Observatory (Pilbratt et al. 2010).

We will use "camera" to identify each optical configuration of the observing instruments, i.e., each different channel or filter arrangement of the instruments will be referred to as different "camera." With this nomenclature, each "camera" has a characteristic optical resolution, spectral response, and point-spread function (PSF).¹⁹ We will refer to the Infrared Array Camera (IRAC), MIPS, PACS, and SPIRE cameras using their nominal wavelengths in microns: IRAC3.6, IRAC4.5, IRAC5.8, IRAC8.0, MIPS24, MIPS70, MIPS160, PACS70, PACS100, PACS160, SPIRE250, SPIRE350, and SPIRE500. For our standard modeling, the PACS and SPIRE data were reduced by HIPE followed by Scanamorphos (see Section 2.2.1 for details). In Appendix E, and Table 6, where we compare PACS data with and without Scanamorphos processing, we use PACS(H)70, PACS(H)100, and PACS(H)160 to denote PACS data that were processed by the HIPE pipeline only.

Table 1 summarizes the optical resolution of the cameras.

2.1. Spitzer

2.2. Herschel

2.3. Atomic and Molecular Gas

NGC 628 and NGC 6946 are in The H i Nearby Galaxy Survey (THINGS; Walter et al. 2008), with ∼6'' resolution H i 21 cm imaging by the NRAO Very Large Array. The H i 21 cm emission is assumed to be optically thin, in which case the H i surface density is directly proportional to the 21 cm line intensity.

The CO J = 2  →  1 transition was observed with ∼13'' angular resolution by the HERA CO Line Extragalactic Survey (HERACLES; Leroy et al. 2009) using the HERA multi-pixel receiver on the IRAM 30 m telescope, with estimated uncertainties of ±20%.

As is usual, the H₂ mass surface density is taken to be proportional to the CO J = 2 → 1 line intensity. X_CO is the ratio of H₂ column density to CO J = 1 → 0 intensity integrated over the line profile. In what follows, we define

$$\begin{equation} X_{{\rm CO},20}\equiv \frac{X_{\rm CO}}{10^{20}\,{\rm H}_2\,{\rm cm}^{-2}\,({\rm K}\,{\rm km\,s}^{-1})^{-1}}, \end{equation} \tag{ 1 }$$

and we assume a J = 2 → 1/J = 1 → 0 antenna temperature ratio R₂₁ = 0.8 (Leroy et al. 2009).

The conversion factor X_CO is uncertain. X_{CO, 20} ≈ 2 is the value normally adopted for molecular gas in the Milky Way (MW; Dame et al. 2001; Okumura et al. 2009). Planck Collaboration et al. (2011b) found X_{CO, 20} = 2.54 ± 0.13 for the MW. Draine et al. (2007) found that X_{CO, 20} ≈ 4 appeared to give the most reasonable dust/H mass ratios for the SINGS galaxy sample. Blitz et al. (2007) found X_{CO, 20} ≈ 4 to be the best overall value for galaxies in the Local Group. Leroy et al. (2011) found X_{CO, 20} = 1.2–4.2 for M 31, M 33, and the Large Magellanic Cloud (LMC), and the very high values X_{CO, 20} = 14 and X_{CO, 20} = 32 for NGC 6822 and the SMC. In the present work, we take X_{CO, 20} = 4 ± 1 to be the best overall value for NGC 628 and NGC 6946. In Section 6 (Figure 3) we show the H gas maps and dust/H mass ratio for X_{CO, 20} = 2, 3, 4.

The H i and H₂ masses are added to generate maps of the total H surface density Σ_H. Including helium would give Σ_gas ≈ 1.38 × Σ_H, but in our discussion will use Σ_H as it is the "observable," avoiding uncertainties in the helium abundance.

Before modeling the spatially resolved spectral energy distributions (SEDs), it is necessary to adjust the images in a number of ways to ensure meaningful results. In the following sections we describe the image analysis steps in detail. These steps include background estimation and subtraction, convolution to a common PSF, and resampling of the convolved images to a common pixel grid. We also correct the final images for missing or bad data in the original images and estimate the uncertainties on the flux in each pixel. The pixel sizes used in the final maps are chosen to Nyquist sample the final-map PSFs.

All camera images are first rotated so that north is up, and then trimmed to a common sky region. We then estimate and remove a background from each image, as described in Appendix B. Following this, we convolve the images to a common PSF, and resample all the images on a common final-map grid. We correct the final images for the bad pixels (or missing data) in the original images, as described in Appendix C. Finally, we use the background pixel dispersion to estimate the pixel flux uncertainties, as described in Appendix D.

3.1. Background Recognition and Subtraction

We use a multistep algorithm to generate a Background Mask, consisting of the area not covered by either the target galaxy or other discrete sources (e.g., recognizable background galaxies or foreground stars). This algorithm is described in detail in Appendix A and avoids overestimating the background in cameras with low S/N.

After the background mask is generated, for each image we estimate the best-fit background "tilted plane," as described in Appendix A. After the final best-fit background "tilted plane" has been found for each camera, it is subtracted from the original images. The background subtraction is performed on each original map independently, using its native pixel grid.

The result is, for each camera, a background-subtracted image with its native pixel grid and PSF. Figures 1 and 2 show the resulting background-subtracted images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Convolution to a Common PSF

3.3. Uncertainty Estimation

The composition of interstellar dust remains uncertain, but models based on a mixture of amorphous silicate grains and carbonaceous grains have proven successful in reproducing the main observed properties of interstellar dust. We employ the dust model of Draine & Li (2007, hereafter DL07), using "MW" grain size distributions (Weingartner & Draine 2001, hereafter WD01). The DL07 dust model has a mixture of amorphous silicate grains and carbonaceous grains, with a distribution of grain sizes, chosen to reproduce the wavelength dependence of interstellar extinction within a few kpc of the Sun (Weingartner & Draine 2001). The silicate and carbonaceous content of the dust grains was constrained by observations of the gas phase depletions in the interstellar medium (ISM).

The bulk of the dust in the diffuse ISM is heated by a general diffuse radiation field contributed by many stars. However, some dust grains will happen to be located in regions close to luminous stars, such as photodissociation regions (PDRs) near OB stars, where the starlight heating the dust will be much more intense than the diffuse starlight illuminating the bulk of the grains. Since our pixels have a large physical size (≈500 pc side for MIPS160 PSF), we will assume that, in each pixel, there is dust exposed to a distribution of starlight intensities.

4.1. Carbonaceous Grains

The carbonaceous grains are assumed to have the properties of PAH molecules or clusters when the number of carbon atoms per grain N_C ≲ 10⁵, but to have the properties of graphite when N_C ≫ 10⁵, with an ad hoc smooth transition between the two regimes.

A carbonaceous particle of equivalent radius a is taken to have absorption cross section:

$$\begin{eqnarray} C_{\rm abs}(a,\lambda) &=& (1-f_{\rm g})N_{\rm C}\sigma _{\rm PAH}({\rm H:C},\lambda)\nonumber \\ &&+ f_{\rm g} C_{\rm abs}({\rm graphite},a,\lambda) , \end{eqnarray} \tag{ 2 }$$

where σ_PAH (H: C, λ) is the absorption cross section per C for PAH material with given H:C ratio, C_abs (graphite, a, λ) is the absorption cross section calculated for randomly oriented graphite spheres of radius a, and the graphite "weight" is taken to be

$$\begin{eqnarray} &&f_{\rm g} = 0.01\;\;{\rm for}\;\; a \le a_c \hspace{96pt} \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&= 0.01 + 0.99[1-(a_c/a)^3] \;\;{\rm for}\;\; a > a_c , \end{eqnarray} \tag{ 4 }$$

where we take the transition radius a_c = 0.0050 μm. Carbonaceous grains with a > a_c are, therefore, treated as having a graphitic component, with the graphite weight f_g → 1 for a ≫ a_c. However, with f_g = 0.01 for a < a_c, even very small PAHs are assumed to have a small continuum opacity underlying the PAH features.

The PAH absorption cross section can be represented as the sum of a number of vibrational features with specified central wavelength, FWHM, and band strength. Li & Draine (2001, hereafter LD01) presented a set of resonance parameters that appeared to be consistent with pre-Spitzer observations. Smith et al. (2007) used the PAHFIT fitting software with the SINGS spectra to improve observational determinations of central wavelengths, shapes, and overall strengths of the PAH emission profiles; DL07 used these results to adjust the LD01 profile parameters. Subsequent modeling of the SINGS nuclear spectra (ARD12) led to some additional small changes in some of the PAH band strengths. In the present study we employ the PAH cross sections from DL07 and ARD12.

Draine & Li (2007) adopted the model put forward by Draine & Lee (1984, hereafter DL84) for the far-infrared properties of graphite. Graphite is a highly anisotropic material, with very different responses for $\vec{E}\parallel \vec{c}$ and $\vec{E}\perp \vec{c}$, where the $\vec{c}$ axis is normal to the "basal plane." DL84 included "free-electron" contributions δ^f_⊥, δ^f_∥ to the dielectric tensor, using a simple Drude model for the free-electron response,

$$\begin{equation} \delta \epsilon ^f(\omega) = \frac{-(\omega _p\tau)^2}{(\omega \tau)^2+i\omega \tau }, \end{equation} \tag{ 5 }$$

where ω = 2πc/λ, and to allow for size effects the "mean free time" τ is taken to be

$$\begin{equation} \tau ^{-1} = \tau _{\rm bulk}^{-1} + v_{\rm F}/a, \end{equation} \tag{ 6 }$$

where τ_bulk is the mean free time in bulk material, v_F is the Fermi speed, and a is the grain radius.

For $\vec{E}\perp \vec{c}$ we continue to use the graphite dielectric function from DL84. However, for $\vec{E}\parallel \vec{c}$, the DL84 free-electron model for δ^f_∥ resulted in an opacity at λ ≳ 100 μm that gave somewhat more emission than observed in the MW cirrus by Finkbeiner et al. (1999). In addition, the free-electron model used by DL84 for $\vec{E}\parallel \vec{c}$ produced an opacity peak near 33 μm that does not give a good match to Infrared Spectrograph (IRS) observations of emission from regions where the grains are hot enough to radiate near 33 μm. To broaden the opacity peak, we now take the free-electron contribution for $\vec{E}\parallel \vec{c}$ to be

$$\begin{equation} \delta \epsilon _\parallel ^f(\omega) = - \sum _{j=1}^2 \frac{(\omega _{p,j}\tau _j)^2}{(\omega \tau _j)^2+i\omega \tau _j}, \end{equation} \tag{ 7 }$$

with ω_{p, 1} = 1 × 10¹⁴ s⁻¹, and ω_{p, 2} = 2 × 10¹⁴ s⁻¹. The τ_j are obtained from Equation (6) with τ_{bulk, 1} = 3.51 × 10⁻¹⁴ s, τ_{bulk, 2} = 0.88 × 10⁻¹⁴ s. This gives a d.c. electrical conductivity ∑_j(ω²_{p, j}τ_j)/4π = 5.62 × 10¹³ s⁻¹ = 62.5 mho cm⁻¹, within the range reported for high-quality graphite crystals at 300 K [∼1 mho cm⁻¹ (Klein 1962) to ∼200 mho cm⁻¹ (Primak 1956)]. This d.c. conductivity is larger than the value 30 mho cm⁻¹ adopted by DL84; the increased conductivity lowers the FIR emission and brings the overall emission spectrum into better agreement with the observed spectrum from Finkbeiner et al. (1999). We take v_F = 3.7 × 10⁶(1 + T/255 K)^1/2 cm s⁻¹. The two-component free-electron form of Equation (7) is not intended to have physical significance. It is adopted because it is analytic, satisfies the Kramers–Kronig relations, gives a reasonable value for the d.c. electrical conductivity, and results in an opacity that is less peaked than the original single-component form (5). The resulting graphite opacity varies as λ⁻² for λ ≳ 200 μm.

4.2. PAH Abundance q_PAH

4.3. Amorphous Silicate Grains

4.4. Dust Heating

Each dust grain is assumed to be heated by radiation with energy density per unit frequency:

$$\begin{equation} u_\nu =U\times u_\nu ^{{\rm MMP83}}, \end{equation} \tag{ 8 }$$

where U is a dimensionless scaling factor and u^MMP83_ν is the interstellar radiation field estimated by Mathis et al. (1983) for the solar neighborhood. We ignore variations in the spectral shape.

Each pixel in our modeling will be larger than 2 × 10⁴ pc², so it will contain ISM in a variety of physical environments. A fraction (1 − γ) of the dust mass is assumed to be heated by starlight with a single intensity U = U_min (i.e., heated by a diffuse ISM radiation field), while the remaining fraction γ of the dust mass is exposed to a power-law distribution of starlight intensities between U_min and U_max with dM/dU∝U^−α.

The starlight heating intensities are thus characterized by four parameters: γ, U_min, U_max, and α, where the fractional dust mass dM_d(U) heated by starlight intensities in (U, U + dU) is

$$\begin{eqnarray} \frac{1}{M_{d,{\rm tot}}} \left(\frac{{\it dM}_{{\rm d}}}{{\it dU}}\right) &=& (1-\gamma) \delta (U-U_{\rm min})+ \gamma \frac{(\alpha -1)}{U_{\rm min}^{1-\alpha }-U_{\rm max}^{1-\alpha }}U^{-\alpha } \nonumber \\ &&\quad\;{\rm for}\; U_{\rm min}\le U \le U_{\rm max} \nonumber \\ \end{eqnarray} \tag{ 9 }$$

for α ≠ 1, and

$$\begin{eqnarray} \frac{1}{M_{d,{\rm tot}}} \left(\frac{{\it dM}_{{\rm d}}}{{\it dU}}\right) &=& (1-\gamma) \delta (U-U_{\rm min})+ \gamma \frac{1}{ \ln \left(U_{\rm max}/U_{\rm min}\right)} U^{-1} \nonumber \\ &&\quad {\rm for}\; U_{\rm min}\le U \le U_{\rm max} \nonumber \\ \end{eqnarray} \tag{ 10 }$$

for α = 1, where $M_{d,{\rm tot}}\equiv {\int ^{U_{\rm max}}_{U_{\rm min}}(dM_{{\rm d}}/dU^{\prime })dU^{\prime }}$. More complicated starlight heating distributions could be contemplated, but we find that the simple four-parameter (γ, U_min, U_max, α) model of Equations (9) and (10) appears able to usually provide an acceptable fit to observed SEDs in star-forming galaxies with near-solar metallicities.

Galliano et al (2011) recently claimed that the emission from the LMC can be reproduced using a starlight distribution function that lacks the "delta function" component of Equations (9) and (10), i.e., fixed γ = 1. However, we show in Section 9 that adding the "delta function" component significantly improves the quality of the fit for the pixels in NGC 628 and NGC 6946.

Many authors choose to fit the λ ≳ 70 μm emission using a blackbody B_ν(T_d) multiplied by a power-law opacity ∝λ^−β. The best-fit value of T_d is closely related to our heating parameters U_min, γ, α. In a subsequent work (G. Aniano & B. T. Draine 2012, in preparation) we show that T_d ≈ 20 U^0.15_min K, when the DL07 SED is approximated by a blackbody multiplied by a power-law opacity.

Given that there may be significant regional variations in the starlight spectrum, U should be interpreted not as a measure of the starlight energy density, but rather as the ratio of the actual dust heating rate to the heating rate for the MMP83 radiation field.

The fraction of dust luminosity emerging in the PAH features does depend on the spectrum of the starlight heating the dust. Draine (2011a) showed that the fraction of the dust emission appearing at 8 μm increases by a factor of 1.57 as the starlight spectrum is changed from MMP83 to a 20 kK blackbody (cutoff at 13.6 eV). If the actual hν < 13.6 eV starlight spectrum is harder (softer) than the MMP83 spectrum assumed in the models, we will overestimate (underestimate) q_PAH.

4.5. Contribution of Direct Starlight

Starlight enters in the dust modeling in two ways: via dust heating (as discussed in Section 4.4) and as a direct starlight component (i.e., direct starlight escaping the observed region). Our main goal in the present work is to study the properties of the dust and the starlight heating the dust, so we adopt a simple model for the direct starlight component. Following Bendo et al. (2006) and Draine et al. (2007), we approximate the λ > 3 μm stellar emission from the galaxy as simply

$$\begin{equation} F_{\star }(\lambda) = \Omega _\star B_\nu (\lambda ,T_\star) , \end{equation} \tag{ 11 }$$

where Ω_⋆ is the solid angle subtended by the stars, B_ν is the blackbody function, and T_⋆ = 5000 K is a representative photospheric temperature to approximate the integrated stellar emission at λ > 3 μm. The direct starlight contribution is thus adjusted by only the parameter Ω_⋆. Direct starlight will only contribute significantly to the IRAC bands, and very marginally to MIPS24. In ARD12 corrections arising from photospheric absorption as well as emission from hot circumstellar dust around asymptotic giant branch stars are studied. These corrections are only important in elliptical galaxies with little ISM, and the results in the present paper would be virtually unchanged if included. We also neglect possible reddening at the wavelengths (λ ⩾ 3.6 μm) in the present study.

Although in principle the direct starlight and heating starlight parameters should be connected, the uncertain and complex distribution of dust and stars within the galaxies make such connection very complex.²⁷ Our simplified treatment of the direct starlight should only be regarded as a way of "removing" the direct starlight component from the near infrared photometry so we can have an estimate of the dust emission.

4.6. Dust Model Emission

NGC 628 and NGC 6946 are well resolved, with each galaxy providing many independent pixels, even at MIPS160 resolution. For each pixel j, we find the model of dust and starlight that best reproduces the observed SED, within the modeling scheme described by DL07.

As discussed in Section 4.5, starlight enters the fitting in two ways: via direct starlight in the pixel and by heating the dust. The direct starlight contribution to pixel j is adjusted by varying only the parameter Ω_{⋆, j}. The heating starlight intensity is characterized by four parameters: γ_j, U_{min, j}, U_{max, j}, and α_j, where the dust mass dM_d heated by starlight intensities in (U, U + dU) is given by Equation (9). For each pixel j we adjust the total dust mass M_{d, j}, the PAH abundance parameter q_{PAH, j} (PAH mass fraction), and the characteristics of the starlight heating the dust in that pixel. If mid-IR photometry is unavailable, one loses the ability to constrain q_PAH, but (adopting some arbitrary value of q_PAH for the modeling), the dust mass estimation itself would be largely unaffected. Fortunately, we do have mid-IR coverage of our galaxies, so we can obtain full q_PAH maps for them. The present modeling assumes the grains to be heated by a standard starlight spectrum (corresponding to the starlight in the local ISM), and to have a fixed balance between PAH neutrals and ions. When fitting to IRAC, MIPS, PACS, and SPIRE photometry, the best-fit value of q_PAH is then essentially proportional to the strength of the (nonstellar) IRAC8.0 band power relative to the total IR power.

We will find that γ_j ≪ 1 in nearly all regions where the dust luminosity surface density $\Sigma _{L_{\rm d},j}>10^7 \,L_{\odot }\,{\rm kpc}^{-2}$: here, U_{min, j} is presumed to represent the diffuse ISM, or the counterpart to the "infrared cirrus" component of the galaxy (Low et al. 1984), accounting for the bulk of the dust mass in pixel j. The small fraction γ_j of the dust mass exposed to starlight intensities U > U_{min, j} is presumed to correspond primarily to dust in star-forming regions.

The model flux density in camera k is

$$\begin{equation} F({\rm model},\lambda _k)=\langle F_{\star }(\lambda)\rangle _k +\frac{M_{\rm d}}{4\pi D^2}\langle p({\rm model})\rangle _k , \end{equation} \tag{ 12 }$$

where 〈F_⋆(λ)〉_k is the direct contribution of starlight given by Equation (11) convolved with the instrumental response function. The dust model is characterized by {Ω_⋆, M_d, q_PAH, γ, U_min , U_max , α}.

In principle, U_{max, j} could be treated as an adjustable parameter. Previous work (Draine et al. 2007) has shown that the quality of the global fit to the SED is relatively insensitive to the choice of U_max. We experimented by allowing U_{max, j} to be fitted²⁹ in the resolved maps, and found that the best value for most of the pixels is U_{max, j} = 10⁷. In the pixels where the resulting best-fit value is not 10⁷, fixing U_{max, j} = 10⁷ did not decrease the quality of the fit significantly, i.e., the total χ² is essentially the same, as shown in Section 9. Allowing U_max to be fitted in the range 10³ ⩽ U_max ⩽ 10⁷ or fixing it to U_max = 10⁷ does not produce appreciable changes in the inferred dust masses. We therefore fix U_{max, j} = 10⁷ in our modeling.

The limits on adjustable parameters are given in Table 2. The allowed range for U_min is determined by the wavelength coverage of the data used in the fit. For the SINGS galaxy sample, it was found that if the photometry extends to λ_max = 160 μm, models with U_min ⩾ 0.6 are well constrained. However, if longer wavelength data are available, we allow the possibility of cooler dust, heated by starlight intensities U < 0.6, down to U_min = 0.06 if λ_max = 250 μm, and down to U_min = 0.01 if λ_max ⩾ 350 μm.

We observe that for a given set of parameters {q_PAH, U_min, U_max, α} the model emission is multi-linear in {Ω_⋆, M_dust, γ}. This allows us to easily calculate the best values of {Ω_⋆, M_dust, γ} for a given parameter set {q_PAH, U_min, U_max, α}. Therefore, when looking for the best-fit model in the seven-dimensional model parameter space {Ω_⋆, M_d, q_PAH, γ, U_min U_max, α}, we only need to do a search over the four-dimensional subspace spanned by {q_PAH, U_min, U_max, α}. In the case of a fixed U_max, the search is performed over a three-dimensional space spanned by {q_PAH, U_min, α}. In any case, for the computed grid of {q_PAH, U_min, U_max, α}, the multidimensional search for optimal parameters can be performed by brute force, rather than needing to rely on a nonlinear minimization algorithm.

With U_max fixed, for each pixel j, the model library is used to search for the model parameter vector ξ_j = {Ω_⋆, M_d, q_PAH, γ, U_min, α} that minimizes

$$\begin{equation} \chi _j^2 \equiv \sum _{k} \frac{[F_{{\rm obs},j}(\lambda _k)_j-F({\rm model},\lambda _k)]^2}{\sigma _{\lambda _k,j}^2} , \end{equation} \tag{ 13 }$$

where F_{obs, j}(λ_k) is the observed flux density and $\sigma _{\lambda _k,j}$ is the 1σ uncertainty in the measured flux density for pixel j at wavelength λ_k (see Appendix D for a detailed discussion on how $\sigma _{\lambda _k,j}$ is obtained).³⁰

The above procedure yields "best-fit" estimates for the model parameter vector ξ_j = {Ω_⋆, M_d, q_PAH, γ, U_min, α} for each pixel j. Each pixel is fitted independently of the remaining pixels in the galaxy. Since the final-map pixel size is chosen to Nyquist sample the final-map PSF, the camera images are smooth on a pixel scale. The fact that we obtain smooth parameter maps is an indication of the stability of the fitting procedure, i.e., even though every pixel is modeled independently of its neighbors, the continuity in the images leads to continuity in the results. The quoted "best-fit" parameter maps arise from fitting the observed flux in each pixel.³¹

5.1. Properties Derived from the Models

The infrared luminosity L_{d, j} for the dust model in the pixel j is

$$\begin{equation} L_{{\rm d},j} = P_0(q_{{\rm PAH},j}) \,M_{{\rm d},j}\, \overline{U}_j\,, \end{equation} \tag{ 14 }$$

where P₀(q_PAH) (which depends only weakly on q_PAH) is the total power radiated per unit dust mass by the model when heated by starlight with intensity U = 1, and $\overline{U}_j$ is the mass-weighted mean starlight heating intensity, given by

$$\begin{eqnarray} &&\overline{U}_j = (1-\gamma _j) U_{{\rm min},j} + \gamma _j \nonumber \\ &&\quad\times\left\lbrace \begin{array}{@{}l l}\left(\dsty\frac{\alpha _j-1}{\alpha _j-2}\right) \left[\dsty\frac{U_{\rm max}^{2-\alpha _j}-U_{{\rm min},j}^{2-\alpha _j}}{U_{\rm max}^{1-\alpha _j}-U_{{\rm min},j}^{1-\alpha _j}} \right] &{\rm for}\;\alpha _j\ne 1,\alpha _j\ne 2\\\ms U_{\rm max}\, \dsty\frac{1-U_{{\rm min},j}/U_{\rm max}}{\ln \left(U_{\rm max}/U_{{\rm min},j}\right)} &{\rm for}\; \alpha _j = 1\\\ms U_{{\rm min},j}\, \dsty\frac{\ln \left(U_{\rm max}/U_{{\rm min},j}\right)}{1-U_{{\rm min},j}/U_{\rm max}} &{\rm for}\;\alpha _j=2 \end{array}{.} \right.\nonumber\\ \end{eqnarray} \tag{ 15 }$$

Star-forming regions have significant starlight power absorbed by dust grains in regions of high starlight intensity, which generally correspond to PDRs. We will refer to the luminosity radiated by dust in regions with U > U_PDR as L_PDR, given by

$$\begin{eqnarray} &&L_{{\rm PDR},j} = P_0(q_{{\rm PAH},j}) \,M_{{\rm d},j} \,\gamma _j \nonumber \\ &&\quad \times \left\lbrace \begin{array}{@{}l l}\left(\dsty\frac{\alpha _j-1}{\alpha _j-2}\right) \left[\dsty\frac{U_{\rm max}^{2-\alpha _j}-U_{\rm PDR}^{2-\alpha _j}}{U_{\rm max}^{1-\alpha _j}-U_{{\rm min},j}^{1-\alpha _j}} \right] &{\rm for}\;\alpha _j\ne 1,\alpha _j\ne 2\\\ms\ms \dsty\frac{U_{\rm max}-U_{\rm PDR}}{\ln \left(U_{\rm max}/U_{{\rm min},j}\right)} &{\rm for}\;\alpha _j = 1\\\ms\ms \dsty\frac{\ln \left(U_{\rm max}/U_{\rm PDR}\right)}{U_{{\rm min},j}^{-1}-U_{\rm max}^{-1}} &{\rm for}\;\alpha _j=2\\ \end{array}. \right. \nonumber \\ \end{eqnarray} \tag{ 16 }$$

We take U_PDR = 10² as a plausible cutoff to select dust in high intensity regions (choosing another cutoff value would change only the inferred L_{PDR, j}, leaving all the remaining dust parameters unaltered). We further define f_{PDR, j} as

$$\begin{equation} f_{{\rm PDR},j} \equiv {L_{{\rm PDR},j} \over L_{{\rm d},j}}\,. \end{equation} \tag{ 17 }$$

The region observed is at a distance D from the observer and Ω_j is the solid angle of pixel j. For each pixel j, the best-fit model vector {Ω_⋆, M_d, q_PAH, γ, U_min, α}_j corresponds to a dust mass surface density:

$$\begin{equation} \Sigma _{M_{\rm d},j} \equiv \frac{1}{D^2\Omega _j} M_{{\rm d},j}. \end{equation} \tag{ 18 }$$

Similarly, we can compute the infrared luminosity surface density $\Sigma _{L_{\rm d},j}$ and the surface density of dust luminosity from regions with U > U_PDR, $\Sigma _{L_{\rm PDR},j}$, as

$$\begin{equation} \Sigma _{L_{\rm d},j} \equiv \frac{1}{D^2\Omega _j} L_{{\rm d},j}, \quad \quad \Sigma _{L_{\rm PDR},j} \equiv \frac{1}{D^2\Omega _j} L_{{\rm PDR},j}. \end{equation} \tag{ 19 }$$

The DL07 dust models used here are consistent with the MW ratio of visual extinction to H column, A_V/N_H = 5.34 × 10⁻²² mag cm²/H, for a dust/H ratio $\Sigma _{M_{\rm d}}/N_{\rm H}m_{\rm H}= 0.010$ (see Table 3 of Draine et al. 2007). The dust surface density corresponds to a visual extinction (through the disk):

$$\begin{equation} A_V = 0.67 \left(\frac{\Sigma _{M_{\rm d}}}{10^5\,M_{\odot }\,{\rm kpc}^{-2}}\right) \,{\rm mag}. \end{equation} \tag{ 20 }$$

5.2. Global Quantities

After the resolved (pixel-by-pixel) modeling of the galaxy is performed, we compute a set of global quantities by adding or taking weighted means (denoted as 〈...〉) of the quantities in each individual pixel of the map. The total dust mass M_{d, tot}, total dust luminosity L_{d, tot}, and total dust luminosity radiated by dust in regions with U > U_PDR, L_{d, tot}, are given by

$$\begin{eqnarray} M_{{\rm d},{\rm tot}} &\equiv & \sum _{j=1}^N M_{{\rm d},j} ,\quad L_{{\rm d},{\rm tot}} \equiv \sum _{j=1}^N L_{{\rm d},j} ,\nonumber \\ L_{{\rm PDR},{\rm tot}} &\equiv & \sum _{j=1}^N L_{{\rm PDR},j} , \end{eqnarray} \tag{ 21 }$$

where the sums extend over all the pixels j that correspond to the target galaxy (see Appendix A for the galaxy segmentation procedure). The dust-mass-weighted PAH mass fraction 〈q_PAH〉, and mean starlight intensity $\langle \overline{U} \rangle$, are given by

$$\begin{eqnarray} &&\langle q_{\rm PAH}\rangle \equiv \frac{\sum _{j=1}^N q_{{\rm PAH},j}\, M_{{\rm d},j}}{\sum _{j=1}^N M_{{\rm d},j}} ,\quad \langle \overline{U}\rangle \equiv \frac{\sum _{j=1}^N \overline{U}_j \, M_{{\rm d},j}}{\sum _{j=1}^N M_{{\rm d},j}} . \nonumber \\ \end{eqnarray} \tag{ 22 }$$

The dust mass-weighted minimum starlight intensity 〈U_min〉 is given by

$$\begin{equation} \langle U_{\rm min}\rangle \equiv \frac{\sum _{j=1}^N (1-\gamma _j)\, U_{{\rm min},j}\,M_{{\rm d},j}}{\sum _{j=1}^N (1-\gamma _j)M_{{\rm d},j}} . \end{equation} \tag{ 23 }$$

The dust-luminosity-weighted value of f_PDR, 〈f_PDR〉, is

$$\begin{equation} \langle f_{\rm PDR}\rangle \equiv \frac{L_{{\rm PDR},{\rm tot}}}{L_{{\rm d},{\rm tot}}} . \end{equation} \tag{ 24 }$$

Alternatively, we can fit the dust model to the global photometry for each galaxy³² (i.e., a single-pixel dust model). The dust parameters obtained from the global-photometry (single-pixel) model fit will be compared to the corresponding resolved modeling global quantities defined in Equations (21)–(24) in Section 8.

Clearly, the emission summed over the map will not correspond to the emission of dust exposed to starlight of the form given by Equations (9) and (10) since different pixels will have different values of U_{min, j}, γ_j, and α_j. For completeness, we could define the dust mass-weighted mass fraction heated by a power-law (U > U_min) component 〈γ〉 as

$$\begin{equation} \langle \gamma \rangle \equiv \frac{\sum _{j=1}^N \gamma _j \, M_{{\rm d},j}}{\sum _{j=1}^N M_{{\rm d},j}}, \end{equation} \tag{ 25 }$$

and an "effective power-law exponent" 〈α〉, defined as the value that satisfies

$$\begin{eqnarray} &&L_{{\rm PDR},{\rm tot}} = \frac{ L_{{\rm d},{\rm tot}}}{\langle \overline{U}\rangle } \, \langle \gamma \rangle \, \left(\frac{\langle \alpha \rangle -1}{\langle \alpha \rangle -2}\right) \left(\frac{U_{\rm max}^{2-\langle \alpha \rangle }-U_{\rm PDR}^{2-\langle \alpha \rangle }}{U_{\rm max}^{1-\langle \alpha \rangle }-\langle U_{\rm min}\rangle ^{1-\langle \alpha \rangle }} \right). \nonumber \\ \end{eqnarray} \tag{ 26 }$$

Unfortunately, 〈γ〉 and 〈α〉 do not turn out to be useful global quantities, and will, in general, differ significantly from the values obtained from the global-photometry (single-pixel) model fit.

5.3. Parameter Uncertainty Estimates

To estimate uncertainties in the derived dust parameters for each pixel j, we simulate data by adding zero-mean random noise δF(λ_k)_{j, r}, for r = 1, 2, 3, ..., N_r, to the observed flux F(λ_k)_j in each band, and fit the simulated noisy data. In Appendix D we describe the statistical construction of the sample F(λ_k)_{j, r} of N_r random realizations. For the fit parameters {a, b, ...} ∈ {Ω_⋆, M_d, q_PAH, γ, U_min, α} we have a set of N_r + 1 values and we calculate the covariance matrix:

$$\begin{equation} V_{ab} = V_{ba} \equiv N_r^{-1}\sum _{r=1}^{N_r}(a_r-a_0)(b_r-b_0) , \end{equation} \tag{ 27 }$$

where a₀, b₀ are the best-fit parameter values for the observed fluxes, and a_r, b_r are the best-fit values for the rth random noise realization. For each model parameter a, the 1σ uncertainty is taken to be

$$\begin{equation} \sigma (a) = V_{aa}^{1/2}. \end{equation} \tag{ 28 }$$

Each random noise realization r = 1, 2, ...N_r produces a global quantity via Equations (21)–(24) and (26). We proceed to calculate uncertainties for global quantities using Equations (27) and (28). In Appendix E we describe the details of the parameter uncertainty estimation.

We construct dust maps with several different angular resolutions. For a given resolution, one can only use cameras with the PSF FWHM smaller than the reference PSF. In a normal star-forming galaxy, most of the dust mass is at temperatures T_d ≈ 15–25 K, with νL_ν peaking at λ ≈ hc/6kT_d ≈ 100–160 μm. Reliable estimates of the dust mass therefore require long-wavelength data, at least out to 160 μm. In the present study we will compare dust models using maps at the resolution of the PACS160 (FWHM = 112), SPIRE250 (FWHM = 182), and MIPS160 (FWHM = 388) cameras. As discussed by Aniano et al. (2011) the MIPS160 PSF cannot be convolved safely into the SPIRE500 PSF (with FWHM = 361). This implies that MIPS160 is the narrowest PSF that can be used if we want to include MIPS160 data into the modeling. A drawback of using the MIPS160 PSF is its extended wings, which can cause power radiated by the bright central regions to contribute significantly to the observed surface brightness of the faint outer regions. However, we will see in Section 7 that this effect does not significantly affect dust mass estimates, and discrepancies between PACS and MIPS photometry make it important to include the MIPS160 camera in the dust modeling.

In order to generate the H surface density maps (used in the dust mass/H mass ratio maps), a value of X_{CO, 20} needs to be chosen. Figure 3 shows the total H maps (first and third rows) for both galaxies, for X_{CO, 20} = 2, 3, and 4, convolved to the MIPS160 PSF. Note that the CO maps cover almost all of both galaxy masks, but do not cover the full field of view, leading to the box-like step in the H mass maps.

Zoom In Zoom Out Reset image size

Using dust maps based on all of the available photometry (with the MIPS160 PSF—see Figures 4(c) and 9(c) below), Figure 3 shows the dust/H mass ratios for the different values of X_{CO, 20}, for NGC 628 (second row) and NGC 6946 (fourth row). The choice X_{CO, 20} = 4 gives the smoothest dust/H mass ratios over the galaxies (outside of the central region of NGC 6946, which is further discussed in Section 6.2.1). As discussed in Section 2.3, we take X_{CO, 20} = 4 for both NGC 628 and NGC 6946.

6.1. NGC 628

6.1.1. Maps of Gas and Dust

NGC 628 (= M 74), at a distance D = 7.2 Mpc, is classified as SAc. With major and minor optical diameters of 10.5 and 9.5 arcmin, it is well resolved even by the MIPS160 camera.

The star formation rate is estimated to be 0.7 ± 0.2 M_☉ yr⁻¹ (Calzetti et al. 2010). Two supernovae have been observed in NGC 628: SN 2002ap (Type Ic) and SN 2003gd (Type II-P).

The total H i mass from 21 cm observations is $M({\rm H\,\mathsc{i}})=(3.7\pm 0.2)\times 10^9\,M_{\odot }$ (Walter et al. 2008) and the total H₂ mass estimated from observations of CO 2–1 is M(H₂) = (2.5 ± 0.6) × 10⁹(X_{CO, 20}/4) M_☉ (Leroy et al. 2009). The fact that the adopted value of X_CO is larger than the value X_{CO, 20} ≈ 2 found for resolved CO clouds in the MW will be discussed in Section 10 below.

Figure 4 shows the resulting maps of H surface density (panels (a)–(c)), dust surface density $\Sigma _{M_{\rm d}}$ (panels (d)–(f)), dust luminosity surface density $\Sigma _{L_{\rm d}}$ (panels (g)–(i)), and dust/H mass ratio (panels (j)–(l)), obtained by fitting photometry with PACS160, SPIRE250, and MIPS160 resolution, using all the compatible cameras in each case. The dust models with PACS160 resolution show clear spiral structure, but the noise in the smaller pixels is such that the dust is only reliably detected at surface densities $\Sigma _{M_{\rm d}}\gtrsim 10^{5.0}\,M_{\odot }\,{\rm kpc}^{-2}$, corresponding to A_V ≳ 0.7 mag. If the mapping is done at the resolution of the SPIRE250 camera, the dust is reliably detected for $\Sigma _{M_{\rm d}} \gtrsim 10^{4.5}\,M_{\odot }\,{\rm kpc}^{-2}$, corresponding to A_V ≳ 0.2 mag; at this resolution, the spiral structure is still visible. Maps made at MIPS160 resolution are the most sensitive, because of the larger pixel size, and the fact that they use data from all of the cameras, including the MIPS160 camera. Unfortunately, the lower resolution of these maps (36'' FWHM) largely washes out the spiral structure which is visible in higher resolution maps of NGC 628. Nevertheless, imaging at MIPS160 resolution allows reliable detection of dust at surface densities as low as 10^4.3 M_☉ kpc⁻², or A_V ≳ 0.14 mag.

Zoom In Zoom Out Reset image size

We note that the dust/H mass ratio maps change significantly when the modeling is done at the different resolution/camera combination (see the last row of Figure 4). These discrepancies are mainly due to the low sensitivity of the PACS cameras in the low-surface brightness areas (compare the MIPS and PACS flux in the outer parts of the galaxies in Figures 1 and 2), and discrepancies in the high surface brightness areas. In Section 7 we discuss the PACS–MIPS photometry discrepancies further.

6.1.2. Maps of q_PAH and Starlight Parameters

Figure 5 shows maps of the PAH abundance q_PAH (panels (a)–(c)), the starlight intensity parameter U_min (panels (d)–(f)), and the PDR fraction f_PDR (panels (g)–(i)), over the "galaxy mask" region where the galaxy is well detected. The GM for NGC 628 has a diameter of ∼016, or 20 kpc at 7.2 Mpc.

Zoom In Zoom Out Reset image size

The PAH abundance parameter q_PAH, shown in Figures 5(a)–(c), is remarkably uniform over the region where it can be reliably estimated. In Figure 5(c), q_PAH varies from a high of ∼0.05 a few kpc from the center down to ∼0.035 near the edge of the galaxy mask. If there is a radial gradient in q_PAH, it is weak, consistent with the weak gradients found for the SINGS sample (including NGC 628) by Muñoz-Mateos et al. (2009).

At PACS160 resolution, the starlight intensity parameter U_min (see Figure 5(d)) varies between ∼0.6 and ∼3 over most of NGC 628, following the galaxy structures, but near the edges of the galaxy mask U_min appears to rise. This is because the reduced S/N results in PACS70/PACS160 ratios that appear to be anomalously high, leading to high inferred values of U_min for some pixels. This is probably the result of low PACS160 fluxes for those pixels, making it appear that the dust is rather warm. We see that when SPIRE250 data are introduced (Figure 5(e)), we have many fewer high values of U_min near the edge of the galaxy mask, and the U_min values in the brighter regions appear well behaved. This continues when the MIPS160, SPIRE350, and SPIRE500 data are brought into the fit in the MIPS160 resolution image (Figure 5(f)).

The bottom row of Figure 5 (panels (g)–(i)) shows f_PDR, the fraction of the dust luminosity coming from dust heated by starlight intensities U > 100, which we expect to be associated with star-forming regions. At PACS160 resolution, most of the galaxy has f_PDR ≈ 0.03 and some small bright regions (for example, the center of aperture 2 in Figure 8) have higher values, up to f_PDR ≈ 0.3. These regions with high values of f_PDR generally coincide with Hα peaks, and also with the regions of the highest dust luminosity per area (see Figure 4(g)). As we shift to coarser modeling (i.e., using PSFs with larger FWHM), the f_PDR peaks are smoothed, and the general galaxy pixels tend to have larger f_PDR values: the dynamic range of values of f_PDR decreases. At MIPS160 resolution, most of the galaxy has f_PDR > 0.05. The overall (dust-luminosity-weighted) mean for the galaxy is 〈f_PDR〉 = 0.116: 11.6% of the total IR power is radiated by dust in regions where U > 100.

Unfortunately, at MIPS160 resolution the arm structure of both galaxies is not clearly resolved (see the MIPS160 images of the galaxies in panel (f) of Figures 1 and 2). Therefore, we cannot reliably study variation of the dust parameters between arm and interarm regions. A subsequent work (L. K. Hunt et al. 2012, in preparation) will study radial variations in the model parameters.

6.1.3. Comparison between Observed and Modeled Flux Densities

How well does the dust model reproduce the SPIRE photometry? Figure 6 shows the ratios of model-predicted intensity to observed intensity at λ = 250, 350 and 500 μm for dust models obtained by fitting photometry with PACS160, SPIRE250, and MIPS160 resolution. In order to make the comparison, we degrade either the observed image or the model-predicted image to a common resolution (i.e., when the modeling is done at PACS160 resolution, we convolve the model-predicted SPIRE250 image to the SPIRE250 PSF, and when the modeling is done at MIPS160 resolution, we convolve the observed SPIRE250 image to the MIPS160 PSF). Except near the edge of the galaxy mask (where the low S/N in the PACS data becomes an issue), the modeling tends to overpredict the SPIRE500 photometry—there is no evidence for a significant mass of very cold dust radiating only at the longer wavelengths.

Zoom In Zoom Out Reset image size

When we attempt to model at PACS160 resolution (using only IRAC, MIPS24, and PACS data to constrain the model), the model predictions at SPIRE250, SPIRE350, and SPIRE500 do not agree very well with observations (see Figures 6(a), (d), and (g)). If we coarsen the modeling to SPIRE250 resolution and add SPIRE250 to the model constraints, we now reproduce the SPIRE250 image (not surprising) and the model predictions at 350 μm and 500 μm are within ±30% and ±50% of the SPIRE observations, respectively, over most of the galaxy mask. If we further coarsen the modeling to the MIPS160 PSF, and use all the data to constrain the model, we find good agreement at all SPIRE bands (see Figures 6(c), (f), and (i)), with only SPIRE500 slightly overpredicted by ≈10% in some regions.

Figures 6(c), (f), and (i) show that when all of the cameras are used to constrain the modeling (and with the improved S/N of the larger MIPS160 pixels), the model emission is generally in good agreement with the SPIRE imaging—for all three SPIRE bands the ratio of model/observation is close to 1 over most of the galaxy, with significant departures only just at the edge of the galaxy mask, where observational errors are likely to be the explanation.

6.1.4. Global SED Fitting

Figure 7 shows global SEDs for NGC 628. The observed photometry is represented by rectangular boxes (Spitzer (IRAC, MIPS) in red; Herschel (PACS(S), SPIRE) in blue). The model convolved with camera response function (i.e., expected camera photometry for the model) is represented by diamonds (Spitzer in red, Herschel in blue). The two rows show the global SED for NGC 628. The observed global photometry for NGC 628 is the same across the two top rows, but the models differ. Each black curve is a fitted model spectrum. The diamonds show the fitted model spectrum convolved with the instrumental response function for each camera, i.e., the expected photometry for the fitted model.

Zoom In Zoom Out Reset image size

The top row shows "single-pixel" models for the galaxy based on (a) IRAC, MIPS24, and PACS data only, (b) IRAC, MIPS24, PACS, and SPIRE250, and (c) IRAC, MIPS, PACS, and SPIRE (all the cameras).

The second row shows the predicted model SED obtained by fitting a dust model to each resolved pixel, and then summing the emission over all the pixels. These predicted SEDs differ from the previous "single pixel" predicted SED because the dust modeling is a nonlinear process. The total dust mass predicted in this way is, however, similar to the "single-pixel" prediction (see Section 8 for details). When the modeling is done at lower resolutions (MIPS160) the "single pixel" and map-averaged quantities are in closer agreement.

6.1.5. Fitting in Selected Apertures

Figure 8 shows SEDs for four circular apertures located on NGC 628, sampling a wide range of surface brightnesses, $\Sigma _{L_{\rm d}}\approx (0.3\hbox{--}9)\times 10^{7}\,L_{\odot }\,{\rm kpc}^{-2}$ (see Figure 4 for aperture locations). Aperture 4 is located partially outside the galaxy mask, where the single-pixel dust modeling is not reliable, but the improved S/N of a large aperture allows reliable determination of the dust model parameters. We fit the DL07 model to the summed flux within each aperture. The top row shows the SED for a 60'' (diameter) circular aperture centered on the galaxy nucleus. The second row shows the SED for a 60'' aperture located on a bright spot on the spiral arms. We observe that the PACS photometry is larger than the corresponding MIPS photometry in the high surface brightness apertures 1 and 2. When both data sets are included in the modeling, the dust model gives values intermediate between PACS and MIPS. The third and fourth rows show SEDs in 80'' and 120'' apertures further from the center. We employ larger apertures in order to obtain reasonable S/N in these fainter regions. In aperture 4, the MIPS and PACS photometry differ significantly. When SPIRE, MIPS70, and MIPS160 are not used (e.g., in the PACS160 resolution modeling, Figure 8(j)) the high PACS70/PACS160 ratio causes the model to infer very high values of U_min and low values of $\Sigma _{M_{\rm d}}$, and hence to underpredict the SPIRE photometry. When we include SPIRE250 (Figure 8(k)) the model can reproduce SPIRE250, but continues to underpredict SPIRE350 and SPIRE500. Finally, when SPIRE350,500 and MIPS70,160 are added as constraints (Figure 8(l)), the modeling improves dramatically in aperture 4 and other low-brightness areas, reproducing most of the data, and making it clear that PACS70 is an outlier.

Zoom In Zoom Out Reset image size

The astute reader will note that the estimated uncertainties for, e.g., the PACS100 global photometry differ between the columns (i.e., for fluxes extracted after convolving to different PSFs). For each PSF, we estimate the noise per pixel based on the pixel statistics in the background region (see Appendix D). We then make a simple assumption concerning the pixel-to-pixel noise correlation. The fact that the uncertainty estimates for the aperture fluxes depend on the PSF is an indication that our assumption about the correlated component of the noise is imperfect. This is only an issue for the faint, low S/N data, such as the PACS fluxes in apertures 3 and 4.

In aperture 4 (Figure 8, last row), the model uses γ = 1, i.e., the dust is heated almost entirely by a power-law U distribution, but with very high values of γ. This corresponds to a broad distribution of starlight intensities in the U ≳ U_min range, with very little power in the high intensity range (PDR). This may be an artifact arising from the large photometric uncertainties.

6.2. NGC 6946

NGC 6946, an active star-forming galaxy, is classified as SABcd. At a distance D = 6.8 Mpc, the total H i mass is $M({\rm H\,\,\mathsc{i}})=(5.5\pm 0.3)\times 10^9\,M_{\odot }$ (Walter et al. 2008), with 42% of the H i falling within the galaxy mask. The H₂ mass within the galaxy mask is M(H₂) = (8.6  ±  0.6) × 10⁹(X_{CO, 20}/4) M_☉ (Leroy et al. 2009). As discussed previously, we adopt X_{CO, 20} = 4 as a global estimate. The star formation rate is estimated to be 4.5 M_☉ yr⁻¹ (Calzetti et al. 2010). NGC 6946 is remarkable for hosting at least nine supernovae over the past century (Prieto et al. 2008).

The variable carbon star V0778 Cyg, located near R.A. = 309.044, decl. = 60.082 (slightly off the bottom left corner of the maps shown in Figures 9–11), saturates the IRAC and MIPS detectors. The associated image artifacts affect the background estimation and subtraction, making the modeling less reliable in the bottom left corner of the maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.2.1. Dust and H Mass Maps

Figure 9 shows maps of the gas in NGC 6946 obtained from the THINGS 21 cm map (Walter et al. 2008) and the HERACLES CO 2–1 map (Leroy et al. 2009), convolved to the resolution of PACS160, SPIRE250, and MIPS160. The second row shows dust maps obtained by the present study. The dust map obtained using the PACS160 PSF has dust clearly visible only out to a radius of ∼250'' (8 kpc at 6.8 Mpc), but the 11'' (500 pc) FWHM of the beam resolves the spiral structure. The arm/interarm contrast in dust surface density appears to be approximately a factor of ∼2 in this image. Introducing SPIRE 250 μm data requires the PSF to be broadened, making the spiral arms less apparent, but allows dust to be detected out to larger radii, notably in the northern spiral arm. The MIPS160 resolution image clearly shows dust in the northern spiral arm out to a distance of ∼15 kpc from the nucleus.

The dust luminosity/area $\Sigma _{L_{\rm d}}$ (Figure 9, bottom row) peaks at ∼10^10.1 L_☉ kpc⁻² in the nucleus (see the PACS160 resolution map), and can be followed down to $\Sigma _{L_{\rm d}}\approx 10^{7.0}\,L_{\odot }\,{\rm kpc}^{-2}$ in the MIPS160 resolution map. Similarly, the IR luminosity/area contributed by PDRs peaks at $\Sigma _{L_{\rm PDR}}\approx 10^9\,L_{\odot }\,{\rm kpc}^{-2}$ near the center, and is reliably measured down to ∼10^5.6 L_☉ kpc⁻².

The dust/H mass ratio shown in Figures 9(j)–(l), calculated assuming X_{CO, 20} = 4, shows a pronounced minimum of ∼0.005 in the central  kpc. In this region the gas is primarily molecular, and the estimated dust/H mass ratio is therefore sensitive to the value assumed for X_CO. NGC 6946 has approximately solar metallicity, and we expect the dust/H mass ratio to be ∼0.01. The surprisingly low dust/H mass ratios found in the center of NGC 6946 indicate that the value of X_CO near the center should be about a factor of ∼2 smaller than the value X_{CO, 20} = 4 adopted for the bulk of the galaxy. This conclusion is consistent with interferometric studies of giant molecular clouds (GMCs) which indicate X_{CO, 20} ≈ 1.25 near the center of NGC 6946 (Donovan Meyer et al. 2012), using virial mass estimates for individual GMCs.

6.2.2. Maps of q_PAH and Starlight Parameters

The PAH abundance q_PAH in NGC 6946 rises from a minimum of ≲ 1% near the nucleus, reaching levels of ∼(4 ± 1)% at distances ∼1–8 kpc from the center, ultimately appearing to decline at the edge of the detectable region. The central minimum in q_PAH is real (it is independent of uncertainties in the appropriate value of X_CO). The local minimum of q_PAH at the center may be the result of destruction of PAHs by processes associated with star formation (there is no evidence of active galactic nucleus activity in NGC 6946). We note that there are a number of other local minima of q_PAH evident in Figures 10(a)–(c); just as for the central minimum, these extranuclear minima tend to coincide with peaks in dust luminosity surface brightness $\Sigma _{L_{\rm d}}$ (see Figures 9(g) and (h)) and peaks in f_PDR (see Figures 10(j), (k), and (l)). These are both indicators of luminous star-forming regions, which can be expected to coincide with H ii regions around hot stars, which appear to destroy PAHs (e.g., Povich et al. 2007). Supernova blast waves are presumed to also destroy PAHs.

Because the outer falloff in q_PAH is occurring as the S/N is falling to low values, it is not certain whether the observed decline is real or an artifact of different levels of background subtraction at different wavelengths. However, the decline in q_PAH in the outer regions persists in the MIPS160 resolution map. These variations in q_PAH will be examined in more detail in future studies.

The diffuse starlight intensity parameter U_min (Figure 10, middle row) shows a general decline with increasing distance from the center. The PACS160 resolution map of U_min is somewhat noisy, particularly in the outer regions, but the MIPS160 resolution map of U_min shows a systematic decline from values as high as U_min ≈ 10 in the central ∼500 pc down to U_min ≈ 0.15 in the outer regions, ∼10 kpc from the center. Even lower values of U_min are estimated for some pixels, but this is seen only at the lowest S/N levels.

The bottom row of Figure 10 shows f_PDR. The behavior of f_PDR in NGC 6946 is similar to NGC 628. Bright complexes coincide with local maxima of f_PDR. Outside very bright complexes f_PDR is quite smooth across the galaxy.

6.2.3. Comparison between Observed and Modeled Flux Densities

The model-predicted/observed flux ratios for NGC 6946 (shown in Figure 11) behave similarly to those of NGC 628 (see Figure 6). When long-wavelength data are not used in the dust model fit, the modeling tends to overpredict the intensity at the longer wavelengths. However, when the SPIRE photometry is included in the model constraints, the overprediction is greatly reduced. In the case of MIPS160 PSF (all bands used), the SPIRE250 photometry is reproduced by the model within 10%, and SPIRE500 tends to be overpredicted by less than 15%.

6.2.4. Global SED Fitting

Figure 12 shows the global SED for the galaxy NGC 6946 (similar to Figure 7 for NGC 628). The top row shows the SED for a "single-pixel" model and the second row shows the predicted model SED obtained by summing over the model SED for individual pixels. Again, we observe small differences between the global modeling (top row) and resolved studies (bottom row) due to the nonlinear behavior of the modeling process. It is seen that when the modeling is done using the MIPS160 PSF, the dust mass estimated from the multi-pixel model is in good agreement with that obtained for a single-pixel model.

Zoom In Zoom Out Reset image size

6.2.5. Fitting in Selected Apertures

Figure 13 shows SEDs for four circular apertures located on NGC 6946, sampling a wide range of surface brightnesses, $\Sigma _{L_{\rm d}}\approx (0.07\hbox{--}13)\times 10^{8}\,L_{\odot }\,{\rm kpc}^{-2}$, similar to Figure 8 for NGC 628. The top row shows the SED for a 60'' diameter circular aperture centered on the galaxy nucleus. The second row shows the SED for a 60'' aperture located on a bright spot on the spiral arms. The third and fourth rows show SEDs in 120'' apertures further from the center. Aperture 4 is located completely outside the galaxy mask, but the S/N of a large aperture allows a reliable determination of the dust model parameters.

Zoom In Zoom Out Reset image size

The nuclear emission is strongly peaked. The Herschel photometry for the central aperture noticeably decreases when the image is convolved to MIPS160 resolution, and the estimated dust luminosity surface brightness in the 60'' extraction aperture drops from 2.2 × 10⁹ L_☉ kpc⁻² (Figures 13(a) and (b)) to 1.3 × 10⁹ L_☉ kpc⁻² (Figure 13(c)).

We observe that the PACS photometry in the high surface brightness nucleus tends to generally be larger than the corresponding MIPS photometry—compare the PACS and MIPS fluxes in Figure 13(c). In apertures 2, 3, and 4, the PACS and MIPS photometry (at common resolution—see Figures 13(f), (i), and (l)), the PACS and MIPS photometry is in fairly good agreement. The disagreement in the PACS versus MIPS photometry is larger for NGC 6946 than for NGC 628 (see Table 6 below), presumably related to the much higher peak surface brightnesses present in NGC 6946. At MIPS160 resolution, $\Sigma _{L_{\rm d}}=9\times 10^7\,L_{\odot }\,{\rm kpc}^{-2}$ in NGC 628 aperture 1 is only 7% of $\Sigma _{L_{\rm d}}$ in NGC 6946 aperture 1. Table 3 summarizes the main properties of the galaxies, and the results, of the modelling.

Table 3. Galaxy Properties

	NGC 628	NGC 6946
Typea	SAc	SABcd
Da (Mpc)	7.2	6.8
Resolution at D	1'' = 34.9 pc	1'' = 33.0 pc
Galaxy optical sizeb (kpc)	11.0 × 9.9	11.4 × 9.7
Galaxy mask area A (kpc2)	274	307
Galaxy mask radius $\sqrt{A/\pi }$ (kpc)	9.4	9.9
$M({\rm H\,\mathsc{i}})$c (M☉), total in map	(3.7 ± 0.2) × 109	(5.5 ± 0.3) × 109
M(H2)d (M☉), total in map	(2.5 ± 0.6) × 109	(9.7 ± 2.7) × 109
$M({\rm H\,\mathsc{i}})$ (M☉), within mask	(1.5 ± 0.1) × 109	(2.1 ± 0.1) × 109
M(H2)d (M☉), within mask	(2.0 ± 0.2) × 109	(8.6 ± 0.5) × 109
MH (M☉), within mask	(3.5 ± 0.3) × 109	(10.7 ± 0.6) × 109
log10(O/H) + 12e	8.27 ± 0.12	8.37 ± 0.06
L(Hα)f (L☉), total	1.88 × 107	1.00 × 108
SFRg (M☉ yr−1), total	0.7 ± 0.2	4.5 ± 1.2
Fν(IRAC3.6) (Jy), within mask	0.83 ± 0.10	3.13 ± 0.52
Fν(IRAC4.5) (Jy), within mask	0.57 ± 0.06	2.22 ± 0.32
Fν(IRAC5.8) (Jy), within mask	0.98 ± 0.23	4.9 ± 1.2
Fν(IRAC8.0) (Jy), within mask	2.61 ± 0.44	13.2 ± 2.4
Fν(MIPS24) (Jy), within mask	3.04 ± 0.33	18.9 ± 2.0
Fν(MIPS70) (Jy), within mask	31.5 ± 9.8	196 ± 55
Fν(MIPS160) (Jy), within mask	104 ± 20	420 ± 127
Fν(PACSS70) (Jy), within mask	39 ± 12	245 ± 58
Fν(PACSS100) (Jy), within mask	73 ± 21	433 ± 122
Fν(PACSS160) (Jy), within mask	110 ± 20	526 ± 126
Fν(SPIRE250) (Jy), within mask	59.5 ± 7.0	247 ± 29
Fν(SPIRE350) (Jy), within mask	27.4 ± 3.3	101 ± 12
Fν(SPIRE500) (Jy), within mask	10.5 ± 1.4	35.4 ± 4.2
Md (M☉), within mask	(2.9 ± 0.4) × 107	(6.7 ± 0.6) × 107
100 × Md/MH, within mask	0.82 ± 0.17	0.63 ± 0.09
Ld (L☉), within mask	(6.8 ± 0.4) × 109	(3.3 ± 0.3) × 1010
LPDR (L☉), within mask	(7.9 ± 1.4) × 108	(4.7 ± 0.8) × 109
〈Umin〉, within mask	1.6 ± 0.3	3.1 ± 0.6
$\overline{U}$, within mask	1.7 ± 0.3	3.6 ± 0.6
〈fPDR〉, within mask	(11.6 ± 1.3)%	(14.2 ± 2.8)%
〈qPAH〉, within mask	(3.7 ± 0.3)%	(3.6 ± 0.4)%

Notes. ^aFrom Kennicutt et al. (2011). ^bMajor and minor radii at μ_B = 25 mag arcsec⁻² isophote, from Kennicutt et al. (2011). ^cWe use the "Natural" (NA) weighting maps (see Walter et al. 2008 for details). ^dFrom Leroy et al. (2009), for X_{CO, 20} = 4 (see the text). ^ePT05 H ii region abundances, unweighted average, from Moustakas et al. (2010). ^fFrom Kennicutt et al. (2008), uncorrected for internal extinction. ^gStar formation rate; from Calzetti et al. (2010), based on Hα+24 μm.

Download table as:  ASCIITypeset image

We investigate the effect of only using certain sets of cameras in our resolved (multi-pixel) modeling, since this situation will arise in galaxies that were only observed with some of the instruments. Dust estimates for NGC 628 and NGC 6946 based on different camera combinations are shown in Figure 14, where we compare the dust mass estimate with our "gold standard," the dust mass estimate using all cameras, at MIPS160 resolution.³³

Zoom In Zoom Out Reset image size

In Figures 14–17 the horizontal axis represents different combinations of PSF and cameras. The nomenclature convention is as follows. The first characters describe the PSF used; P160, S250, S350, S500, M70, and M160 refer to PACS160, SPIRE250, SPIRE350, SPIRE500, MIPS70, and MIPS160, respectively.

The first group of three digits refers to the MIPS24, MIPS70, and MIPS160 cameras, with one indicating usage. The next set of three digits similarly indicates usage of the three PACS bands, and the final three digits correspond to usage of the SPIRE cameras. For example, S350 110 111 110 uses MIPS24, MIPS70, PACS70, PACS100, PACS160, SPIRE250, and SPIRE350 cameras at SPIRE350 PSF. We always use the Scanamorphos data reduction for PACS.

The 10 right-most columns correspond to the MIPS160 PSF, and the remaining eight left-most columns have different PSFs (sorted by ascending FWHM left to right). The modeling done at MIPS160 resolution uses different camera combinations, allowing one to examine, at fixed angular resolution, the effects of different wavelength coverage and PACS/MIPS discrepancies.

Figure 14 shows that the dust mass estimates for these two galaxies appear to be quite good when we use the three SPIRE bands: the total dust mass estimated without using MIPS160 (at MIPS160 or SPIRE500 PSF) is off by only ∼10%. It therefore appears that resolved maps of dust in nearby galaxies can be reliably obtained at the resolution of SPIRE500. The dust mass estimates appear to be good provided at least two SPIRE bands (SPIRE250 and SPIRE350) are employed: the total dust mass estimated without using SPIRE500 and MIPS160 is off by 30% for NGC 628, and 22% for NGC 6946. The S250 100 111 100 dust map yields a global dust mass that exceeds the best estimate by ∼38%. For some purposes this ∼38% loss of accuracy will be acceptable, as it makes possible a dust map with FWHM = 182.

In order to map the dust at higher angular resolution (PACS160 PSF), we can only use a subset of all the cameras available (IRAC, MIPS24, PACS). Dust maps made at the resolution of PACS160 are less reliable, for three reasons: (1) with smaller pixels, the S/N of the image is poor in the faint regions, (2) there seem to be significant systematic uncertainties associated with the PACS photometry, (3) the dust models are unconstrained at λ > 160 μm, and (4) our PACS maps are less sensitive than the MIPS counterparts. The nonlinearity of the dust model fitting procedure can result in large errors in dust mass in the low S/N regions. Dust masses estimated using maps at PACS160 resolution (using IRAC, MIPS24, and PACS data only) can overestimate the dust masses by factors ≈1.4–1.8. If greater accuracy is required, PACS160 resolution should only be used in the highest surface brightness regions, where S/N issues should be less critical.

It has sometimes been argued that dust masses estimated using only data at wavelengths λ ⩽ 160 μm (e.g., those based on MIPS photometry only) are suspect because the observations are insensitive to "cold" dust at temperatures T_d ≲ 12 K. The SINGS sample included 17 galaxies with SCUBA 850 μm, and Draine et al. (2007) found that reliable dust mass estimates could be obtained with the DL07 dust model without using data longward of 160 μm: these galaxies did not appear to contain significant masses of "cold dust." The SPIRE photometry for NGC 628 and NGC 6946 confirms this: these galaxies do not appear to contain substantial masses of cold dust. In fact, we find that models fitted to the PACS photometry alone tend to predict more emission at long wavelengths than is observed. This is evident in Figures 6(a) and (d), and Figures 11(a) and (d), which show that DL07 dust models based only on λ ⩽ 160 μm data (IRAC, MIPS24, and PACS) tend to overpredict the actual fluxes at 250 and 500 μm in the central 6 arcmin (12 kpc). This is also seen in Figures 7(d) and 12(d), which compare the observed global SED with the sum of the SEDs from a multi-pixel model for each galaxy. The model (the diamonds) overpredicts the SPIRE data by factors of ∼1.1, 1.2, and 1.2 at 250, 350, and 500 μm. This is rather good agreement, and indicates that there is no substantial amount of very cold dust present in either galaxy. In fact, when the SPIRE data are used to constrain the fit, the best-fit models actually raise the starlight intensities slightly in order to match both the PACS, MIPS, and SPIRE photometry, reducing the dust mass estimate. Note the very good agreement between observed and model SEDs in Figures 7(f) and 12(f), where all data are used to constrain the multi-pixel fit.

The main difference in the estimated dust masses in the different PSF/cameras combination is due to the different values of U_min found when we do not use all the cameras available. Figure 15 shows 〈U_min〉 (defined in Equation (23)), the dust-mass-weighted value of the starlight intensity parameter U_min relative to the best estimate for 〈U_min〉—the estimate obtained using all the data at MIPS160 resolution.

Zoom In Zoom Out Reset image size

The KINGFISH sample consists of galaxies that are near enough to be well resolved, enabling us to make maps of dust and starlight properties. However, because the dust modeling used here is a nonlinear procedure, the total dust mass estimate, for example, will depend on what resolution is employed. In the preceding sections, we usually made dust maps using the best angular resolution possible with each chosen camera set used, with the limiting resolution therefore determined by the longest wavelength data used in the fit (or the inclusion of MIPS160).

Herschel is frequently used to observe distant galaxies that are unresolved at the longest wavelengths. For the KINGFISH galaxies, we will therefore also estimate dust masses using the global photometry (the flux from within the galaxy mask), to see how dust mass estimates based on global photometry compare with the best estimates for a well-resolved galaxy.

Figure 16 shows the dust mass estimated using the DL07 dust models to fit the global SED, divided by the "gold standard," i.e., the sum of the dust mass obtained from pixel-by-pixel modeling of the resolved galaxy in the MIPS160 PSF using all the cameras. Recall that the DL07 model assumes that most of the dust is exposed to a single starlight intensity U_min. A multi-pixel model with U_{min, j} varying from pixel to pixel cannot be exactly reproduced by a model that assumes a single value of U_min for the entire galaxy, hence fitting the summed photometry with a single-pixel model will in general give a dust mass that will differ from that obtained by summing over a multi-pixel model.

Zoom In Zoom Out Reset image size

Figure 16 shows that if all (IRAC, MIPS, PACS, SPIRE) data are used, the dust mass obtained by fitting the global SED is very close to that obtained from a multi-pixel fit of the resolved galaxy: the global dust mass is within 18% of the best estimate for NGC 628, and within 3% of the best estimate for NGC 6946. If data from some cameras are not used, the single-pixel estimate for the dust mass will of course change, but we see that for NGC 628 and NGC 6946 the resulting dust masses vary surprisingly little provided we have SPIRE data or MIPS160. If we use only PACS for λ ⩾ 70 μm, the dust mass estimate for NGC 628 is high by ∼10%. The M160$\_$111$\_$000$\_$000 column in Figure 16 compares the global dust mass estimated from a single-pixel fit using only Spitzer (IRAC and MIPS) data versus our "gold standard" best estimate (multi-pixel model using all cameras). The dust masses estimated using only the global photometry from Spitzer are within ∼30% of our present best estimates. Thus, at least for normal-metallicity star-forming galaxies such as NGC 628 and NGC 6946, dust mass estimates based on IRAC and MIPS photometry, e.g., the dust masses estimated for the SINGS sample (Draine et al. 2007) appear to be reliable.

Figure 17 compares the global dust mass estimated using global fluxes to the dust masses estimated summing the dust masses over each map, for different sets of cameras. Essentially, it is a measure of the nonlinear behavior of the dust modeling. For each camera and PSF combination, performing resolved modeling allows the dust to be colder in some regions, and this introduces more dust mass. Modeling based on global photometry underpredicts the dust masses by 0%–20% in most cases.

Zoom In Zoom Out Reset image size

In principle, using the MIPS160 PSF results as our "gold standard" best estimate for the dust parameters may also suffer from the above-mentioned nonlinear behavior: each MIPS160 pixel covers a rather large area of the galaxy, in which the dust properties (e.g., q_PAH) or radiation field (e.g., U_min) may have variations. Unfortunately, the discrepancy of PACS160 and MIPS160 fluxes makes it important to retain MIPS160 photometry, and, therefore, to use the broader MIPS160 PSF.

Table 4 summarizes the main parameters obtained by single-pixel and multi-pixel modeling of the galaxies at MIPS160 resolution, using all the cameras.

Table 4. Multi-pixel and Single-pixel Best-fit Model Parameters at MIPS160 Resolution

	NGC 628		NGC 6946
Parameter	Resolveda	Globalb	Resolveda	Globalb
Md (M☉)	2.87  ×  107	2.33  ×  107	6.74  ×  107	6.62 × 107
Ld (L☉)	6.83  ×  109	7.08  ×  109	3.31  ×  1010	3.20 × 1010
Umin	1.46	2.00	3.05	3.00
$\overline{U}$	1.75	2.23	3.61	3.55
100 × fPDR	11.6	10.2	14.2	15.3
100 × qPAH	3.67	3.70	3.55	3.30
100 × Md/MHc	0.82	0.66	0.63	0.61

Notes. ^aTotal or mean values are defined in Equations (21)–(24). ^bValues are from a model fit to the global photometry. ^cFor X_{CO, 20} = 4.

Download table as:  ASCIITypeset image

Recently Galliano et al. (2011), in a resolved study of the dust in the LMC, advocated a starlight heating intensity distribution function different from what is employed in the current work. They recommended using the starlight distribution function of Dale & Helou (2002), with a power-law distribution between two adjustable parameters U_min and U_max, with adjustable power-law exponent α. We note that Galliano et al. (2011) use a slightly different nomenclature: U_max ≡ U_min + ΔU. The distribution used in the present work (Equations (9) and (10)) incorporates an additional "delta function" component at U_min, but fixes the value of U_max = 10⁷.

Galliano et al. (2011) claimed that satisfactory fits can be obtained without the "delta function" component. We note that while we fix U_max to a given value in our model, Galliano et al. (2011) use U_max as a free parameter. Both models have thus the same number of degrees of freedom. Galliano et al. (2011) based their claim on maps of the LMC using IRAC, MIPS, and SPIRE. Here we test this claim using our maps of NGC 628 and NGC 6946.

Figure 18 shows maps of χ² (goodness of fit) obtained for three different starlight distribution schemes, all done at MIPS160 resolution, using all 13 cameras. The top row images correspond to NGC 0628 and the bottom row images correspond to NGC 6946. In the left column the starlight distribution is a power law between U_min and U_max with power-law exponent α (all adjusted) without the delta function contribution at U_min, as proposed by Galliano et al. (2011). It has six adjustable parameters (Ω_⋆, M_dust, q_PAH, U_min, U_max, α), and therefore 13 − 6 = 7 degrees of freedom per pixel. The fit for this starlight distribution is poor, giving χ² in the range 10–20 in the bright regions of the galaxies. In the middle column the starlight distribution is a power law between U_min and fixed U_max = 10⁷, with power-law exponent α (only U_min and α are adjusted), plus a delta function contribution at U_min. This fit produces excellent results and is the one recommended in the present work. We note that this fit has the same number of degrees of freedom as the ones in the left column panels; the adjustable parameter U_max is replaced by the normalization of the "delta function" via the parameter γ. In the right column the starlight distribution is a power law between U_min and U_max with power-law exponent α (all adjusted), plus a delta function contribution at U_min. Although the right panels have one extra adjustable parameter (U_max) with respect to the recommended fit, the fit quality is similar to the middle panels. If U_max were a relevant parameter, one would expect χ² to decrease significantly (of order 1). In fact, in the right column fit, for most of the bright pixels the best-fit U_max value is U_max = 10⁷ and thus the value of χ² is for those pixels exactly the same as in the center column fit.

Zoom In Zoom Out Reset image size

We found that the χ² values of the distribution used in the current work are significantly smaller than the ones found using the Galliano et al. (2011) distribution. Therefore, we recommend adding the "delta function" component at U_min. Since the fit does not improve significantly by allowing U_max to be fitted, we recommend fixing it to U_max = 10⁷, leading to a simpler, more robust fit. It should be noted that the Galliano et al. (2011) study is based on the dust in the LMC, and in our present work we use NGC 628 and NGC 6946, two galaxies with metallicities closer to the MW metallicity. Additionally, due to the proximity of the LMC, each pixel covers a physical area much smaller than our modeling pixels.

10.1. Dust Mass Estimation

Dust mass estimates are, of course, model dependent. Above we have estimated the dust mass using a specific dust model, the DL07 silicate-graphite-PAH model, with an assumed parametric form for the starlight intensities heating the dust. The physical dust model and the ansatz for the distribution of starlight intensities together appear to successfully reproduce the observed SEDs, and to give reliable estimates of the dust masses.

As discussed in Section 4, the far-infrared opacity of the amorphous silicate originally put forward by DL84 was subsequently adjusted slightly (Li & Draine 2001) to improve agreement with the far-infrared and submillimeter emission observed for the local high-latitude dust. Thus the DL07 model uses dust properties that can reproduce the observed FIR-submillimeter emission from the MW cirrus with a single starlight heating intensity—no "cold dust" is needed for the MW cirrus. In this dust model, the opacities are assumed to be independent of the dust temperatures.

Here we see that the same dust model, with suitable adjustment of the starlight assumed to be heating the dust, is able to reproduce the emission from NGC 628 and NGC 6946 out to 500 μm, without introducing any "cold dust" component. The modeling performed at MIPS160 resolution, using all the IRAC, MIPS, PACS, and SPIRE cameras, in addition to being able to reproduce the observed SED, gives dust masses that are in line with the expected dust/H mass ratios for these galaxies: the dust/H mass ratio images in Figures 5 and 10, at MIPS160 resolution, are smooth, and the global dust/H mass ratios are 0.0082 ± 0.0017 and 0.0063 ± 0.0009 for NGC 628 and NGC 6946, respectively.

Observed depletions in the local MW indicate a dust/H mass ratio M_dust/M_H = 0.0091 ± 0.0006 (Draine 2011b, Table 23.1). Thus, a galaxy with a similar fraction of interstellar heavy elements in dust would be expected to have

$$\begin{equation} \frac{M_{\rm dust}}{M_{{\rm H}}} \approx 0.0091 \times 10^{[{\rm O}]} , \end{equation} \tag{ 29 }$$

where [O] ≡ log₁₀[(O/H)/(O/H)_☉]. Thus galaxies with heavy element abundances (and depletions) similar to the local MW should have M_dust/M_H ≈ 0.009.

NGC 628 and NGC 6946 appear to be mature star-forming galaxies that would be expected to have interstellar metallicities similar to the local MW. If so, then the values of M_dust/M_H found by fitting a dust model to the observations appear to be in excellent agreement with expectations. We note that the H ii region oxygen abundances found by Moustakas et al. (2010) together with (A_O)_☉ = 12 + log₁₀(O/H)_☉ = 8.73 (Asplund et al. 2009) yield [O] = −0.46 and −0.36 for NGC 628 and NGC 6946, and Equation (29) would predict M_d/M_H = 0.0032 and 0.0040 for NGC 628 and NGC 6946, respectively. However, we suspect that the PT05 oxygen abundance estimates may be biased low: Moustakas et al. (2010) list PT05 H ii region oxygen abundances for 38 galaxies; the highest oxygen abundance is A_O = 8.59 ± 0.11 for NGC 4826. It seems unlikely that none of the galaxies in their sample have oxygen abundances that are solar or supersolar.

10.2. Dust Opacity in Molecular Gas and the Value of X_CO

10.3. Dust Mass Estimates from Single-temperature Fits

10.4. Radial Gradients in Dust and Starlight

10.5. Interpretation of the Starlight Heating Parameters α and γ

As discussed in Section 4.4, a fraction γ of the dust mass is taken to be heated by a power-law distribution of starlight intensities between U_min and U_max with dM/dU∝U^−α, and the remaining fraction (1 − γ) of the dust mass is heated by a radiation field with intensity U_min. Dust heated by this simple parameterization of the starlight intensities is quite successful in reproducing the observed SED (see the SEDs in Figures 7(c) and (f); 8(c), (f), (i), and (l); 12(c) and (f); and 13(c), (f), (i), and (l)).

Dale et al. (2001) discussed two ideal cases that have analytical predictions for the value of α: a single star in a homogeneous diffuse medium, and a dark cloud. The first case, with U∝r⁻², has dM_d/dU∝U^−2.5, i.e., α = 2.5. In this case, U_max would be very large, ≳ 10¹⁰, the heating rate at which dust grains would vaporize. In the second case, a slab where the heating intensity is primarily attenuated by dust absorption, one would get dM_d/dU∝U⁻¹, i.e., α = 1, and U_max would be set by the starlight intensity at the edge of the slab. This case might correspond to a PDR, the edge of a dark cloud. Different clouds in the pixel would have different values of U_max, so their superposition could be approximated as a power law with a slightly larger value of α. We expect dust in a real galaxy to have 1 ≲ α ≲ 2.5.

Consider for the moment the fraction γ of dust with dM/dU∝U^−α. The power radiated is dL ∝ U dM ∝ U^{2 − α} dlog U. A value of α = 2 would have uniform dust luminosity per unit interval in log U. A value α < 2 would concentrate the dust luminosity L in the high radiation field regions (U ≈ U_max), and a value α > 2 would have L dominated by dust with U ≈ U_min.

Figure 19 shows the starlight power-law component parameters α and γ. The power-law index α has a preferred value α ≈ 1.6 in the bright regions. In these regions, the power-law distribution is representing the heating of dust in high-U regions, e.g., PDRs near OB stars. There are also regions with α ≈ 2.4. In these regions, most of the dust luminosity is concentrated in regions with U ≈ U_min, and the power-law component is essentially making the U = U_min component broader, i.e., is allowing for variations in the starlight intensities in the diffuse ISM.

Zoom In Zoom Out Reset image size

We do not attach great physical significance to the parameters γ and α. The more physically meaningful quantities are the mass-weighted mean starlight intensity $\overline{U}$, and f_PDR, the fraction of the dust luminosity that originates in regions with U > 10². In Figures 5 and 10, one sees that f_PDR > 0.05 over most of the galaxy mask for both galaxies.

Using Spitzer and Herschel data, we perform a resolved study of the dust physical parameters, and of the starlight heating the dust, in the nearby galaxies NGC 628 and NGC 6946. We employ the DL07 dust model, having amorphous silicate grains and carbonaceous grains, including PAHs (see Section 4). The model has a distribution of grain sizes, and we allow for a distribution of starlight intensities heating the dust within each pixel. The model assumes a frequency-dependent opacity that reproduces the observed emission spectrum of high-latitude dust.

In order to perform resolved studies of the galaxies, it is important to convolve all the images into a common PSF. This is achieved using the convolution kernels generated by Aniano et al. (2011). We perform our dust modeling using all appropriate combinations of PSFs and cameras. Table 1 lists some of the resolutions and cameras used in this work.

Our principal findings are as follows.

• 1.  
  The dust model is quite successful in reproducing the observed SED over the full wavelength range 6–500 μm where dust emission dominates (see the SEDs in Figures 7(c) and (f); 8(c), (f), (i), and (l); 12(c) and (f); and 13(c), (f), (i), and (l); and the model/observed maps in Figures 6(c) and (f); 11(c) and (f)). The dust model reproduces the emission out to 500 μm without introduction of a "cold dust" component, and with no allowance for possible temperature-dependent dust opacities. The DL07 dust opacities therefore appear to be consistent with both the observed emission from local "cirrus" at high galactic latitudes, and the observed SED from 500 pc sized regions of the ISM in NGC 628 and 6946. There is no indication that T-dependent opacities are needed to reproduce these data.
• 2.  
  Maps of the dust/H mass ratio show it to be relatively uniform over both galaxies, outside of the nucleus (see Figures 4(l) and 9(l)). The derived dust/H mass ratio for these galaxies is consistent with that expected if the interstellar abundances in NGC 628 and NGC 6946 are close to solar. As near-solar abundances seem likely, this is strong support for the quantitative accuracy of the dust masses obtained by modeling the IR and FIR emission.
• 3.  
  Figure 14 shows how dust mass estimates depend on the camera and PSF used. The "gold standard" is taken to be dust modeling done at MIPS160 resolution, using all cameras, with Scanamorphos (Roussel 2012) processing of the PACS data. Relative to this standard, mass estimates done with smaller pixels (giving up the lowest resolution cameras) tend to be biased slightly high. At SPIRE250 resolution, the bias is ∼38%. Working at SPIRE350 resolution provides a good compromise between resolution and accuracy, with ∼30% errors for the total dust mass estimated for NGC 628 and NGC 6946.
• 4.  
  Even though the dust modeling process is nonlinear, resolved modeling is compatible with modeling using the global photometry of the galaxy. Differences in the inferred parameters are in most cases under 20%.
• 5.  
  In NGC 6946 the dust/H mass ratio calculated with a single value of X_{CO, 20} = 4 appears to have a minimum at the center (see Figures 10(a)–(c)). This minimum is interpreted as an artifact due to overestimation of the molecular mass. The dust surface densities found here therefore support the finding by Donovan Meyer et al. (2012) of X_{CO, 20} ≈ 1.2 in the center of NGC 6946.
• 6.  
  The present fitting procedures employ six adjustable parameters for each pixel: Ω_⋆, M_dust, q_PAH, U_min, α, and γ (we fix U_max = 10⁷). For NGC 628 and NGC 6946 we find that this approach provides substantially better fits than the approach advocated by Galliano et al. (2011) who treat Ω_⋆, M_dust, q_PAH, U_min, α, and U_max, as adjustable parameters and fix γ = 1: compare the χ² maps in Figure 18(a) with 18(b), and 18(d) with 18(e).
• 7.  
  The starlight heating the dust within a single pixel is characterized by three adjustable parameters: U_min, γ, and α (see Equations (9) and (10)). U_min is interpreted as the intensity of the diffuse starlight, responsible for the bulk of the dust heating. Maps of U_min (see Figures 2 and 7) show significant structure: there is a tendency of U_min to be higher in spiral arms, as well as a significant radial decline in U_min. The overall mass-weighted mean value of U_min = 1.5 for NGC 628 and U_min = 3.1 for NGC 6946, but U_min declines to values of 0.3 or lower in the outer regions of both galaxies.
• 8.  
  The parameter f_PDR is the fraction of the dust heating power that is contributed by starlight intensities U > 100. The overall value is 〈f_PDR〉 = 0.12 for NGC 628 and 〈f_PDR〉 = 0.14 for NGC 6946, but in both galaxies f_PDR peaks at local maxima of the dust surface brightness, which are associated with star-forming regions. At SPIRE250 resolution, f_PDR reaches values as high as 0.25 in both galaxies.
• 9.  
  Current PACS photometry disagrees with MIPS photometry (see Figures 27 and 28). PACS70 photometry is up to 80% higher than MIPS70, and PACS160 is up to 50% higher than MIPS70 photometry in the bright areas of the galaxies.
• 10.  
  The PACS–MIPS photometry disagreement induces a bias when trying to model dust at high spatial resolutions (PACS160 PSF), where MIPS70 and MIPS160 cannot be used. We do not recommend modeling dust at PACS160 resolution in low surface brightness areas. Modeling done at SPIRE250 resolution is more reliable than modeling at PACS160, but the inferred dust masses still disagree with our "gold standard" (i.e., modeling done at MIPS160 resolution using all the IRAC, MIPS, PACS, and SPIRE cameras) by up to ≈30%.

Subsequent work (G. Aniano et al. 2012, in preparation) will extend this study to all 61 galaxies in the full KINGFISH sample (Kennicutt et al. 2011).

We are grateful to R. H. Lupton for availability of the SM graphics program, and to the anonymous referee for helpful suggestions.

This research was supported in part by JPL grants 1329088 and 1373687, and by NSF grant AST1008570.

Given an image I expressed in its native pixel grid (x, y), and a (potentially empty) set NB_Std of "non-background" sky regions expressed in a standard 4'' × 4'' grid (a grid common to all the images), we iteratively construct a sequence of "background" pixels masks B_n, n = 0, 1, 2, ..., and a sequence of tilted planes plane_n, n = 0, 1, 2, hellip; that approximate the image values over the sets B_n, n = 0, 1, 2, hellip;. The sets B_n, n = 0, 1, 2, hellip; will exclude the "non-background" NB_Std regions, i.e., it is assumed that the pixels in NB_Std are not background-candidate pixels.

Throughout the algorithm, the given "non-background" regions will be masked and not considered in the fitting. We start by bi-linearly interpolating or averaging (depending on the grid relative pixel sizes) the set NB_Std into the image I grid (x, y) to obtain NB_Prev.

We first smooth the original image by convolving it with a (normalized) Gaussian kernel ∝exp [ − (r/r₀)²], with r₀ equal to 0.56 times the width of the native pixel for each camera. The factor 0.56 is chosen so that each pixel contributes half of its flux to the corresponding pixel in the convolved image. This step provides a less noisy image, in which we can more robustly identify the galaxy and the external sources.

We proceed as follows. We start by setting B₀ to be all pixels in the grid (x, y), except for the pixels in NB_prev (if any). We fit a plane plane₀ to the image values over the B₀ pixels. This is achieved by finding the plane₀ that minimizes ₀ defined as

$$\begin{equation} \epsilon _n \equiv \sum _{(x,y) \in B_n} {[I(x,y)-{\rm plane}_n(x,y)]^2}, \end{equation} \tag{ B1 }$$

where n = 0. We compute the dispersion θ₀ defined as

$$\begin{equation} \theta _n \equiv \sqrt{\frac{\epsilon _n}{\#\lbrace (x,y)\in B_n\rbrace -3}}, \end{equation} \tag{ B2 }$$

where n = 0, so the sum is over the B₀ (original "background" candidate) pixels.

We take the first (n = 1) set of "background pixels"³⁵ B₁ to be all pixels that do not belong to NB_prev, or lie more than κ θ₀ above plane₀ (κ is a constant, and we use κ ≡ 1.8, see below for details):

$$\begin{equation} (x,y)\in B_1 \iff \left\lbrace \begin{array}{l}(I(x,y)- {\rm plane}_0)< \kappa \, \theta _o \\ \ms (x,y) \notin {\rm NB}_{{\rm prev}} \end{array}. \right. \end{equation} \tag{ B3 }$$

In the first step, these pixels typically correspond to the brightest point sources and central pixels of the galaxy.

We now proceed iteratively by

• 1.  
  Fitting a tilted plane plane_n to all the candidate "background pixels" B_n, by minimizing _n defined by Equation (B1).
• 2.  
  Computing θ_n, defined by Equation (B2), where the sum extends over the B_n pixels.
• 3.  
  Identifying a new set of "background pixels" B_{n + 1} to be those pixels which lie inside the interval (plane_n − κ θ_n, plane_n + κ θ_n), and do not belong to NB_prev:

  $$\begin{equation} (x,y)\in B_n \iff \left\lbrace \begin{array}{l}|I(x,y)- {\rm plane}_n |< \kappa \, \theta _o \\ \ms (x,y) \notin {\rm NB}_{{\rm prev}} \end{array}. \right. \end{equation} \tag{ B4 }$$

This iteration will generate a sequence of tilted planes plane_n, and a (generally decreasing) set of "background pixels" B_n. The iteration will converge, i.e., it reaches a step M where B_{M + 1} = B_M, typically in M = 20–250 iterations.

Once the iteration converges, we have set of "background pixels" B_M, and a fitting plane plane_M (all the computed sets B_n and planes plane_n will depend on the choice of κ, see below for a discussion). We bi-linearly interpolate or average B_M to the standard 4'' × 4'' grid. This will generate a mask of the "background regions" in the standard grid, so it can be compared with the results of the algorithm from images of different cameras. We can generate standard mask of "non-background regions" by taking the complement of the "background regions" mask. The left columns of Figures 22 and 23 show the resulting residuals for NGC 628 and NGC 6946, respectively, over B_M, for κ = 1.8. Rows 1–4 correspond to IRAC4.5, MIPS24, PACS160, and SPIRE350, respectively. It can be seen that the algorithm masks out all the point sources, extended sources, and several background regions. Residuals slightly positive are found in the edges of the extended emission areas, although most of the extended emission areas are recognized and excluded.

After the iteration converges to a set of "background pixels" B_M, if we desire to perform a background subtraction to the images, we run the iterative steps again. In the second iteration, we use the original (non-smoothed) image, and the complement of the set B_M found as the starting set of "non-background pixels" NB_prev, i.e., the pixels which were recognized as having emission in the first iteration are masked out. The second run of the iteration will converge to a new set of "background pixels" B_N, and a new fitting plane plane_N, typically in N = 20–250 iterations. Clearly, by construction, we always have B_N⊆B_M. Finally, the last plane plane_N is subtracted from I.

The reason to perform the second iteration, with the original image and a larger set of starting "non-background pixels," is that we can remove noisy pixels from the background pixels, recognize the edges of the emission areas more precisely, and perform a more precise background level estimation. Moreover, if we try to perform a single iteration (i.e., by going directly to the second iteration), then the larger image noise will prevent us from recognizing the dimmer sources and extended emission areas. The center columns of Figures 22 and 23 show the resulting residuals for NGC 628 and NGC 6946, respectively, over the sets B_N.

It can be shown that if there exist an infinite set of background points B₀ in which the intensities are independent noise drawn from a distribution with probability density ρ, then the algorithm will asymptotically converge: $B_n \mathop{\longrightarrow}\limits_{n\rightarrow \infty }{} B_{N} \subseteq B_0$, that will depend on κ. We define f(κ), the fraction of the original background-candidate pixels that are asymptotically considered background:

$$\begin{equation} f(\kappa) \equiv \#\lbrace B_{N}\rbrace /\#\lbrace B_0\rbrace. \end{equation} \tag{ B5 }$$

For noise distribution function densities ρ which are analytical around their expectation values, it can be shown that if $\kappa < \sqrt{3}$ then B_n will converge to an empty set ∅; and B_n will converge to a non-empty set provided that $\kappa > \sqrt{3}$, i.e., $f(\kappa < \sqrt{3})=0,\,\,\,\,f(\kappa > \sqrt{3})>0$. In particular, for a top hat distribution (ρ(x) = 1/2⇔|x| ⩽ 1), B_N = B₀ if $\kappa > \sqrt{3}$, i.e., $f(\kappa > \sqrt{3})=1$. If the noise has a Gaussian distribution with standard deviation σ, we can numerically compute f(κ), and the actual σ/θ (the ratio between the distribution standard deviation σ and the computed final dispersion θ). Figure 21 shows the values of f(κ) (dotted green line) and θ/σ (solid blue line) for 1.6 < κ < 2.6 for a Gaussian distribution. The vertical red lines correspond to the threshold value $\kappa = \sqrt{3}$, and the used value κ = 1.8.

Zoom In Zoom Out Reset image size

Our choice of κ is dictated by (1) the need to have $\kappa > \sqrt{3}$ in order for the procedure to converge to a non-empty set, (2) the desire to minimize inclusion of outliers that may be due to rare events (e.g., cosmic rays, image artifacts), and (3) the desire to be as sensitive as possible to dim sources. We choose κ = 1.8 as a good compromise. If the background noise were Gaussian, for our choice κ = 1.8, we have f(κ = 1.8) = 0.5506, θ/σ = 0.4202 (i.e., the final iteration should fit only ≈55% of the background points, and their measured dispersion should be 0.4203σ).

How well does the iterative procedure behave with real images? Typically the iterative procedure converges in 20–250 steps. The right columns of Figures 22 and 23 show the resulting histograms of the background-subtracted intensities in B_N for NGC 628 and NGC 6946, respectively. The solid red curve is a Gaussian distribution function with σ = θ/0.4202, i.e., with the standard deviation expected from the measured dispersion. The fit of the Gaussian curve to the histograms in Figures 22 and 23 is very good, supporting the validity of the presented algorithm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This background recognition algorithm has several advantages over other approaches. First, the algorithm is very robust, automatic, and quick, and produces superior results when dealing with all the ≈1400 images of the KINGFISH sample. Second, by keeping only a subset (≈ 55%) of the background pixels, the last computed dispersion is small enough that even dim sources are easily recognized and masked. This allows us to exclude virtually all detectable foreground or background sources from the subset of points used to determine the background of the images. Third, this algorithm also recognizes areas with dim extended emission. Fourth, image artifacts over the background areas are easily identified, allowing us to minimize bias induced by them in the background estimation. This allows us to avoid over- or under-subtraction of the background in each image, which is critical in order to model the galaxy low surface brightness areas. Table 5 summarizes the background and image information for the different cameras.

Some of the original images include bad pixels due to detector saturation, incomplete scanning, or other problems. This effect is not important when present in the pixels far away from the target galaxy, but may have undesired consequences if not treated appropriately when present near or within the target galaxy. A similar situation arises when observations are made over a small part of a galaxy, for example, with the IRS instrument on board the Spitzer Space Telescope. We develop a method for correcting this problem.

In order to estimate (and correct) the influence of bad pixels of a camera A, we start by choosing a second camera B. We require that the sky coverage of camera B be larger than camera A, or at least that does not present bad pixels in the same areas as A. From all the possible candidates for being the correcting camera B, we choose the camera with wavelength most similar to A. For example, we typically correct the MIPS160 camera with the SPIRE250 camera.³⁶

We divide the original image grid of the camera A, OG_A, into two sets of pixels: OG_{A, good}, the pixels of OG_A with good A data, and OG_{A, bad}, the pixels of OG_A in which the camera A has bad data.

We interpolate the background-subtracted image B from its original grid OG_B into the OG_A grid, obtaining B_interp A, which looks morphologically similar to A. We construct two images B_good and B_bad, both defined on OG_A, as follow:

$$\begin{equation} B_{{\rm good}}(x,y) \in {\rm OG}_A = \left\lbrace \begin{array}{lcc}B_{{\rm interp}\,A} &\iff &(x,y) \in {\rm OG}_{A,{\rm good}} \\ \ms 0 &\iff &(x,y) \in {\rm OG}_{A,{\rm bad}} \end{array} \right. \end{equation} \tag{ C1 }$$

$$\begin{equation} B_{{\rm bad}}(x,y) \in {\rm OG}_A = \left\lbrace \begin{array}{lcc}0 &\iff &(x,y) \in {\rm OG}_{A,{\rm good}} \\ \ms B_{{\rm interp}\,A} &\iff &(x,y) \in {\rm OG}_{A,{\rm bad}} \end{array}. \right. \end{equation} \tag{ C2 }$$

We process the images B_good and B_bad exactly as we do with the image A; we rotate them, convolve with the kernel that corresponds to the image A, and interpolate them to the final grid (FG). A given output pixel in FG will have flux flux_goodB, and flux_badB in the final stage of processing the B_good and B_bad images, respectively. We define F_{A, B}(x, y) ∈ FG, the fraction of the total flux in the FG coming from areas with good A data, as

$$\begin{equation} F_{A,B}={{\rm flux}_{{\rm good}}B\over {\rm flux}_{{\rm good}}B +{\rm flux}_{{\rm bad}}B}. \end{equation} \tag{ C3 }$$

F_{A, B} provides an estimate of how much flux we are missing in each output pixel due to the missing data in the original image A. In the pixels for which F_{A, B} > 0.65 (i.e., those where most of the flux is actually coming from regions where we know A) we can further compensate the final stage of the A image processing by multiplying the value obtained by 1/F. In the points F_{A, B} < 0.65 a significant amount of the flux (more than 35%) is coming from regions where A is not known, so we mask and ignore them.

This correction algorithm would be exact if the color A/B were constant, and the PSFs of both cameras coincided. We test the algorithm by purposely removing some part of an image, and correcting it with itself or another image. In the former case, the results are exact within numerical noise. In the latter case, because we are correcting an image A with a different one B, the method is not exact, but in most cases still provides a quite accurate result.

The χ² maps also shed light on the estimated flux uncertainties. The value of

$$\begin{equation} \chi ^2 \equiv \sum _{k} \frac{[F(\lambda _k)-F({\rm model},\lambda _k)]^2}{\sigma _{\lambda _k}^2} \end{equation} \tag{ D1 }$$

depends not only on the ability of the model to fit the data, but also on the adopted uncertainties $\sigma _{\lambda _k}$. For a given data set and fitted model, underestimation of the uncertainties $\sigma _{\lambda _k}^2$ would lead to higher vales of χ². Recall that the $\sigma _{\lambda _k}$ were estimated by studying the pixel-to-pixel variations in the background regions outside the galaxy mask. If part of the background variations actually comes from real variations in emission from foreground dust, or dusty background galaxies, then our dust model procedure may actually be able to partially fit the background variations, leading to a low value of χ². On the other hand, image artifacts and departures of the real data from the models will lead to higher values of χ². It is clear, for example, that the differences in the PACS and MIPS photometry must lead to an increase in the χ² values, since the model cannot fit simultaneously two different flux values.

If we have N measures from an ideal model with M adjustable parameters, each measure having independent Gaussian noise with variance $\sigma _{\lambda _k}^2$, then the χ² of the fit will have a chi-square distribution with degrees of freedom (dof), where dof = N − M. The χ²/dof map will have mean 1, and dispersion $\sqrt{2/{\rm dof}}$.

In our background noise estimation procedure, after we recognize and mask all the background sources and subtract a "tilted plane," in each image there remains a zero-mean background due to unresolved dim sources. We use the dispersion of the background pixels (with respect to the best-fit plane) to estimate the stochastic background uncertainty. This provides a valid uncertainty estimate since the same background structures will be present over the galaxy, so our galaxy fluxes will be uncertain (at least) at this level. However, this background dispersion is not only due to random noise in the detectors, it also includes a contribution from foreground and background sources. Figure 24 shows the correlation of the dispersion of the background pixels from the MIPS160 (horizontal axis) and SPIRE350 (vertical axis) cameras. Although the residual distributions have zero mean, there is a strong correlation between the MIPS160 and SPIRE350 flux in the background pixels, demonstrating a component coming from real astrophysical sources. Other cameras show similar correlations. The modeling can fit these correlated departures from zero by adding positive or negative amounts of real dust, reproducing the departures and therefore decreasing the χ² of the fit. We therefore expect to have χ²/dof < 1 in areas where the flux uncertainties have an important contribution from unresolved dim sources (i.e., in the background pixels or low surface brightness areas).

Zoom In Zoom Out Reset image size

Figure 25 shows the contribution of different cameras to χ² for NGC 628. The fit is done at MIPS160 resolution, using all 13 cameras available. We employ the recommended dust model with six adjusted parameters (Ω_⋆, M_dust, q_PAH, γ, U_min, α) leading to a fit with dof = 7. The upper left panel (a) shows the χ²/dof map. If the data were artifact-free, perfectly modeled, with independent Gaussian noise, the map would have mean 1, and dispersion $\sqrt{2/7}\approx 0.53$. The actual χ²/dof values are slightly above 1 in the bright regions and below 1 in the outer parts of the galaxy. In the outer parts of the galaxies, the correlation of the noise between the cameras drives the χ²/dof below the expected value 1. In the bright spots of the galaxy, the differences in PACS and MIPS photometry decrease the quality of the fit. We estimate the PACS and MIPS photometry uncertainties in a way that their contribution to the χ² maps is controlled and consistent, but the model will try to fit some intermediate (probably non-physical) value which may impair its ability of accurately fitting the remaining cameras.

Zoom In Zoom Out Reset image size

In order to remedy this situation, uncertainties for the PACS and MIPS photometry are estimated by comparing the PACS and MIPS images (see below). The uncertainties are slightly overestimated, thus giving more weight to the remaining artifact-free cameras. Figure 25(b) shows the PACS(S)70/MIPS70 ratio map. Stripes aligned with the MIPS sky scanning directions near the galaxy center suggest the presence of MIPS image artifacts. Discrepancies in the low surface brightness areas (outer parts of the galaxy) suggest differences in the ability of the telescopes to measure the low surface brightness areas correctly. The top right panel (c) shows the assigned MIPS70 S/N map. The procedure described in Appendix D.1 assigns large uncertainties to the areas where the PACS and MIPS images have discrepancies, leading to reduced S/N. This is a critical step toward getting χ² maps that really measure the goodness of the model fit, not only the data discrepancies.

The lower row of Figure 25 shows the contribution Δχ² of selected cameras to the total χ² for NGC 628. The lower left panel (d) shows the contribution of MIPS24 camera. We note that this contribution is extremely low, under 0.05 for most of the galaxy. This is because the dust model can adjust the 24 μm flux relatively independently of the other observed bands by adjusting γ and α. The 24 μm band is unique in lying a factor of three in wavelength away from the nearest other bands (8 μm and 70 μm). Figure 25(e) is the contribution of MIPS70 μm to the total χ². Even though the MIPS70 camera has severe artifacts, its Δχ² near the center is controlled by the large values of σ assigned. Figure 25(f) shows Δχ² for SPIRE350. The camera behaves exactly as expected, contributing Δχ² ≈ 0.6 near the galaxy nucleus and smaller Δχ² in the low surface brightness areas. The behavior for NGC 6946 is similar.

D.1. Uncertainty Estimate Procedure

In order to estimate uncertainties for the different cameras, we proceed as follows.

For each camera of wavelength λ_k, after the image processing (rotation to R.A.–decl., background subtraction, convolution to a common PSF, correction for boundary effects and missing data, and resampling to the FG) the flux in each final pixel is a (known) linear combination of the flux of (in principle) all the original pixels of the camera. If the statistical properties of the uncertainties in the original pixel fluxes were known, it would be possible to infer the uncertainties (and their statistical properties) in each final pixel. The original maps oversample the beam, and have artifacts that extend over several pixels, so realistic statistical properties of the uncertainties are difficult to determine. We therefore estimate the uncertainties directly in the final (post-processed) image.

We always assume that the pixel noise is not correlated between different cameras, so in what follows we work in each image independently, and omit the subindex identifying the camera in the following discussion. For each pixel, we will consider two independent sources of uncertainties δ_sto and δ_sys. The term δ_sto will be the uncertainty due to stochastic pixel noise and the existence of (unresolved) background dim sources. We estimate it by measuring the flux fluctuations in the background pixels. The term δ_sys will include systematic and calibration uncertainties.

In order to estimate the term δ_sto, the stochastic pixel noise, we proceed as follows. We re-grid the BM to the final-map grid to define a set of N_bg "background pixels." For each camera k, we compute the background dispersion $\delta _{{\rm sto},\lambda _k}$ as

$$\begin{equation} \delta _{{\rm sto},\lambda _k} = \sqrt{\frac{1}{(N_{\rm bg}-3)}\sum _{{\rm (x,y)\in {\rm (BM)}}}{\left[I_{\lambda _k}^{{\rm obs}}(x,y)\right]^2}}, \end{equation} \tag{ D2 }$$

where $I_{\lambda _k}^{{\rm obs}}$ is the background-subtracted flux of the camera of wavelength λ_k in the pixel (x, y) in the last stage of the image processing (i.e., after the image has been rotated to R.A.–decl., background plane subtracted, convolved to a common PSF, corrected for boundary effects and missing data, and resampled to the FG). The background dispersion $\delta _{{\rm sto},\lambda _k}$ includes the propagation of original pixel noise and artifacts into the background pixels, possible non-ideal background subtraction and the contributions of unidentified faint background sources. $\delta _{{\rm sto},\lambda _k}$ does not contain information on the correlation of the noise among the different pixels or calibration uncertainties.

For each camera, we consider the pixel-dependent uncertainty δ_{sys, j} as

$$\begin{equation} \delta _{{\rm sys},j} =\max \lbrace \epsilon \times I_j^{{\rm obs}},K_j\rbrace , \end{equation} \tag{ D3 }$$

where is the calibration uncertainty for the camera; K is 0 for IRAC, MIPS24, and SPIRE cameras and an estimate of systematic uncertainties for MIPS70,160 and PACS cameras (where such estimation is possible). The main difficulty in estimating is that the cameras are typically calibrated using point sources. We are performing resolved studies of an extended sources, so we have additional uncertainties due to the extended source corrections.

For the IRAC cameras, the standard calibration uncertainty is = 0.05 for point sources (Dale et al. 2005), but we will adopt a larger value to account for the uncertainties in the extended emission correction. The IRAC images are calibrated for point sources. We multiplied the intensities by the asymptotic (infinite radii) value of the aperture correction factors C = 0.91, 0.94, 0.66, 0.74 for the 3.6 μm, 4.5 μm, 5.8 μm, and 8.0 μm bands, respectively. This correction would lead to a correctly calibrated image if the source had constant surface brightness. Since galaxies contain nonuniform emission, the aperture correction factor should, in principle, be different in each region. We multiply the entire image by the factors C, and arbitrarily assign 1/3 of the correction as uncertainty. We use = 0.05 + 1/3 × (1/C − 1) = 0.083, 0.071, 0.221, 0.167 for the 3.6 μm, 4.5 μm, 5.8 μm, and 8.0 μm bands, respectively. These corrections will overestimate the uncertainties in the diffuse regions and underestimate the uncertainties near bright point sources.

For the MIPS, PACS, and the SPIRE cameras we adopt = 0.10. This value is larger than the estimated 4% calibration uncertainty for point sources for MIPS24 (Engelbracht et al. 2007), ≈4% for PACS (Müller et al. 2011), and 7% for SPIRE,³⁷ and comparable to estimates for MIPS70 (Gordon et al. 2007) and MIPS160 (Stansberry et al. 2007). For the MIPS70,160 and PACS cameras, however, δ_sys will be dominated by K in most of the pixels.

For MIPS70 and MIPS160, we can estimate the systematic uncertainties K by comparing the PACS and MIPS resolved photometry. In most of the pixels, the observed differences are larger than expected due to the differences in the camera spectral responses or calibration uncertainties. For PACS70 and PACS160 we can estimate the systematic uncertainties K by comparing the PACS and MIPS resolved photometry For the PACS100 camera we can estimate the systematic uncertainties K by extrapolating the uncertainties estimated for the PACS70 and PACS160 cameras. We proceed as follows.

D.1.1. Systematic Uncertainties for MIPS70 and MIPS160 Cameras

Since for the 70 μm and 160 μm bands the PACS FWHM are smaller than the corresponding MIPS FWHM (see Table 1), in all the common-PSF maps where we employ MIPS data we have PACS data as well. Therefore, we can compare PACS and MIPS images after the image processing. The comparison is straightforward since both images are in the same final-map grid and PSF.

We consider the difference

$$\begin{equation} D_{{\rm PM},j} \equiv \big| I_{{\rm PACS(S)},j}^{{\rm obs}} - I_{{\rm MIPS},j}^{{\rm obs}}\big| , \end{equation} \tag{ D4 }$$

where I^obs_{MIPS, j} and I^obs_{PACS(S), j} are the observed MIPS and PACS (Scanamorphos) flux, in the same band (70 μm or 160 μm). We consider the Scanamorphos data reduction more reliable than the HIPE, so we only employ the Scanamorphos images in the PACS–MIPS comparison. The term D_{PM, j} captures both image artifacts and differences in camera performance and calibration. Unfortunately, artifacts in the PACS cameras will induce uncertainties in the MIPS fluxes and vice versa, but there is no robust, automatic way of isolating the image artifacts. We take

$$\begin{equation} K_{{\rm MIPS},j} = D_{{\rm PM},j}. \end{equation} \tag{ D5 }$$

D.1.2. Systematic Uncertainties for PACS70 and PACS160 Cameras

Based on our dust modeling, we consider the Scanamorphos data reduction to be more reliable than the HIPE-only data reduction. When the MIPS and PACS cameras for a band (70 μm or 160 μm) are used in a final-map PSF (e.g., MIPS70 and PACS70 can be used if the final-map PSF is SPIRE350, SPIRE500, or MIPS160) we proceed in a similar way as for the MIPS bands, taking

$$\begin{equation} K_{{\rm PACS},j} = D_{{\rm PM},j} . \end{equation} \tag{ D6 }$$

When the MIPS camera for a band (70 μm or 160 μm) is not used in a given final-map PSF (e.g., MIPS70 cannot be used if the modeling is done at PACS160 PSF), the comparison between PACS and MIPS is more complex. We start by convolving the PACS image into the corresponding MIPS PSF, and re-binning it into the MIPS native pixel grid, denoted as I^obs_PACS@M. We can compare I^obs_PACS@M and I^obs_MIPS since both images are expressed at the same grid and have the same PSF. We construct _PM, the fractional PACS–MIPS difference (with respect to the PACS image), at the MIPS PSF, as

$$\begin{equation} \epsilon _{{\rm PM},J} = {\big| I_{{\rm PACS(S)@M},J}^{{\rm obs}} - I_{{\rm MIPS},J}^{{\rm obs}}\big| \over \big| I_{{\rm PACS(S)@M},J}^{{\rm obs}} \big|}. \end{equation} \tag{ D7 }$$

We re-bin _PM into the final-map grid, and we take

$$\begin{equation} E_{{\rm PM},j}^{\prime } = \epsilon _{{\rm PM},J} \times \big| I_{{\rm PACS(S)},j}^{{\rm obs}}\big| , \end{equation} \tag{ D8 }$$

where J is the pixel in the PACS@M image grid that contains pixel j. This procedure produces an uncertainty image E_{PM, j} at the final-map grid and PSF, that if degraded into the corresponding MIPS PSF would be ∼|I^obs_PACS(S)@M − I_MIPS^obs|. In this situation, we finally take

$$\begin{equation} K_{{\rm PACS},j} =E_{{\rm PM},j}. \end{equation} \tag{ D9 }$$

D.1.3. Systematic Uncertainties for PACS100 Camera

In the case of PACS100, we do not have any other camera available to perform a direct comparison. We therefore estimate the uncertainties by extrapolating the uncertainties previously calculated for PACS70 and PACS160 as follows. We therefore will assume

$$\begin{equation} K_{{\rm PACS100},j} = \max \lbrace K_{70,j}^{\prime },K_{160,j}^{\prime } \rbrace , \end{equation} \tag{ D10 }$$

where

$$\begin{eqnarray} K_{70,j}^{\prime } &=&{K_{{\rm PACS70},j} \over \big| I_{{\rm PACS(S)70},j}^{{\rm obs}}\big| } \times \big| I_{{\rm PACS(S)100},j}^{{\rm obs}}\big| \end{eqnarray} \tag{ D11 }$$

$$\begin{eqnarray} K_{160,j}^{\prime } &=&{K_{{\rm PACS160},j} \over \big| I_{{\rm PACS(S)160},j}^{{\rm obs}}\big| } \times \big| I_{{\rm PACS(S)100},j}^{{\rm obs}}\big|. \end{eqnarray} \tag{ D12 }$$

Figure 26 shows the resulting S/N estimates for MIPS70 and PACS70 observations of NGC 6946.

Zoom In Zoom Out Reset image size

In order to estimate the model parameter uncertainties, we generate a set r = 1, 2, ..., N_r realizations $I_{\lambda _k,j,r}$ including random noise, for each pixel j, modeling the noise $I_{\lambda _k,j,r}-I_{\lambda _k,j}^{{\rm obs}}$ in pixel j as a sum of three components:

$$\begin{equation} I_{\lambda _k,j,r}=I_{\lambda _k,j}^{{\rm obs}} + \delta I_{\lambda _k,j,r,{\rm ind}}+\delta I_{\lambda _k,r,{\rm cor}}+\delta I_{\lambda _k,j,r,{\rm sys}}. \end{equation} \tag{ E1 }$$

$\delta I_{\lambda _k,j,r,{\rm ind}}$ represents the non-correlated noise component of the pixel j, while $\delta I_{\lambda _k,r,{\rm cor}}$ and $\delta I_{\lambda _k,j,r,{\rm sys}}$ represent the correlated noise component and systematic/calibration uncertainties, assumed to be (completely) correlated among the different pixels of the camera. $\delta I_{\lambda _k,j,r,{\rm ind}}$ is drawn from an independent Gaussian distribution for each pixel, while $\delta I_{\lambda _k,r,{\rm cor}}$ and $\delta I_{\lambda _k,j,r, {\rm sys}}$ are drawn from a global Gaussian distribution, properly normalized for each pixel.

For each camera, we generate a set of N × N_r "pixel" Gaussian variables, with zero mean and variance 1, η_{j, r}, j = 1, 2, ...N, r = 1, 2, ..., N_r, where N is the total number of pixels in the image, and a set of 2N_r "global" Gaussian variables, with zero mean and variance 1, $\psi _1,\psi _2\ldots ,\psi _{N_{r}},\phi _1,\phi _2,\ldots ,\phi _{N_{r}}$.

We set

$$\begin{eqnarray} \delta I_{\lambda _k,j,r,{\rm ind}}&= &\sqrt{1/2} \times \delta _{{\rm sto},\lambda _k} \times \eta _{j,r} \end{eqnarray} \tag{ E2 }$$

$$\begin{eqnarray} \delta I_{\lambda _k,r,{\rm cor}} &=& \sqrt{1/2} \times \delta _{{\rm sto},\lambda _k} \times \psi _r \end{eqnarray} \tag{ E3 }$$

$$\begin{eqnarray} \delta I_{\lambda _k,j,r,{\rm sys}} &=& \delta _{{\rm sys},\lambda _k} \times \phi _r , \end{eqnarray} \tag{ E4 }$$

where $\delta _{{\rm sto},\lambda _k}$ and $\delta _{{\rm sys},\lambda _k}$ were defined by Equations (D2) and (D3). The division of the dispersion $\delta _{{\rm sto},\lambda _k}$ into equal correlated and uncorrelated components $\delta I_{\lambda _k,j,r,{\rm ind}}$ and $\delta I_{\lambda _k,r,{\rm cor}}$ is arbitrary, but seems to capture the main features of the noise. A more sophisticated treatment would use the statistic of the background region to characterize the correlation of the noise in pixels as a function of their separation, but this is beyond the scope of the present work.

For each pixel, the 1σ uncertainty in the measured flux density $\sigma _{\lambda _k,j}$ used to calculate the χ²_j of the model (see Equation (13)) is given by

$$\begin{equation} \sigma _{\lambda _k,j} = \sqrt{ (\delta _{{\rm sto},\lambda _k})^2+ (\delta _{{\rm sys},\lambda _k})^2}. \end{equation} \tag{ E5 }$$

Uncertainties for the inferred (pixel) model parameters are computed using Equation (28). In order to compute the uncertainties in the global quantities, we use Equation (28) on the global quantities, i.e., for each random realization we compute the global inferred quantity and finally we compute the dispersion on the inferred global quantities. This approach preserves the correlation in the pixel noise consistently into the global uncertainties.

In Section 7 we show that maps made at PACS160 resolution (using IRAC, MIPS24, and PACS data only) are less reliable. Furthermore, maps using IRAC and PACS data at MIPS160 resolution provide less accurate³⁸ results than maps computed with IRAC and MIPS data using the same PSF. We therefore proceed to analyze the differences in the PACS and MIPS photometry, to shed light onto the roots of such discrepancies.

Figures 27 and 28 compare the two different PACS data reduction schemes (Scanamorphos and HIPE) with the MIPS camera for NGC 628 and NGC 6946, respectively. All images were convolved into a Gaussian PSF with 64'' FWHM, a PSF broad enough that any camera PSF mismatch should not produce any artifact. The left column shows the PACS/MIPS ratio maps over the galaxy mask. The center columns show an intensity–intensity scatter plot for the galaxy pixels (the pixels in the left column). The right column has a zoom of the −5–30 MJy sr⁻¹ region of the center column plots. We compare PACS70 Scanamorphos and MIPS70, PACS70 HIPE, and MIPS70, PACS160 Scanamorphos and MIPS160, and PACS160 HIPE and MIPS160. In the center and right column, the horizontal and vertical red lines correspond to the 1σ error estimates. By the inclusion of the PACS(S)–MIPS difference term in the uncertainty estimates (see Appendices D.1.1 and D.1.2), the 1σ error estimates will intersect the PACS=MIPS line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Scanamorphos and HIPE pipelines give similar results in the high surface brightness areas, having slightly larger intensities than MIPS (see Figures 27 and 28). The main difference between the two data reduction schemes is in the low surface brightness areas: Scanamorphos intensities are larger than HIPE intensities.

Even though the PACS and MIPS cameras have different spectral response, smooth spectra appropriate to star-forming galaxies should result in nearly the same measured intensity for both MIPS70 and PACS70, and likewise for MIPS160 and PACS160. For our best-fit SEDs for NGC 628 and NGC 6946, we expect PACS70/MIPS70 ∼1.085, and PACS160/MIPS160 ∼1.022 (see Table 6). However, the measured PACS global flux densities do not agree so well with MIPS values, as seen in Table 6. The situation is more critical in resolved studies. Even though the global PACS and MIPS fluxes agree within ≈20%, the ratio maps shows strong discrepancies in the three data sets: the band ratios differ from ≈1 in most of the galaxy.

Figure 29 shows the 70 μm/160 μm color inferred from the different data sets. Again, all images were convolved into a Gaussian PSF with 64'' FWHM. The top row corresponds to NGC 628 and the bottom row to NGC 6946. The left column is the PACS-Scanamorphos color, PACS HIPE is shown in the center column, and, MIPS is shown in the right column. The 70 μm/160 μm colors inferred with the different cameras are also quite different, especially in the outer parts of the galaxies.

Zoom In Zoom Out Reset image size

The colors inferred from PACS(S)70/PACS(S)160 and MIPS70/MIPS160 are similar for most of the galaxy, and thus the inferred dust parameters (U_min) obtained at PACS160 resolution are similar to the ones obtained at MIPS160 resolution. The large PACS(S)160/MIPS160 flux ratio makes the modeling at PACS160 resolution add more dust over all the galaxy, compared to our modeling at MIPS160 resolution, when MIPS constraints are present.

The colors inferred from PACS(H)70/PACS(H)160 and MIPS70/MIPS160 are very different in the outer parts of the galaxy, potentially giving different best-fit parameters in the dust modeling. The HIPE pipeline uses a high-pass filter in the time line observing series before constructing the two-dimensional image map. In studies of HIPE performance using simulated data with realistic noise statistics, it was shown that HIPE intensities are always equal to or smaller than the "true" intensities in the extended low surface brightness regions (M. Sauvage 2011, private communication). When only one SPIRE band (SPIRE250) is included, then the extremely low flux of PACS(H)160 in the low surface brightness regions of the galaxy makes the dust appear to be very cold, and thus increases the amount of dust needed to reproduce the SPIRE250 flux in those regions. We therefore do not recommend to use the HIPE data reduction pipeline for the PACS cameras.

We do not recommend working at PACS160 resolution, unless there is reason to think that high S/N has been achieved. We note, however, that even in regions of moderately high surface brightness, the PACS and MIPS photometry often disagree at the 30% level (see Figures 27 and 28), making dust modeling risky if PACS provides the only data longward of 24 μm. If SPIRE250 data are available, dust modeling can be done at SPIRE250 resolution with results that are reasonably reliable. We note in Figure 14 that total dust masses obtained using IRAC, MIPS24, PACS, and SPIRE250 as constraints are in error by only ≈35%, provided the PACS-Scanamorphos data are used.

We further note that NGC 6946 is a particularly challenging galaxy to model: it has over four decades of dynamic range in the PACS160 band, including an extremely high surface brightness nucleus.