ANDROMEDA'S DUST

The Andromeda galaxy, M31, is the nearest large spiral galaxy. At a distance D = 744 kpc ⁷ (Vilardell et al. 2010), M31 provides an opportunity to study the dust and gas in an external star-forming galaxy with spatial resolution that is surpassed only for the Magellanic Clouds. The structure of the stellar spheroid, disk, and halo of M31 is the subject of ongoing investigations, now being carried out using photometry of large numbers of individual stars (e.g., Dalcanton et al. 2012).

The isophotal major radius is R₂₅ = 95' = 20.6 kpc at 744 kpc (de Vaucouleurs et al. 1991). M31 is classified as a SA(s)b spiral (de Vaucouleurs et al. 1991), but the gas and dust do not conform to a regular spiral pattern. Images of both H i (Braun et al. 2009; Chemin et al. 2009) and infrared emission (Haas et al. 1998; Barmby et al. 2006; Gordon et al. 2006; Fritz et al. 2012; Smith et al. 2012) show structure that appears as much ring-like as spiral in character. The centers of the rings are often offset significantly from the dynamical center of M31.⁸ The off-center ring-like structure has been attributed to a nearly head-on collision with M32 (Block et al. 2006).

Previous studies of far-infrared (FIR) emission from the dust in M31 include maps made with IRAS (Habing et al. 1984; Devereux et al. 1994), the Infrared Space Observatory (ISO) (Haas et al. 1998), and Spitzer (Gordon et al. 2006). The total infrared luminosity was well measured, but the limited angular resolution of IRAS 60 and 100 μm (105'' FWHM), ISO 175 μm (110'' FWHM), and Spitzer 160 μm (39'' FWHM) allowed only a relatively coarse image of the dust distribution.

The present study takes advantage of the high sensitivity and angular resolution of Herschel Space Observatory (Pilbratt et al. 2010) to characterize the dust in M31 on angular scales as small as ∼25'' (=90 pc at 744 kpc), and at wavelengths as long as 500 μm, thereby capturing the emission from whatever cold dust may be present. The Herschel Exploitation of Local Galaxy Andromeda (HELGA; Fritz et al. 2012; Smith et al. 2012) recently employed Herschel imaging to study the dust distribution in M31, using single-temperature modified blackbody fits to the 70–500 μm emission from each pixel. The present study differs from HELGA in two ways. First, we use an independent set of Herschel observations (Groves et al. 2012; O. Krause et al., in preparation), somewhat deeper than those obtained by HELGA. Second, we use a physical dust model (Draine & Li 2007) to model the spectral energy distribution (SED) from 6–500 μm, and to estimate the dust mass surface density, intensities of starlight heating the dust, and the polycyclic aromatic hydrocarbon (PAH) abundance, using methods recently developed by Aniano et al. (2012) for studying the galaxies in the KINGFISH sample (Kennicutt et al. 2011).

The organization of the paper is as follows: the observational data are described in Section 2, and the dust-model fitting methods are outlined in Section 3. Results are presented in Sections 4–8, where we estimate the total dust luminosity and mass, the dust/gas ratio, the metallicity as a function of radius, the characteristics of the starlight heating the dust, and the PAH abundance. Evidence for variation of the dust properties is discussed in Section 9. The results are summarized in Section 10.

Appendix A examines the inconsistency between PACS and MIPS photometry of M31. The starlight contribution from the M31 bulge stars is calculated in Appendix B.

M31 has been mapped using the IRAC (Fazio et al. 2004) and MIPS (Rieke et al. 2004) cameras on Spitzer Space Telescope (Werner et al. 2004). More recently, maps have been obtained by the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) cameras on Herschel (Pilbratt et al. 2010).

The present analysis uses IRAC data from Barmby et al. (2006)⁹ and MIPS data from Gordon et al. (2006).¹⁰ The PACS and SPIRE imaging was done in parallel mode at medium scan speed (20'' s⁻¹) for a total time of ∼24 hr (O. Krause et al. 2013, in preparation), and are the same data as used by Groves et al. (2012), except for small changes in calibration. We use the most recent calibrations of PACS and SPIRE.¹¹ We refer to the different images by the camera name and nominal wavelength in microns: IRAC3.6, IRAC4.5, IRAC5.8, IRAC8.0, MIPS24, MIPS70, MIPS160, PACS70, PACS100, PACS160, SPIRE250, SPIRE350, and SPIRE500.

The center of M31 is located at α₂₀₀₀ = 10685, δ₂₀₀₀ = 41269 (Crane et al. 1992), or Galactic coordinates ℓ = 121°, b = −22°. Because M31 is only 22° from the Galactic plane, removal of the Galactic foreground "cirrus" emission is challenging, particularly in view of the large angular extent of M31. We are helped by the high inclination i ≈ 78° of M31, which raises the surface brightness, and increases the contrast with foreground and background emission.

Subtraction of foreground and background emission has been carried out following methods described in Aniano et al. (2012), with automatic identification of background pixels and fitting of a "tilted plane" background model (with three parameters—zero point, tilt, and tilt orientation) for all bands except IRAC5.8 and IRAC8.0. For IRAC5.8 and IRAC8.0, it was found that the simple "tilted plane" background model left obvious large-scale residuals, presumably a consequence of the large angular extent of M31. For IRAC5.8 and IRAC8.0, it was found necessary to use more complex curved surfaces, rather than simple tilted planes, to model the background. Even so, the background estimation for IRAC5.8 and IRAC8.0 seems to be less successful than for other bands, for reasons that are further discussed in Section 8.

The foreground cirrus has significant structure on ∼1° scales (see, e.g., the IRAS 100 μm or LAB 21 cm map of the area around M31), and the simple tilted plane model (or low-order curved surface in the case of IRAC5.8 and IRAC8.0) will not remove all of the foreground cirrus. This will be problematic in the low surface brightness outer regions of M31. However, after the foreground- and background-subtracted FIR images have been fitted by dust models, comparison of the resulting dust map with maps of the H i at the radial velocity of M31 will allow us to assess the effects of imperfect foreground subtraction. Except for IRAC5.8 and IRAC8.0, the present tilted plane foreground model appears adequate for R ≲ 25 kpc (see Section 5).

For modeling the dust in resolved systems, it is essential that multiband imaging be convolved to a common point-spread function (PSF). The present study is carried out at two angular resolutions: that of the SPIRE350 camera (FWHM = 249 = 90 pc at 744 kpc, henceforth referred to as S350 resolution) and that of the MIPS160 camera (FWHM = 39'' = 141 pc at 744 kpc, henceforth referred to as M160 resolution). Image convolutions were carried out using kernels obtained as described by Aniano et al. (2011). For studies at S350 resolution, we are unable to use SPIRE500 or MIPS160 imaging, because the PSF of those cameras is too broad. The convolved images at S350 resolution were sampled with 10'' × 10'' pixels. Images convolved to M160 resolution were sampled with 18'' × 18'' pixels.

Figure 1 shows SPIRE350 images of M31; the left image is at S350 resolution; the image on the right is convolved to M160 resolution. The lower resolution of M160 is evident when the images are compared, but most details of the image that are seen at S350 resolution are also identifiable at M160 resolution. Analysis of M31 at M160 resolution—benefiting from the improved signal/noise of the larger M160 pixels as well as being able to use both SPIRE500 and MIPS160 data to constrain the models—will, therefore, not lose important structural features.

Zoom In Zoom Out Reset image size

M31 is highly inclined. Using DRAO H i 21cm data, Chemin et al. (2009) estimated an average inclination i = 743 and position angle (P.A.) = 377 for the disk at R = 6–27 kpc. From WSRT H i 21 cm observations, Corbelli et al. (2010) obtain an inclination i = 777 and P.A. = 377 for R = 10–13 kpc. Corbelli et al. (2010) showed that the H i kinematics can be described using circular rings with centers offset from the dynamical center of M31, and with inclinations between 75° and 79° for 10 < R < 25 kpc).

We will assume the disk to have a single inclination i = 777 and P.A. = 377. Figure 1 shows ellipses with major radius 5500'' corresponding to R = 19.8 kpc at D = 744 kpc, with i = 777 and P.A. = 377. Radial trends will be studied by averaging in annuli centered on the location of the presumed supermassive black hole, α_J2000 = 106847, δ_J2000 = 412690 (Crane et al. 1992).

The S350 PSF provides excellent spatial resolution (FWHM = 90 pc along the major axis). We will compare results obtained at S350 resolution with results obtained using the M160 PSF. While the dust models at M160 resolution smooth out dust structures smaller than ∼140 pc along the major axis (or ∼700 pc along the minor axis), working at M160 resolution permits (1) extending the wavelength coverage to 500 μm (by including SPIRE500 photometry), and (2) using MIPS160 photometry, which is somewhat deeper than the PACS160 imaging. In addition, the larger pixels provide improved signal/noise in low surface brightness regions. At S350 resolution, the models are constrained by photometry in 11 bands (4 IRAC, 2 MIPS, 3 PACS, 2 SPIRE). At M160 resolution, MIPS160 and SPIRE500 are added, giving a total of 13 photometric constraints for each pixel. Comparison of model results at S350 and M160 resolution will provide insight into the reliability of the method.

3.1. Dust Model

3.2. Renormalization of the DL07 Dust Mass

3.3. Dust Masses and Starlight Intensities

The present study seeks to estimate dust masses in M31 by modeling the observed FIR and submillimeter emission. The dust is assumed to be heated by starlight, and the modeling has the freedom to adjust the starlight intensities such that the grains are heated to temperatures such that their emission spectrum is consistent with the observed shape of the SED. Then, with the dust temperatures set by the starlight heating rates, the dust mass is proportional to the observed emission. Note that there is no single temperature characterizing the dust—even in a single radiation field, grains of different size and composition have different temperatures, and the very small grains undergo temperature fluctuations due to single-photon heating. Additionally, a single "pixel"—which may be several hundred parsecs in transverse dimension—may include subregions with different starlight intensities.

For a given starlight intensity U, grain size a and composition, we solve for the temperature probability distribution (dP/dT)_{U, a, comp}. For large grains dP/dT can be approximated by a δ-function, but for small grains dP/dT can be very broad, and must be solved for as described by Li & Draine (2001). The time-averaged emission spectrum for a grain is then

$$\begin{equation} (p_\nu)_{U,a,{\rm comp}} = 4\pi \int dT \, (dP/dT)_{U,a,{\rm comp}} \, (C_{{\rm abs},\nu })_{a,{\rm comp}} B_\nu (T), \end{equation} \tag{ 1 }$$

where (C_{abs, ν})_{a, comp} is the absorption cross section at frequency ν. The emission per unit mass for a dust mixture exposed to starlight U is

$$\begin{equation} \left(\frac{dL_\nu }{dM_d}\right)_U = \frac{\sum _{\rm comp} \int da (dn/da)_{\rm comp} (p_\nu)_{U,a,{\rm comp}}}{\sum _{\rm comp}\int da (dn/da)_{\rm comp} (4\pi /3)a^3 \rho _{\rm comp}}, \end{equation} \tag{ 2 }$$

where ρ_comp is the solid density, and (dn/da)_comp is the size distribution. The luminosity of a region j is obtained by summing over the starlight distribution:

$$\begin{equation} L_{\nu,j} = \int dU \frac{dM_{{\rm d},j}}{dU} \left(\frac{dL_\nu }{dM_d}\right)_U. \end{equation} \tag{ 3 }$$

For a region j of solid angle Ω_j, with dust mass M_{d, j} = Σ_{Md, j}Ω_jD², the dust mass dM_{d, j} exposed to starlight intensities in [U, U + dU] is assumed to be given by the simple parameterization proposed by DL07:

$$\begin{eqnarray} \frac{dM_{{\rm d},j}}{dU} &=& \Sigma _{M{\rm d},j}\Omega _j D^2 \nonumber \\ &&\times \left[ (1-\gamma _j)\delta (U-U_{{\rm min},j}) + \gamma _j \frac{(\alpha _j-1)U^{-\alpha _j}}{U_{{\rm min},j}^{1-\alpha _j} -U_{\rm max}^{1-\alpha _j}} \right],\nonumber \\ \end{eqnarray} \tag{ 4 }$$

where Σ_{Md, j} is the total dust mass surface density in region j, δ is the Dirac δ-function, and α_j > 1 is a power-law index characterizing the distribution of starlight intensities. A fraction (1 − γ_j) of the dust mass is heated by starlight intensity U = U_{min, j}, with the remaining fraction γ_j exposed to starlight with intensities U_{min, j} < U ⩽ U_max, with a power-law distribution $dM_{\rm d}/dU\propto U^{-\alpha _j}$. For NGC 628 and NGC 6946, Aniano et al. (2012) found that the parameter U_max could be fixed at U_max = 10⁷ without significantly degrading the quality of the fits, hence we also fix U_max = 10⁷. Thus, for each region j, we have five adjustable parameters characterizing the dust: {Σ_{Md, j}, q_{PAH, j}, U_{min, j}, α_j, γ_j}.¹³ We require Σ_{Md, j} ⩾ 0, 0 ⩽ γ_j ⩽ 1, and 1 < α_j < 3. We use one additional parameter to characterize the contribution of direct starlight (taken to have the spectrum of a 5000 K blackbody) to the photometry for pixel j. Thus each region has 6 adjustable parameters, and either 11 or 13 data, depending on whether we use S350 or M160 resolution. Because MIPS70 and PACS70 cover essentially the same wavelengths, and likewise for MIPS160 and PACS160, the effective number of constraints is 10 (at S350 resolution) or 11 (at M160 resolution).

The fitting procedure assumes the DL07 model dust to be heated by starlight with the solar-neighborhood spectrum found by Mathis et al. (1983). In Sections 7 and 8 we examine the effects of changes in the starlight spectrum when the bulge stars make a significant contribution to the dust heating.

3.4. Systematic Uncertainties in Dust Mass Estimation

3.5. Single Pixels versus Annuli

Determination of the dust mass requires that the shape of the dust SED be well-determined, so that the dust temperatures can be constrained. With the signal-to-noise properties of the present observations, we do not attempt to estimate the dust mass in a single 10'' × 10'' pixel unless the dust luminosity per unit area (on the sky plane) Σ_Ld > Σ_{Ld, min} = 1.5 × 10⁶ L_☉ kpc⁻² (IR intensity I > 4.8 × 10⁻⁵ erg s⁻¹ cm⁻² sr⁻¹). The irregular contour in Figure 1 and the other images bounds the contiguous region with Σ_Ld > Σ_{Ld, min} (there are additional pixels with Σ_Ld > Σ_{Ld, min} outside this contour).¹⁴ For Σ_Ld ≲ Σ_{Ld, min} the single-pixel dust mass estimates are quite uncertain, especially at S350 resolution where the pixels are smaller, and SPIRE500 and MIPS160 cannot be used. At M160 resolution the larger (18'' × 18'') pixels and inclusion of SPIRE500 and MIPS160 help, but even at M160 resolution the single-pixel dust mass estimates near Σ_{Ld, min} = 1.5 × 10⁶ L_☉ kpc⁻² have substantial uncertainties.

To study radial dependences, we will sometimes take the dust properties obtained by single-pixel modeling and average them over annuli with widths 10'' along the minor axis, and ΔR = 469 (169 pc) along the major axis. While these annuli are not resolved along the minor axis (the MIPS160 FWHM = 39''), they are resolved along the major axis.¹⁵

In order to study the dust at R > 17 kpc, where a large fraction of the pixels have Σ_Ld < Σ_{Ld, min} (see Figure 2(c)), we use the background-subtracted images to extract the flux in each band from larger concentric annuli with widths 40'' along the minor axis, and ΔR = 188'' (677 pc) along the major axis. These ΔR = 677 pc annuli are fully resolved, even by MIPS160. The dust mass and starlight intensity distribution in each annulus are then estimated by fitting the SED of the annulus. By integrating over the many pixels in each annulus, the random noise is substantially suppressed. We will see that reliable estimates of the mean dust mass surface density are possible out to R ≈ 25 kpc.

Zoom In Zoom Out Reset image size

As described in Aniano et al. (2012), the background-subtracted multiwavelength images are first used to estimate Σ_Ld, the total dust luminosity per projected area on the plane of the sky. The luminosity is obtained by fitting a DL07 dust model to the observed fluxes from the pixel, and then integrating over the model SED. While the model dust mass in individual pixels is unreliable for Σ_Ld ≲ Σ_{Ld, min}, the pixel luminosity estimate is reliable down to considerably lower surface brightness. Figure 2(a) shows the dust luminosity surface density Σ_Ld at M160 resolution. From the image it can be seen by eye that the dust luminosity is detectable down to Σ_Ld ≈ 10^5.9 L_☉ kpc⁻².

The mean luminosity per area as a function of radius is found by averaging the single-pixel luminosities over the ΔR = 177 pc annuli. Figure 2(b) shows Σ_Ldcos i, the dust luminosity surface density on the disk plane, averaged over the ΔR = 177 pc annuli, and plotted against the annular radius R. Triangles show annular averages of single-pixel results obtained at M160 resolution (using all cameras) and at S350 resolution (not using SPIRE500 or MIPS160). Circles show results of modeling the SED of ΔR = 677 pc annuli. The smooth decline in Σ_Ld with increasing R out to R ≈ 23 kpc is an indication that background subtraction has been reasonably successful at least down to surface brightnesses Σ_Ld ≈ 4 × 10⁵ L_☉ kpc⁻² (Σ_Ldcos i ≈ 8 × 10⁴ L_☉ kpc⁻²).

Figure 2(c) shows the fraction f of the pixels at each R that have Σ_Ld > Σ_{Ld, min}. The fraction f = 1 for R < 15 kpc, but at R ≈ 15.6 kpc the fraction begins to fall. At R ≈ 24 kpc the fraction has fallen to f ≈ 0.1; some fraction of these pixels may be raised above Σ_{Ld, min} by unresolved background galaxies.

The dust luminosity surface density Σ_Ld derived from the SEDs of ΔR = 677 pc annuli, shown in Figure 2(b), is in good agreement with Σ_Ld obtained from pixel-based modeling. From the annular SEDs, the dust luminosity surface density for 17 kpc < R < 25 kpc is

$$\begin{eqnarray} \Sigma _{L{\rm d}}\cos i &\approx & 8.0\times 10^5\exp [-(R-17\,{\rm kpc})/3.1\,{\rm kpc}]\, L_{\odot }\,{\rm kpc}^{-2}.\nonumber \\ \end{eqnarray} \tag{ 5 }$$

For purposes of estimating the integrated dust properties from the pixel-based modeling, we sum over the M160 pixel-based modeling for R < 17 kpc, and use the annular SED modeling for R > 17 kpc (see Table 4). We estimate the R < 25 kpc dust luminosity to be L_d = 4.26 × 10⁹ L_☉.

Because νL_ν peaks near 100 μm, the global dust luminosity of M31 was already reasonably well-determined by the low-resolution observations by IRAS (60, 100 μm) and COBE (140, 240 μm). Adding 175 μm imaging by ISO, Haas et al. (1998) obtained L_d ≈ 3.8 × 10⁹ L_☉.¹⁶ Using MIPS photometry, Gordon et al. (2006) found L_d ≈ 4.0 × 10⁹ L_☉, consistent, within the uncertainties, with the present value of L_d = 4.26 × 10⁹ L_☉.

There is considerable structure in the dust luminosity image in Figure 2(a), with close similarities to the SPIRE 350 μm image in Figure 1. While the structure is not purely "circular" (on the plane of the galaxy), the dust luminosity density extracted in (inclined) circular annuli centered on the dynamical center, shown in Figure 2(b), exhibits a strong central peak (due to small amounts of relatively warm dust), two clear rings (at R = 5.6 kpc and R = 11.2 kpc) and indications of a third ring (at R ≈ 15.1 kpc). The 11.2 kpc ring was evident in the original IRAS imaging (Habing et al. 1984), and all three rings were noted in the ISO 175 μm imaging by Haas et al. (1998).

5.1. Dust Mass

Before the launch of Herschel Space Observatory, the mass of dust in M31 had been estimated by Xu & Helou (1996), Haas et al. (1998), and Gordon et al. (2006) (see Table 5). With Herschel imaging out to 500 μm, we are now able to obtain improved estimates for the total dust mass, and to map the dust out to ∼20 kpc. Figure 3(a) shows the dust mass surface density Σ_Md on the plane of the sky, obtained by modeling at M160 resolution.

Zoom In Zoom Out Reset image size

Table 5. Estimates for the Dust Mass in M31

Md(107 M☉)a	Data used	Reference
2.2 ± 0.7	IRAS, extinction, H i	Xu & Helou (1996)
3.5 ± 1.0	IRAS100, COBE-DIRBE 140,240, ISO 175 μm	Haas et al. (1998)
4	IRAS, COBE-DIRBE, ISO, MIPS	Gordon et al. (2006)
5.05 ± 0.45	PACS 100,160; SPIRE 250,350,500; annuli	HELGA I Fritz et al. (2012)
2.6	PACS 100,160; SPIRE 250,350,500; S500 pixels	HELGA II Smith et al. (2012)
5.4 ± 1.1	MIPS, PACS, SPIRE; M160 pixels and annuli	This work

Note. ^aFor D = 744 kpc.

Download table as:  ASCIITypeset image

In addition to constructing maps showing the dust properties in every pixel, we also average the dust properties over rings that are circular on the disk plane in order to study radial trends. For each annulus, we calculate the mean dust mass surface density obtained by summing the dust masses in pixels with Σ_Ld > Σ_{Ld, min}, dividing by the total area of the annulus. For R < 15.6 kpc all pixels have Σ_Ld > Σ_{Ld, min}, but this procedure will underestimate the actual mean surface density for R > 16 kpc. Figure 3(b) shows the dust mass surface density projected onto the M31 disk plane, as a function of galactocentric radius R. Triangles show annular averages of single-pixel modeling at S350 resolution (MIPS160 and SPIRE500 not used) and at M160 resolution (using all the data). Dust mass surface densities obtained at S350 resolution (not using MIPS160 and SPIRE500) agree to within ∼10% with dust mass estimates obtained at M160 resolution, using all cameras.

Dust mass surface densities Σ_Md estimated by fitting SEDs of ΔR = 677 pc annuli are shown (open circles) in Figure 3(b), and are seen to agree to within ∼10% with the results of the pixel-based modeling for R < 15 kpc where the signal-to-noise is high. For R > 16 kpc, where a growing fraction of pixels has Σ_Ld < Σ_{Ld, min}, the annular SED is the best way to estimate the dust mass. The dust mass surface density at R > 17 kpc is seen to be approximated by the broken line in Figure 3(b):

$$\begin{eqnarray} \Sigma _{M{\rm d}}\cos i \approx 2.2&\times & 10^4\exp [-(R-17\,{\rm kpc})/3.8\,{\rm kpc}]\,M_{\odot }\,{\rm kpc}^{-2}\nonumber \\ & & {\rm for}\quad R>17\,{\rm kpc}. \end{eqnarray} \tag{ 6 }$$

The dust mass surface density in Figure 3(b) shows two distinct peaks, corresponding to rings at R = 5.6 kpc and R = 11.2 kpc, with a third ring at R = 15.1 kpc. The 11.2 kpc ring coincides with a peak in the H i 21 cm emission, a ring of star formation (seen in Hα; Devereux et al. 1994), and a ∼40% overdensity of stars with ages >1 Gyr (Dalcanton et al. 2012). The image of the dust mass surface density in Figure 3(a) shows a conspicuous deficiency of dust between ∼16 and 20 kpc on the SW side of the disk. As will be discussed further below, this is not an artifact—the H i gas shows a similar deficiency in the same region.

The DL07 dust model has¹⁷

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

At R ≈ 20 kpc, the dust mass surface density (projected onto the disk of M31) Σ_Mdcos i ≈ 9 × 10³ M_☉ kpc⁻², corresponds to a visual extinction A_V ≈ 0.07 mag normal to the disk of M31, or A_V ≈ 0.3 mag along our line-of-sight. The peak dust mass surface density occurs at R = 11.2 kpc with a dust surface density, projected onto the M31 disk, Σ_Mdcos i ≈ 7 × 10⁴ M_☉ kpc⁻², corresponding to A_V ≈ 0.5 mag normal to the disk (A_V ≈ 2.4 mag along our line-of-sight). Many pixels in Figure 3 have higher dust surface densities, reaching Σ_Md ≈ 10^6.0 M_☉ kpc⁻², corresponding to A_V ≈ 7 mag. In these pixels the dust is presumably distributed inhomogeneously. Hubble Space Telescope observations of M31 (Dalcanton et al. 2012) may allow measurement of reddening toward many individual stars in M31; it will be of great interest to compare the stellar reddening values with the present maps of dust surface density.

We consider our best estimate for the dust mass to be the dust mass obtained using all cameras (including MIPS160). For R < 17 kpc we use the M160 resolution single-pixel estimates for Σ_Md, while for R > 17 kpc we use the annular SEDs to estimate the dust mass in each annulus. We find a total dust mass M_d = 5.4 × 10⁷ M_☉ within R = 25 kpc (see Table 4).

It is difficult to estimate objectively the uncertainty in the estimate for M_d. Calibration uncertainties in the photometry itself are at least 10% for each of the MIPS, PACS, and SPIRE bands. The inconsistencies between MIPS and PACS at 70 μm and 160 μm are larger than expected (see Appendix A), and, in addition, there are difficult-to-assess errors arising from assumptions made in the modeling about the physical properties of the dust, and simplified treatments of the starlight intensity distribution. Overall, we adopt a tentative uncertainty estimate of 20% for the global dust mass given in Tables 4 and 5.

Figure 4 shows the cumulative dust luminosity and dust mass as a function of R. We find that half of the dust mass in M31 lies at R > 12.3 kpc, and >10% lies beyond 19 kpc.

Zoom In Zoom Out Reset image size

The HELGA survey (Fritz et al. 2012; Smith et al. 2012) obtained PACS and SPIRE imaging of M31 at a scan speed of 60'' s⁻¹ (three times faster than the scan speed for the PACS and SPIRE imaging used in the present study). HELGA I (Fritz et al. 2012) measured the fluxes in a central circular aperture and five concentric annuli extending out to 20 kpc. Except for the central circle, the annular boundaries were ellipses. The λ ⩾ 100 μm SED of each annulus was fit with a modified blackbody F_ν∝ν^βB_ν(T), with fixed β = 1.8. The total dust mass M_d = (5.05 ± 0.45) × 10⁷ M_☉ found within 20 kpc by HELGA I is close to the present estimate M_d(R < 20 kpc) = (5.0 ± 1.0) × 10⁷ M_☉.

HELGA II (Smith et al. 2012) restricted their study to pixels satisfying their criterion of 5σ detection in each of six bands (MIPS70, PACS100, PACS160, and the SPIRE bands). The MIPS70 photometry was used only as an upper limit: a modified blackbody Aν^βB_ν(T) was used, with A, β, and T adjusted to fit the 100–500 μm photometry for each pixel; models for the cool dust were acceptable provided they did not exceed the MIPS70 photometry. In the regions included in their study, HELGA II found a total dust mass M_d = 2.6 × 10⁷ M_☉, but Smith et al. (2012) note that the pixels satisfying their 5σ criterion accounted for only ∼50% of the global 500 μm emission from M31.

5.2. Dust-to-gas Ratio and Metallicity

Braun et al. (2009) mapped the H i 21 cm emission out to ∼25 kpc from the center of M31. After correcting for self-absorption, they obtained a map of the H i surface density, which we use here. Integrating this map out to R = 25 kpc yields an H i mass $M({{\rm H\,\scriptsize{I}}})=6.38\times 10^9\,M_{\odot }$ (88% of the total H i in their map). CO (1–0) has been mapped out to ∼12 kpc, with the associated H₂ mass estimated¹⁸ to be M(H₂) ≈ 3.4 × 10⁸ M_☉ (Nieten et al. 2006). We take the total H surface density to be $\Sigma _{\rm H}\approx \Sigma ({\rm H\,\scriptsize{I}})+\Sigma ({\rm H}_2)$, where Σ(H₂) is obtained from the observed CO 1–0 emission assuming a constant X_CO = 2 × 10²⁰ H₂ cm⁻²( K km s⁻¹)⁻¹.

We may be underestimating Σ_H. H₂ that is "CO dark" (i.e., not associated with CO 1–0 emission) is assumed to be a small fraction of the total H₂ mass. More importantly, we do not include H ii gas in our estimate for Σ_H. The ionized gas mass in bright H ii regions is small, but the mass in low-density diffuse H ii may not be negligible. For our Galaxy, it is estimated that diffuse H ii accounts for ∼23% of the total ISM mass at R < 20 kpc (Draine 2011b). The center of M31 has Σ_H ≈ 1.8 × 10⁶ M_☉ kpc⁻² averaged over R < 1 kpc. The extinction-corrected surface brightness of Hα ∼2 × 10⁻⁵ erg cm⁻² s⁻¹ sr⁻¹ (Devereux et al. 1994; Tabatabaei & Berkhuijsen 2010) corresponds to an H ii surface density $6\times 10^6 (\,{\rm cm}^{-3}/\langle n_e\rangle) T_4^{0.94} \,M_{\odot }\,{\rm kpc}^{-2}$, where 〈n_e〉 is the electron density within the emitting regions. The density-sensitive [S ii]6716/[S ii]6731 line ratio indicates 〈n_e〉 > 1.2 × 10² cm⁻³ for R ≲ 0.2 kpc (Ciardullo et al. 1988), hence the mass in ∼10⁴ K plasma appears to be small compared to that in neutral gas. X-ray observations indicate that R ≲ 1 kpc is filled with hot gas with T ≈ 4 × 10⁶ K and a mean surface density ∼2 × 10⁵ M_☉ kpc⁻² (Bogdán & Gilfanov 2008), only ∼10% of the R < 1 kpc H i mass.

It is also possible that we may have overestimated Σ_H in the center—Leroy et al. (2011) find that X_CO in the center (R < 2 kpc) of M31 is lower by about a factor two. However, there is little CO emission at R ≲ 2 kpc (Nieten et al. 2006), hence the total gas mass there is insensitive to uncertainties in X_CO.

Figure 5 shows Σ_H at S350 and M160 resolution. Note the deficiency of gas in the SW for R > 15.6 kpc; the same region is also seen (see Figure 2) to be deficient in emission from dust. The reason for the deficiency in dust and gas in the SW corner of M31 (between R = 16 and R ≈ 20 kpc) is not known, but this part of the M31 disk has the appearance of having been affected by some recent event. Block et al. (2006) argued that a nearly head-on collision with M32 ∼210 Myr ago can account for the observed offset of the center of the 11 kpc ring, but such an encounter would not seem likely to produce the observed deficiency of H i and dust in the SW at R ≈ 16–20 kpc. It may be the result of an encounter with another member of the Local Group within the past gigayears. Lewis et al. (2013) discuss the Giant Stellar Stream, extending from the SW side of the disk toward M33. H i 21 cm observations also show a diffuse gaseous filament connecting M31 and M33, including a feature extending from the SW side of the M31 disk toward M33. The origin of the Giant Stellar Stream and the H i filamentary structure is uncertain, but strongly suggestive of a close passage of M33, possibly accounting for the deficiency of interstellar matter on the SW side of the disk, ∼18 kpc from the center of M31.

Zoom In Zoom Out Reset image size

Figure 6(a) is a map of the dust-to-H mass ratio at M160 resolution. Some azimuthal structure is evident, but the main feature is a conspicuous radial trend, with the dust-to-H ratio peaking at the center and declining with R. Figure 6(b) shows results for every pixel (at M160 resolution) satisfying the condition Σ_Ld > Σ_{Ld, min}. Because of the low signal/noise ratio (S/N) in individual pixels, the scatter in the derived M_d/M_H is quite pronounced for R ≳ 17 kpc. However, the mean dust/H ratio in each annulus exhibits a clear trend, decreasing with increasing R. The central dust-to-H ratio is M_d/M_H ≈ 0.028, declining to ∼0.005 at R ≈ 20 kpc: a factor of ∼5 change in the dust-to-H ratio moving from the center to R ≈ 20 kpc.

Zoom In Zoom Out Reset image size

Figure 7 (triangles) shows the radial profile of the dust-to-H ratio estimated at M160 resolution. Each triangle in Figure 7 is obtained by summing the M160 resolution dust and gas within rings for R < 23 kpc, including only pixels satisfying the criterion Σ_Ld > Σ_{Ld, min}, and calculating the ratio of (total dust)/(total gas) for that ring. The dust/gas ratio therefore applies only to the pixels satisfying the Σ_Ld > Σ_{Ld, min} cut. From Figure 2(b) we see that essentially 100% of the pixels at R < 16 kpc satisfy the surface brightness cut. However, at R = 20 kpc, ∼50% of the pixels in the annulus do not satisfy the cut. This incompleteness is due in part to the general radial decline in surface brightness, but in part it reflects the conspicuous deficit of both gas and dust in the SW corner of M31. Also shown in Figure 7 (circles) is the dust-to-H mass ratio estimated for ΔR = 677 pc annuli, where the dust mass is obtained by modeling the SED of each annulus, and the gas mass is the total mass of gas within the annulus. The dust-to-gas ratio is well-behaved out to 25 kpc. For R < 18 kpc the dust-to-H mass ratio estimates for the ΔR = 677 pc annuli at R = 16–23 kpc are in good agreement with the dust-to-H mass ratios determined only for the pixels with Σ_Ld > Σ_{Ld, min}. For R > 18 kpc the dust/gas ratios for the pixels with Σ_Ld > Σ_{Ld, min} is somewhat higher than the result from the annular photometry. This bias is attributable to the fact that for Σ_Ld just below Σ_{Ld, min}, noise can raise the pixel above the threshold, leading to overestimation of the dust mass. In this regime, dust mass estimates from the annular photometry and modeling should be more reliable.

Zoom In Zoom Out Reset image size

We find that the dust/H ratio for 0–25 kpc declines monotonically with increasing R, with

$$\begin{equation} \frac{M_{\rm d}}{M_{\rm H}} \approx \left\lbrace \begin{array}{@{}l l}0.0280 \exp (-R/8.4\,{\rm kpc}) & R< 8\,{\rm kpc}\\ \ms 0.0165 \exp (-R/19\,{\rm kpc}) & 8\,{\rm kpc}< R < 18\,{\rm kpc}\\ \ms 0.0605 \exp (-R/8\,{\rm kpc}) & 18\,{\rm kpc}< R \lesssim 25\,{\rm kpc},\\ \end{array} \right. \end{equation} \tag{ 8 }$$

shown as a dashed line in Figure 7. Also shown in Figure 7(a) is the dust/gas ratio estimated by HELGA I (Fritz et al. 2012) for each of 5 radial zones. The HELGA I results are in fair agreement with our findings interior to ∼8 kpc, but at larger radii tend to exceed our dust mass estimates.

HELGA II (Smith et al. 2012) conclude that the dust/H ratio follows an exponential profile in M31, with M_dust/M_H ≈ 0.049exp (− R/8.7 kpc); this is plotted in Figure 7(a). Our central dust/H ratio ∼0.028 is only 60% of the value found by HELGA II.

If depletions are similar to the local Milky Way, we expect

$$\begin{equation} \frac{M_{\rm d}}{M_{\rm H}} \approx 0.0091 \frac{Z}{Z_\odot }, \end{equation} \tag{ 9 }$$

where Z is the mass fraction of elements other than H and He. This allows us to estimate the metallicity from our measured dust/H mass ratio (8):

$$\begin{equation} \frac{Z}{Z_\odot } \approx \left\lbrace \begin{array}{@{}l l}3.08 \exp (-R/8.4\,{\rm kpc}) & R< 8\,{\rm kpc}\\ \ms 1.81 \exp (-R/19\,{\rm kpc}) & 8\,{\rm kpc}< R < 18\,{\rm kpc}\\ \ms 6.65 \exp (-R/8\,{\rm kpc}) & 18\,{\rm kpc}< R \lesssim 25\,{\rm kpc}. \end{array} \right. \end{equation} \tag{ 10 }$$

This relation is plotted in Figure 7(b).

Zurita & Bresolin (2012, hereafter ZB12) measured elemental abundances in M31 H ii regions using "direct" methods between 8 and 16 kpc. When they allow for depletion of oxygen into grains and a bias against H ii regions with high oxygen abundances, they estimate that

$$\begin{eqnarray} {\rm (O/H)} &\approx & 1.8 {\rm (O/H)}_\odot \exp (-R/19\,{\rm kpc}). \end{eqnarray} \tag{ 11 }$$

However, it is important to note that the H ii region abundance determinations have appreciable uncertainties, and do not agree well with metallicity determinations in the atmospheres of B supergiants. Our result (Equation (10)) for the metallicity is in excellent agreement with the ZB12 metallicity over the 8 kpc < R < 16 kpc range where the ZB12 H ii region metallicities were based on direct determinations of the gas temperature, and are therefore most reliable.

According to Equation (10), the ISM in M31 has supersolar abundances for R ≲ 11 kpc. This is consistent with the high WC/WN Stellar ratio observed in M31 (Neugent et al. 2012).

If the distribution of both stars and dust were known, the intensity and spectrum of the starlight heating the dust could be obtained from the equations of radiative transfer (see, e.g., Popescu & Tuffs 2013). In dusty star-forming galaxies this is a formidable problem, because of the complex and correlated spatial distributions of both stars and dust.

To model the infrared emission from dust, we take a much simpler approach, empirical approach to the starlight heating. Within a single "pixel" (which may include ∼10⁴  pc² of disk area) the dust may be exposed to a wide range of starlight intensities, ranging from the general starlight background in the diffuse ISM, to high intensities found in star-forming regions. The DL07 model adopts a parameterized distribution function for the starlight heating rate: the dust within a given pixel is assumed to be subject to starlight heating rates ranging from U = U_min to a peak value U_max = 10⁷, with an intensity distribution given by Equation (4). When fitting the dust model to the data, we in effect use the dust grains as photometers to determine the intensity of starlight in the regions where dust is present. The parameter U_min is interpreted as being the starlight heating rate in the diffuse ISM. The mean (weighted by dust mass) starlight heating rate within a pixel is 〈U〉.

Figure 8(a) is a map of U_min for M31. The U_min parameter is strongly peaked at the center. Figure 8(b) shows a generally smooth decline of U_min with increasing galactocentric radius R, declining from U_min ≈ 25 in the central 200 pc (at S350 resolution) to U_min ≈ 0.2 at R = 16 kpc. Beyond 16 kpc, U_min obtained from single-pixel modeling (triangles) begins to differ from U_min obtained from annular photometry (circles).

Zoom In Zoom Out Reset image size

For starlight with the interstellar radiation field spectrum of Mathis et al. (1983, hereafter MMP83) the dust heating rate parameter U is

$$\begin{equation} U = \frac{u_\star }{8.6\times 10^{-13}\,{\rm erg}\,{\rm cm}^{-3}}, \end{equation} \tag{ 12 }$$

where u_⋆ is the starlight energy density. For the DL07 model, the characteristic grain temperature (of the grains dominating the emission at λ > 100 μm) is related to the heating rate parameter U as

$$\begin{equation} T_{\rm d,char}\approx 18\, U^{1/6} \,{\rm K}. \end{equation} \tag{ 13 }$$

This is only a representative temperature—dust grains of different sizes and compositions illuminated by a single radiation field have different steady-state temperatures, and very small grains undergo temperature fluctuations due to the quantized heating by stellar photons. Figure 9(b) shows T_{d, char} as a function of radius.

Zoom In Zoom Out Reset image size

The mean starlight heating rate 〈U〉 is shown in Figure 9(a). At S350 resolution, the center has 〈U〉 ≈ 32 (see Figure 9), corresponding to T_{d, char} ≈ 33 K. At R = 15 kpc we find U ≈ 0.3, or T_{d, char} ≈ 15 K, and at R = 20 kpc we find U ≈ 0.25, or T_{d, char} ≈ 14 K.

Near the center, our T_{d, char} ≈ 33 K is in close agreement with T_d ≈ 35 K estimated by Groves et al. (2012), who modeled the 100–500 μm SED using a modified blackbody, F_ν∝ν²B_ν(T_d) in ΔR = 230 pc annuli for 0 < R < 15 kpc. At R = 15 kpc, our T_{d, char} = 15 K is somewhat higher than the value 12.5 K found by Groves et al. (2012).

How does the starlight heating rate inferred from the dust IR emission compare with what is expected? At R > 5 kpc, the observed mean heating rate 〈U〉 in Figure 9(a) is presumably primarily due to starlight from disk stars. In the star-forming parts of the disk, the radiation field heating the dust is from a mixture of young and old stars, modified by dust attenuation; 〈U〉_disk includes the contribution of both young and old stars. Estimating the distribution of U values seen by the dust would be a challenging radiative transfer problem even if we knew the three-dimensional distributions of stars and dust. The observed 〈U〉 in Figure 10(a) shows an approximately exponential decline with increasing R out to 20 kpc, which can be approximated by

$$\begin{equation} \langle U\rangle _{\rm disk} \approx 1.0 \exp (-R/12\,{\rm kpc}). \end{equation} \tag{ 14 }$$

This is plotted in Figures 10; we will use it to estimate the contribution of disk starlight to 〈U〉 in the central regions.

Zoom In Zoom Out Reset image size

Near the center of M31, the radiation from the stellar bulge population becomes dominant. The distribution of bulge luminosity, and the resulting energy density of starlight from the bulge, is discussed by Groves et al. (2012), and elaborated further in Appendix B. From a three-dimensional model for the stellar bulge, the energy density of bulge starlight is estimated to be

$$\begin{equation} u_{\star,{\rm bulge}}^{(0)}(R) = 5.1\times 10^{-11} I(R/r_b) \,{\rm erg}\,{\rm cm}^{-3}, \end{equation} \tag{ 15 }$$

where I(x) is given by Equation (B4), r_b = 0.58 kpc is the core radius of the bulge (see Appendix B)¹⁹ and the superscript (0) indicates that it is a theoretical estimate, with dust extinction neglected. For the starlight spectrum of the bulge, we estimate

$$\begin{equation} U_{\rm bulge}^{(0)} = \frac{u_{\star,{\rm bulge}}^{(0)}}{1.7\times 10^{-12}\,{\rm erg}\,{\rm cm}^{-3}} = 30\,I(R/r_b). \end{equation} \tag{ 16 }$$

Equation (16) gives $U_{\rm bulge}^{(0)}(R=1\,{\rm kpc})=7.9$, whereas our estimated heating rate from the IR SED at R = 1 kpc is 〈U〉 ≈ 7.5 (see Figure 10), which presumably includes a contribution from disk stars, which is estimated from Equation (14) to be ∼0.92. Thus at R = 1 kpc we might estimate the heating contribution of the bulge stars to be ∼6.5, 80% of our theoretical estimate $U_{\rm bulge}^{(0)}$ at this radius. At smaller radii we find that $U_{\rm bulge}^{(0)}$ exceeds the inferred 〈U〉 by ∼30%.

There are a number of possible explanations for the discrepancy between the theoretical estimate $U_{\rm bulge}^{(0)}$ and 〈U〉 estimated from modeling the IR SED: (1) perhaps the bulge starlight has simply been overestimated by ∼30%. (2) The dust grains that we are using as "photometers" are likely to be mainly located in clouds distributed in a thin disk; internal extinction could lower the starlight heating rate below the optically thin estimate (Equation (16)). (3) Some of the dust is presumably above and below the disk plane; thus the observed emission from the center will include emission from dust that is actually at larger radii, projected onto the center. (4) To estimate 〈U〉 from the infrared observations, we have assumed the dust grains to have the properties of dust in the DL07 model. The dust near the center of M31 may well differ from the solar-neighborhood dust on which the DL07 is based. If the dust grains at R ≲ 2 kpc have a lower ratio of (optical absorption cross section)/(FIR absorption cross section) than the DL07 model grains, a given radiation field will heat them to a lower temperature than the DL07 dust. A 30% reduction in the ratio (optical absorption cross section)/(FIR absorption cross section) would be sufficient to remove the discrepancy. Such a reduction in the ratio of optical absorption to FIR emission cross section would occur if the grain radii in the center were larger by ∼30%.

To approximate our observed heating rate 〈U〉, we will adopt the spatial profile expected for the bulge starlight, but will scale down the heating rate by a factor 0.7. To this we add our empirical estimate (Equation (14)) for the heating due to the disk stars. Thus

$$\begin{eqnarray} \langle U\rangle &=& U_{\rm bulge} + \langle U\rangle _{\rm disk\,}, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} U_{\rm bulge} &=& 0.7 U_{\rm bulge}^{(0)} = 21\, I(R/r_b), r_b=0.58\,{\rm kpc}\,, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \langle U\rangle _{\rm disk} &\approx & 1.0\exp (-R/12\,{\rm kpc}). \end{eqnarray} \tag{ 19 }$$

This estimate of U is plotted in Figures 9 and 10. The agreement is good, except at R < 0.4 kpc (see Figure 10) where the observed dust heating rates fall somewhat below Equation (17). The deviation at R < 0.4 kpc is likely due to the limited spatial resolution of the observations. Even at S350 resolution, the FWHM of the PSF corresponds to ΔR = 0.4 kpc along the minor axis. In addition, as already noted above, out-of-plane dust, projected along the line-of-sight, can also act to lower the "observed" dust temperatures near the center.

Using modified blackbody fits, Groves et al. (2012) found dust temperatures in the center that are ∼15% lower than the T_{d, char} values found here. Their lower dust temperature estimates may be due to use of an earlier version of the data reduction pipeline and calibration factors, and possibly also to use of a modified blackbody rather than the multicomponent physical grain model used here.

To summarize, we find relatively good agreement between the observed dust emission spectrum and that which would be expected for heating by a combination of the bulge starlight and a disk heating component. The observed dust temperatures are only slightly below what would be expected. The inferred heating rate from the bulge stars is only ∼30% lower than predicted for a simple model of the bulge starlight. The 30% discrepancy in heating rate corresponds to only a 5% discrepancy in grain temperature. Given that there are a number of effects that could account for such a discrepancy, this agreement is gratifying.

The dust temperature T depends on the starlight heating rate parameter U and on the ratio 〈C_abs〉_⋆(a)/〈C_abs(a)〉_T, where 〈C_abs(a)〉_⋆ is the dust absorption cross section averaged over the spectrum of the illuminating starlight for a grain of radius a, and 〈C_abs(a)〉_T is the Planck-averaged absorption cross section for grain temperature T. The fact that the observed dust temperature is within 5% of the predicted dust temperature indicates that the actual values of 〈C_abs〉_⋆/〈C_abs〉_T are close to the values in the DL07 grain model. This builds confidence in our grain model, and hence in the dust mass estimates, which are proportional to ρa/〈C_abs(a)〉_T. Given the extreme environmental differences, it is remarkable that the dust properties near the center of M31 appear to be so similar to values inferred for dust in the solar neighborhood.

The intensity in the IRAC bands includes both direct starlight and emission from dust population. The observed intensities at λ ⩾ 3.6 μm are modeled as the sum of a stellar component (modeled as a 5000 K blackbody) plus a nonstellar component (F_ν)_ns contributed primarily by PAHs. In practice, the stellar component is determined by the IRAC3.6 and IRAC4.5 photometry; subtraction of the stellar component typically leaves a positive "nonstellar" residual in IRAC5.8 and IRAC8.0 that can be reproduced by varying the PAH abundance in the dust model.

The parameter q_PAH in the DL07 model is defined to be the fraction of the total dust mass contributed by PAHs containing fewer than 10³ C atoms. In the model fitting, q_PAH is essentially proportional to the ratio of the power in the 6.2 and 7.7 μm PAH emission features divided by the total power radiated by the dust. When only IRAC photometry is available for λ < 20 μm, the q_PAH parameter is, in practice, proportional to the nonstellar contribution to IRAC8.0 divided by the aggregate 70–350 μm luminosity. For a fixed starlight spectrum, F(8 μm) is approximately proportional to q_PAH and to the total infrared power, F_TIR:

$$\begin{equation} (\nu F_\nu)_{\rm IRAC8,ns}\approx A_\star q_{\rm PAH}F_{\rm TIR}\,, \end{equation} \tag{ 20 }$$

where (νF_ν)_{IRAC8, ns} is the nonstellar contribution to the IRAC 8.0 μm band, and the dimensionless coefficient A_⋆ depends on the spectrum of the illuminating starlight. Draine & Li (2007) show that A_MMP83 ≈ 4.72 for the MMP83 spectrum of the starlight in the solar neighborhood.

Varying the spectrum of the starlight can cause A_⋆ to change, for two reasons. (1) The PAH absorption cross section depends on wavelength differently from the absorption of overall grain mixture (∝F_TIR), hence the fraction of the starlight power that is absorbed by PAHs will depend on the spectrum. (2) The PAH emission spectrum is the result of single-photon heating, hence the fraction of energy absorbed by PAHs that is reradiated in the IRAC8.0 band depends on the illuminating spectrum. As an example of the dependence of A_⋆ on the illuminating spectrum, Draine (2011a) showed that A_20 kK ≈ 7.1 = 1.5A_MMP83 for a 20 kK blackbody cut off at 13.6 eV.

B. T. Draine & Li (in preparation) calculated the emission from the DL07 dust model for illumination by starlight with the spectrum of the M31 bulge population, finding A_bulge ≈ 1.95. For a mixed spectrum, where the overall dust heating rate is U = U_bulge + U_MMP83, the effective value is the dust luminosity-weighted mean:

$$\begin{eqnarray} A_\star &=& \frac{A_{\rm bulge}U_{\rm bulge}+A_{\rm MMP83}U_{\rm MMP83}}{U_{\rm bulge} + U_{\rm MMP83}} \nonumber\\ \ms &=& \frac{1.95 U_{\rm bulge}+4.72 U_{\rm MMP83}}{U_{\rm bulge}+U_{\rm MMP83}}. \end{eqnarray} \tag{ 21 }$$

To correct the estimate of q_PAH made assuming the MMP83 radiation field, we will take

$$\begin{equation} \left(q_{\rm PAH}\right)_{\rm corr} = \frac{4.72}{A_\star } \times \left(q_{\rm PAH}\right)_{\rm MMP83}, \end{equation} \tag{ 22 }$$

where (q_PAH)_MMP83 is the value of q_PAH estimated assuming the dust to be heated by starlight with the MMP83 spectrum, with U_bulge and U_MMP83 given by Equations (18) and (19).

A map of (q_PAH)_MMP83 at S350 resolution is shown in Figure 11(a). Figure 11(b) shows the radial profile of (q_PAH)_MMP83 and (q_PAH)_corr. Interestingly, q_PAH appears to peak in the 11 kpc ring, attaining a value (q_PAH)_corr ≈ 0.049 that is close to the value q_PAH ≈ 0.045 estimated for the diffuse ISM in the solar neighborhood.

Zoom In Zoom Out Reset image size

In Figure 11, (q_PAH)_MMP83 declines as one approaches the center, reaching a value (q_PAH)_MMP83 ≈ 0.02 in the central regions. However, this decline is largely an artifact of assuming that the spectrum of the illuminating starlight is independent of radius, which is incorrect—in the central regions the starlight is dominated by light from the bulge stars, which is much redder than the MMP83 spectrum.

Figure 11(b) shows that when we allow for the starlight being increasingly dominated by an old stellar population as we move to the center, (q_PAH)_corr shows only limited variation as we move from R = 11 kpc to the central kiloparsec. Evidently the balance between PAH formation and destruction in the ISM remains relatively constant from the central kiloparsec out to R = 20 kpc. There does appear to be a systematic radial decline in q_PAH for R > 11 kpc, but this again might be an artifact of a radial gradient in the spectrum of the starlight as one moves from the 11.2 kpc ring—where star formation is active—to outer regions where there appears to be little contemporary star formation.

At R ≳ 20 kpc q_PAH estimated from the ΔR = 677 pc annuli appears to rise. However, we suspect this to be an artifact of imperfect background subtraction in the 5.8 and 8.0 μm images. Although background subtraction works well in other bands out to ∼25 kpc, it appears less successful for the IRAC 5.8 μm and 8.0 μm bands. The difficulty with background subtraction may be due to systematic effects on the IRAC detectors, including an effect referred to as "banding," multiplexor "bleeding," and scattered light (Hora et al. 2004). The SAGE–SMC survey (Gordon et al. 2011) was able to minimize these problems by combining images taken with very different roll angles together with custom processing techniques, but we are simply using the M31 images from Barmby et al. (2006). It is also possible that the Galactic foreground has structure arising from variations in PAH abundance or ionization state on ∼05 scales, which would not be identified by the background subtraction procedures used here.

In the present work, we attempt to reproduce the observed SED in each pixel using the DL07 dust model and a parameterized distribution of starlight intensities. Above we have examined the values of the dust modeling parameters, such as the dust mass, q_PAH, and properties of the starlight intensity distribution. Here we compare the models with observations to see how well the model actually reproduces the data, and whether any systematic deviations are present that indicate systematic problems with the modeling.

Aside from minor effects associated with variation in the PAH abundance when q_PAH is allowed to vary, the composition (amorphous silicate and graphitic carbon) is assumed to be the same everywhere in the DL07 models used to fit the infrared emission. Further, the dust opacity is assumed to be independent of temperature. The DL07 dust model used here has an opacity with a fixed dependence on frequency at long wavelengths (λ ≳ 50 μm). Here we compare the model to the observed emission from M31 to look for residuals that might be indicative of differences between the model and the actual dust in M31.

We will focus on the behavior of the dust opacity in the 250–500 μm region. Figure 12 shows the ratio of model to observation at SPIRE250 and SPIRE500. The first impression is that the model is good: the ratio of model to observation is generally between 0.86 and 1.16, which seems good in view of noise in the observations, uncertainties in calibration, and the general uncertainties in the adopted dust opacities. However, systematic trends are evident: in the central regions of M31 (excluding the center itself) the model tends to be low at SPIRE250, and high at SPIRE500.

Zoom In Zoom Out Reset image size

Let β ≡ ln [κ(250 μm)/κ(500 μm)]/ln 2 be the effective power-law index of the dust opacity between 250 and 500 μm. For the DL07 model (see Table 2) this ratio is β_DL07 = 2.08. If the fitted dust temperatures were left unchanged, the 500/250 flux ratio could be brought into agreement with observations if β were changed to

$$\begin{eqnarray} &&\beta _{\rm obs} = \beta _{\rm DL07}\nonumber \\ &&\quad + \frac{\ln ([I_\nu (250\,{\mu {\rm m}})/I_\nu (500\,{\mu {\rm m}})]_{\rm obs}/ [I_\nu (250\,{\mu {\rm m}})/I_\nu (500\,{\mu {\rm m}})]_{\rm model})}{\ln (500/250)}. \nonumber \\ \end{eqnarray} \tag{ 23 }$$

Figure 13 shows β_obs as a function of R for R < 22 kpc, using only M160 resolution pixels with Σ_Ld > Σ_{Ld, min}. For 1 < R < 10 kpc, β_obs is larger than β_DL07, indicating that the opacity between 250 and 500 μm should fall more rapidly than ν^2.08.

Zoom In Zoom Out Reset image size

The dust opacity ratio κ(250 μm)/κ(500 μm) could vary with location because the dust composition is varying, or it could conceivably result from variations of the grain temperature, if the dust opacities are temperature dependent. The DL07 model assumes the dust opacities to be independent of temperature. However, some materials do exhibit temperature-dependent opacities in the laboratory (e.g., Mennella et al. 1998; Boudet et al. 2005; Coupeaud et al. 2011), and such behavior is expected in some models of amorphous solids (Meny et al. 2007). A number of studies have claimed that interstellar dust opacities are temperature dependent (Dupac et al. 2003; Paradis et al. 2010, 2011; Liang et al. 2012), although apparent β–T correlations can arise from both observational noise and line-of-sight temperature variations (Shetty et al. 2009a, 2009b; Kelly et al. 2012). The "two level system" (TLS) model (Meny et al. 2007; Paradis et al. 2011) predicts that increasing dust temperature T should lead to a lower value of β. However, in M31 it appears that the 1–10 kpc regions where the dust is warmer than the outer disk have a larger value of β than the value in the outer disk, and the center—where the dust is hottest—has essentially the same value of β as the relatively cold dust at R ≈ 20 kpc. The observed variations in β do not seem to be consistent with what would be expected for the TLS model, unless one allows for substantial radial variations in the TLS model parameters themselves. Rather than attributing the apparent changes in β to temperature, it would appear instead that the dust composition must be varying with radius in M31.

The HELGA II collaboration (Smith et al. 2012) reported a sharp change in the dust properties at R ≈ 3 kpc, based on the value of the dust emissivity index β obtained from their modified blackbody fits. Interior to 3 kpc, the inferred dust temperature decreased from T = 27.5 K at the center to T = 16.8 K at R = 3.0 kpc, with β simultaneously rising from 2 to 2.5; beyond R = 3 kpc, they found the temperature to be rising and β falling with increasing R, reaching T = 18 K and β = 1.5 at R = 15 kpc. The decrease in β at large R was invoked to account for an apparent "500 μm excess."

The HELGA II result for β, shown in Figure 13, is similar to the present study for R ≲ 7 kpc, but we find very different behavior for R ≳ 7 kpc. We see no evidence of a 500 μm excess at large R—quite the contrary, the DL07 model tends to slightly overpredict SPIRE500 in the outer regions of M31.

The HELGA II analysis was based on observations of M31 that were shallower than the data used here, and also used an earlier calibration of the three SPIRE bands. It is possible that the differences between the HELGA II results and those of the present study may be due in part to differences in signal/noise, calibration or data reduction pipeline. In the outer regions results are also sensitive to background subtraction.

If we were over- or under-subtracting the IR backgrounds, we would obtain very low or very high dust/gas ratios. The fact that we obtain sensibly behaved dust/H ratios out to R = 25 kpc (see Figures 6 and 7) is evidence that the automatic background subtraction algorithm used here (Aniano et al. 2012)—which is based only on the IR imaging, and makes no use of H i 21 cm or CO emission—is working well, at least out to R = 25 kpc.

Our overall conclusion is that the DL07 dust properties appear to provide a generally good match to the observations of M31, except for the need for a modest increase in the opacity index β in the 2–6 kpc region. It would be of interest to apply this same approach to other galaxies (e.g., M33) to see if similar variations in β_obs are found.

The principal conclusions of this work are as follows.

• 1.  
  Consistent with previous observations, we find that the dust mass surface density in M31 peaks in two rings, at R = 5.6 kpc and R = 11.2 kpc, with a third ring seen at R ≈ 15.1 kpc.
• 2.  
  We find a total dust mass M_d = (5.4 ± 1.1) × 10⁷ M_☉ within R = 25 kpc. 95% of this dust mass lies within R = 21 kpc, with the dust surface density peaking at R = 11.2 kpc.
• 3.  
  The dust/H mass ratio exhibits a smooth radial decline with increasing R, from ∼0.027 at the center to ∼0.0027 at R = 25 kpc.
• 4.  
  The dust/H mass ratio parallels measurements of O/H in H ii regions, consistent with a constant fraction of the refractory elements Mg, Si, and Fe being in dust. Based on our estimated dust/H mass ratio, we infer that the metallicity Z/Z_☉ varies from ∼3 at R = 0 to ∼0.3 at R = 25 kpc—see Figure 7(b) and Equation (10).
• 5.  
  The starlight heating rate parameter 〈U〉 shows a nearly monotonic decline with galactocentric radius, from 〈U〉 ≈ 50 at the center (at S350 resolution) to 〈U〉 ≈ 0.2 at R ≈ 20 kpc.
• 6.  
  We confirm the finding of Groves et al. (2012) that the dust heating in the central 2 kpc is dominated by light from the stellar bulge. We find that the starlight heating rates inferred from the observed IR emission are consistent with the heating rates expected from the bulge starlight. It is remarkable that the dust properties near the center of M31 appear to be similar to the properties of dust in the solar neighborhood.
• 7.  
  After taking into account variation in the spectrum of the starlight, we find the PAH abundance q_PAH to be approximately constant from the center of M31 out to ∼20 kpc. There is some indication of decline with R for R > 11 kpc. The global value of q_PAH = 0.039 for R < 17 kpc.
• 8.  
  While the DL07 dust model generally provides a good fit to the observed SED, there are some systematic deviations. The dust at R ≈ 1–6 kpc appears to have an opacity spectral index 2.2 < β < 2.33 in the 250–500 μm wavelength range, whereas the dust at R > 7 kpc has β ≈ 2.08, consistent with the dust in the DL07 model. We are in approximate agreement with the radial variation of β found by HELGA II (Smith et al. 2012) in the central ∼7 kpc, but do not confirm their finding of low β values for R > 8 kpc.
• 9.  
  At large radii R ≳ 10 kpc the DL07 model, with fixed opacity, is consistent with the observed photometry, and returns dust masses that are consistent with the observations of the gas and metallicity.

We thank Edvige Corbelli, Stephen Eales, Jeremy Goodman, Matthew Smith, and the anonymous referee for helpful comments. This work was supported in part by NSF grant AST 1008570. G.A. acknowledges support from European Research Council grant ERC-267934.

MIPS and PACS have two wavelengths in common: 70 μm and 160 μm. While the filter response functions are not identical, they are similar enough that we would expect MIPS and PACS to measure very similar flux densities for smooth SEDs. However, comparisons of MIPS and PACS imaging often shows discrepancies that are much larger than expected—see, e.g., the cases of NGC 628 and NGC 6946 (Aniano et al. 2012).

Figure 14 shows the ratio of PACS/MIPS photometry of M31 in the 70 and 160 μm bands, after first convolving each image to the common resolution of the MIPS160 PSF.

Zoom In Zoom Out Reset image size

The PACS70/MIPS70 comparison is particularly striking: while there are a few regions near the center with PACS70/MIPS70 ∼ 0.7, many regions have PACS70/MIPS70 >2. The fact that this occurs over extended regions makes it clear that this is not simply a consequence of random noise. The fact that it is non-uniform over the image makes it clear that it is not simply a result of calibration. Background subtraction is of critical importance, but we have employed the same background estimation procedures (Aniano et al. 2012) for both MIPS and PACS. It is now thought that MIPS70 suffers from sublinear behavior in high surface brightness regions, but such high surface brightnesses are found only at the center of M31, whereas we see high values of PACS70/MIPS70 occurring across the disk. The reason for the photometric discrepancy is unclear, and we therefore opt to use both MIPS70 and PACS70 data in our model-fitting.

The PACS160/MIPS160 comparison is much more satisfactory than the PACS70/MIPS70 comparison, but there are still many regions—particularly along the 11 kpc ring—where PACS160/MIPS160 >1.4. This difference is much larger than the claimed uncertainties in the observed intensities at these locations. MIPS160 is now thought to have sublinear response for I_ν ≳ 50 MJy sr⁻¹ (Paladini et al. 2012), but most of M31 is below this value. This again emphasizes the value of being able to include the MIPS160 data in the analysis, possible only if one degrades all the other imaging to the MIPS160 PSF (FWHM 39'').

Groves et al. (2012) have discussed the contribution of M31's bulge stars to the heating of dust. For the adopted bulge stellar luminosity density profile (Geehan et al. 2006)

$$\begin{equation} \rho _\star (R) = \frac{L_{\rm bulge}}{2\pi r_b^3} \frac{1}{(R/r_b)(1+R/r_b)^3}, \end{equation} \tag{ B1 }$$

the starlight energy density due to bulge stars (neglecting possible absorption or scattering by dust) is

$$\begin{eqnarray} u_{\star,\rm bulge}(R) &=& \frac{L_{\rm bulge}}{4\pi r_b^2 c} I(R/r_b),\\ \ms I(x) &\equiv & \frac{1}{x}\int _0^\infty \frac{\ln \left(\frac{x+y}{|x-y|}\right)}{(1+y)^3}dy\,\\ \ms &=& \frac{1}{x^2-1} + \frac{\ln (1+x)}{2x(1+x)^2} - \frac{2\ln x}{(x^2-1)^2}. \end{eqnarray} \tag{ B2 }$$

I(x) is logarithmically divergent for x → 0. Near x = 1,

$$\begin{equation} I(x=1+\epsilon) = \frac{1}{2} + \frac{\ln 2}{8} -0.778(x-1) + O((x-1)^2). \end{equation} \tag{ B5 }$$

The projected luminosity interior to radius R is

$$\begin{eqnarray} L({<}R) &=& \int _0^R \rho _\star (r) 4\pi r^2 dr + \int _R^\infty \rho _\star (r) \nonumber \\ &&\times [1-\sqrt{1-(R/r)^2}] 4\pi r^2 dr\\ \ms &=& L_{\rm bulge} \left[1 - 2\alpha \int _0^{\pi /2}\frac{\cos ^2\theta d\theta }{(1+\alpha \sin \theta)^3} \right] \alpha \equiv \frac{r_b}{R}\qquad \\ \ms &=& 0.4711 L_{\rm bulge}\quad {\rm for}\ \alpha =0.61. \end{eqnarray} \tag{ B6 }$$

We obtain L_bulge by integrating the average of the intrinsic and reddened spectrum obtained by Groves et al. (2012) for L(< 1 kpc) over the wavelength interval [0.0912 μm, 4 μm]. Groves et al. (2012) assumed D = 780 kpc and r_b = 0.61 kpc. Corrected to D = 744 kpc, we obtain L_bulge(0.0912–4 μm) = 1.60 × 10¹⁰ L_☉, with r_b = 0.58 kpc. Thus, the 0.0912 μm < λ < 4 μm energy density due to the bulge stars is

$$\begin{equation} u_{\star,{\rm bulge}}(R) = 5.07\times 10^{-11} I(R/r_b) \,{\rm erg}\,{\rm cm}^{-3}. \end{equation} \tag{ B9 }$$

This energy density of bulge starlight gives a theoretical estimate for the dust heating rate, normalized to the heating rate for the MMP83 radiation field,

$$\begin{equation} U_{\rm bulge}^{(0)}(R) = 30 I(R/r_b). \end{equation} \tag{ B10 }$$

For R = 1 kpc this gives $U_{\rm bulge}^{(0)}(1\,{\rm kpc}) = 7.9$. This is our zeroth order estimate for the dust heating rate in the central few kiloparsecs of M31. In Section 7 we show that this theoretical estimate for the dust heating rate agrees well with the dust heating rate inferred from the IR SED.