TThe Astrophysical Journal , , 172 (2014 Jan. 10) Andromeda’s Dust
B. T. Draine , , G. Aniano , , Oliver Krause , Brent Groves , Karin Sandstrom , Robert Braun ,Adam Leroy , Ulrich Klaas , Hendrik Linz , Hans-Walter Rix , Eva Schinnerer , AnikaSchmiedeke , and Fabian Walter ABSTRACT
Spitzer Space Telescope and
Herschel Space Observatory imaging of M31 is used,with a physical dust model, to construct maps of dust surface density, dust-to-gas ratio,starlight heating intensity, and polycyclic aromatic hydrocarbon (PAH) abundance, outto R ≈
25 kpc. The global dust mass is M d = 5 . × M (cid:12) , the global dust/H massratio is M d /M H = 0 . (cid:104) q PAH (cid:105) = 0 . R = 5 . R = 11 . R ≈ . M d /M H ≈ . ∼ . R ≈
25 kpc. From the dust/gas ratio, we estimate theinterstellar mediu (ISM) metallicity to vary by a factor ∼
10, from
Z/Z (cid:12) ≈ R = 0to ∼ . R = 25 kpc. The dust heating rate parameter (cid:104) U (cid:105) peaks at the center, with (cid:104) U (cid:105) ≈
35, declining to (cid:104) U (cid:105) ≈ .
25 at R = 20 kpc. Within the central kiloparsec, thestarlight heating intensity inferred from the dust modeling is close to what is estimatedfrom the stars in the bulge. The PAH abundance reaches a peak q PAH ≈ .
045 at R ≈ . q PAH for the dust in the central kiloparsec is similar to the overall value of q PAH in thedisk. The silicate–graphite–PAH dust model used here is generally able to reproduce theobserved dust spectral energy distribution across M31, but overpredicts 500 µ m emissionat R ≈ R = 2–6 kpc, the dust opacity varies more steeplywith frequency (with β ≈ . µ m) than in the model. Subject headings: dust, extinction – infrared: galaxies – infrared: ISM Princeton University Observatory, Peyton Hall, Princeton, NJ 08544-1001, USA; [email protected] Osservatorio Astrofisico Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy Institut d’Astrophysique Spatiale, F-91405 Orsay, France; [email protected] Max-Planck-Institut fur Astronomie, Konigstuhl 17, D-69117 Heidelberg, Germany CSIRO – Astronomy and Space Science, PO Box 76, Epping, NWS 1710, Australia National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA a r X i v : . [ a s t r o - ph . C O ] J a n
1. Introduction
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 externalstar-forming galaxy with spatial resolution that is surpassed only for the Magellanic Clouds. Thestructure 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 (cid:48) = 20 . I (Braun et al. 2009; Chemin et al. 2009) andinfrared emission (Haas et al. 1998; Barmby et al. 2006; Gordon et al. 2006; Fritz et al. 2012; Smithet al. 2012) show structure that appears as much ring-like as spiral in character. The centers ofthe rings are often offset significantly from the dynamical center of M31. The off-center ring-likestructure 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),
Infrared Space Observatory ( ISO ) (Haas et al.1998), and
Spitzer (Gordon et al. 2006). The total infrared luminosity was well-measured, but thelimited angular resolution of
IRAS
60 and 100 µ m (105 (cid:48)(cid:48) FWHM),
ISO µ m (110 (cid:48)(cid:48) FWHM), and
Spitzer µ m (39 (cid:48)(cid:48) FWHM) allowed only a relatively coarse image of the dust distribution.The present study takes advantage of the high sensitivity and angular resolution of
HerschelSpace Observatory (Pilbratt et al. 2010) to characterize the dust in M31 on angular scales as smallas ∼ (cid:48)(cid:48) (= 90 pc @ 744 kpc), and at wavelengths as long as 500 µ m, thereby capturing the emissionfrom 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 dustdistribution in M31, using single-temperature modified blackbody fits to the 70 − µ m emissionfrom each pixel. The present study differs from HELGA in two ways. First, we use an independentset of Herschel observations (Groves et al. 2012, O. Krause et al. 2013, in preparation), somewhatdeeper than those obtained by HELGA. Secondly, we use a physical dust model (Draine & Li 2007)to model the spectral energy distribution (SED) from 6 − µ m, and to estimate the dust masssurface 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 thegalaxies in the KINGFISH sample (Kennicutt et al. 2011).The organization of the paper is as follows: the observational data are described in Section All radii, luminosities and masses in this paper have been corrected to D = 744 kpc. For example, in model P1 of Corbelli et al. (2010) the center of the ring at 10–12 kpc is offset from the dynamicalcenter by 0 .
2, and the dust-model fitting methods are outlined in Section 3. Results are presented in Sections4–8, where we estimate the total dust luminosity and mass, the dust/gas ratio, the metallicity asa 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 summarizedin Section 10.Appendix A examines the inconsistency between PACS and MIPS photometry of M31. Thestarlight contribution from the M31 bulge stars is calculated in Appendix B.
2. Observational Data
M31 has been mapped using the IRAC (Fazio et al. 2004) and MIPS (Rieke et al. 2004) camerason
Spitzer Space Telescope (Werner et al. 2004). More recently, maps have been obtained by thePACS (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 Gordonet al. (2006). The PACS and SPIRE imaging was done in parallel mode at medium scan speed(20 (cid:48)(cid:48) s − ) for a total time of ∼
24 hr (O. Krause et al. 2013, in preparation), and are the samedata as used by Groves et al. (2012), except for small changes in calibration. We use the mostrecent calibrations of PACS and SPIRE. We refer to the different images by the camera nameand 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 α J2000 = 10 . ◦ , δ J2000 = 41 . ◦ (Crane et al. 1992), orGalactic coordinates (cid:96) = 121 ◦ , b = − ◦ . Because M31 is only 22 ◦ from the Galactic plane, removalof the Galactic foreground “cirrus” emission is challenging, particularly in view of the large angularextent of M31. We are helped by the high inclination i ≈ ◦ of M31, which raises the surfacebrightness, and increases the contrast with foreground and background emission.Subtraction of foreground and background emission has been carried out following methodsdescribed 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 IRAC images in bands 1–4 were multiplied by extended source calibration factors 0.91, 0.94, 0.66, 0.74 (Reachet al. 2005). MIPS images were generated by the Mips enhancer v3.10 pipeline on 2007 Jul 3. The PACS and SPIRE images were processed by HIPE v9, and the Level 1 HIPE images were then processedby Scanamorphos v18.0 (Roussel 2013). We used the calibration files in HIPE v9 (version 42 for PACS, and version10.0 for SPIRE). Intensities in the SPIRE bands were obtained by dividing the HIPE v9 flux density per beam byeffective beam solid angles Ω = (1 . , . , . × − sr for SPIRE250, 350, and 500, as recommended by Griffinet al. (2013). Fig. 1.—
Background-subtracted 350 µ m images of M31. Left: native resolution of SPIRE350. Right: SPIRE350image convolved to the MIPS160 PSF. It is seen that most of the structure visible at S350 resolution remains at M160resolution, with improved signal/noise in low surface brightness regions. The contour delineates the “galaxy mask”:the region within which the signal-to-noise ratio (S/N) is high enough that dust modeling can be done pixel-by-pixel(see text). The ellipse shown is a tilted circle of radius R = 5500 (cid:48)(cid:48) = 19 . i = 77 . ◦ and P . A . = 37 . ◦ (Corbelli et al. 2010). 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 ofthe large angular extent of M31. For IRAC5.8 and IRAC8.0, it was found necessary to use morecomplex curved surfaces, rather than simple tilted planes, to model the background. Even so, thebackground 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 ∼ ◦ scales (see, e.g., the IRAS µ m or LAB21 cm map of the area around M31), and the simple tilted plane model (or low-order curved surfacein the case of IRAC5.8 and IRAC8.0) will not remove all of the foreground cirrus. This will beproblematic in the low surface brightness outer regions of M31. However, after the foreground- andbackground-subtracted FIR images have been fitted by dust models, comparison of the resultingdust map with maps of the H I at the radial velocity of M31 will allow us to assess the effectsof imperfect foreground subtraction. Except for IRAC5.8 and IRAC8.0, the present tilted planeforeground model appears adequate for R (cid:46)
25 kpc (see Section 5).For modeling the dust in resolved systems, it is essential that multiband imaging be convolvedto a common point spread function (PSF). The present study is carried out at two angular res-olutions: that of the SPIRE350 camera (FWHM=24 . (cid:48)(cid:48) = 90 pc @ 744 kpc, henceforth referred toas S350 resolution ) and that of the MIPS160 camera (FWHM=39 (cid:48)(cid:48) = 141 pc @ 744 kpc, henceforthreferred to as
M160 resolution ). Image convolutions were carried out using kernels obtained asdescribed by Aniano et al. (2011). For studies at S350 resolution, we are unable to use SPIRE500or MIPS160 imaging, because the PSF of those cameras is too broad. The convolved images atS350 resolution were sampled with 10 (cid:48)(cid:48) × (cid:48)(cid:48) pixels. Images convolved to M160 resolution weresampled with 18 (cid:48)(cid:48) × (cid:48)(cid:48) pixels.Figure 1 shows SPIRE350 images of M31; the left image is at S350 resolution; the image onthe right is convolved to M160 resolution. The lower resolution of M160 is evident when the imagesare compared, but most details of the image that are seen at S350 resolution are also identifiable atM160 resolution. Analysis of M31 at M160 resolution – benefiting from the improved signal/noise ofthe larger M160 pixels as well as being able to use both SPIRE500 and MIPS160 data to constrainthe models – will, therefore, not lose important structural features.M31 is highly inclined. Using DRAO H I i = 74 . ◦ and P . A . = 37 . ◦ for the disk at R = 6 −
27 kpc. From WSRT H I i = 77 . ◦ and position angle P . A . = 37 . ◦ for R = 10 −
13 kpc. Corbelli et al. (2010) showed that the H I kinematics can be described usingcircular rings with centers offset from the dynamical center of M31, and with inclinations between75 ◦ and 79 ◦ for 10 < R <
25 kpc).We will assume the disk to have a single inclination i = 77 . ◦ and position angle P . A . = 37 . ◦ .Figure 1 shows ellipses with major radius 5500 (cid:48)(cid:48) corresponding to R = 19 . D = 744 kpc,with i = 77 . ◦ and P . A . = 37 . ◦ . Radial trends will be studied by averaging in annuli centered onthe location of the presumed supermassive black hole, α J2000 = 10 . ◦ , δ J2000 = 41 . ◦ (Craneet al. 1992).The S350 PSF provides excellent spatial resolution (FWHM = 90 pc along the major axis). Wewill compare results obtained at S350 resolution with results obtained using the M160 PSF. Whilethe dust models at M160 resolution smooth out dust structures smaller than ∼
140 pc along themajor axis (or ∼
700 pc along the minor axis), working at M160 resolution permits (1) extendingthe wavelength coverage to 500 µ m (by including SPIRE500 photometry), and (2) using MIPS160photometry, which is somewhat deeper than the PACS160 imaging. In addition, the larger pixelsprovide improved signal/noise in low surface brightness regions. At S350 resolution, the models areconstrained 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 reliabilityof the method.
3. Modeling the Dust and Starlight Heating3.1. Dust Model
We model the dust in M31 using the dust models of Draine & Li (2007, hereafter DL07). Thedust is assumed to be a mixture of carbonaceous grains and amorphous silicate grains with a sizedistribution that is consistent with studies of the wavelength-dependent extinction and infraredemission produced by dust in the diffuse interstellar medium (ISM) in the solar neighborhood(Weingartner & Draine 2001, hereafter WD01). We treat the PAH abundance as variable, propor-tional to the parameter q PAH , where q PAH is defined to be the fraction of the total dust mass thatis contributed by PAHs containing fewer than 10 C atoms.The optical properties of the dust mixture are based on an adopted grain size distribution,dielectric functions for the graphite and silicate materials, and adopted absorption cross sections perC atom for the PAHs. Heat capacities are also required to calculate the temperature distributionfunction for the smaller carbonaceous and silicate grains (Draine & Li 2001).By fitting the dust model to the extinction curve per H in the solar neighborhood, the grainvolume per H is constrained. We assume the size distribution for the local Milky Way dust with R V = 3 . V ( cm / H) M d /M H M d /M H DL07 DL07 Renormalized DL07Amorphous silicate 3 . × − a b Carbonaceous 2 . × − c d Total 5 . × − a,c b,d Observed toward ζ Oph 0.0091 ea For ρ sil = 3 . − . b For ρ sil = 3 . − . c For ρ carb . = 2 . − . d For ρ carb . = 1 . − . e From Table 23.1 of Draine (2011b).
The dust volume/H in the DL07 dust model is determined by fitting the observed extinctionin the solar neighborhood; the dust mass/H can then be calculated if solid densities are assumedTable 2: Some Properties of Renormalized DL07 Dust ModelProperty Silicates Carbonaceous
Total A V / Σ M d (mag/( M (cid:12) pc − )) 5.403 13.39 7.394 κ (60 µ m) ( cm g − ) 102.5 106.3 103.4 κ (70 µ m) ( cm g − ) 73.62 72.24 73.28 κ (100 µ m) ( cm g − ) 34.59 29.02 33.20 κ (160 µ m) ( cm g − ) 13.11 10.69 12.51 κ (250 µ m) ( cm g − ) 5.263 3.981 4.941 κ (350 µ m) ( cm g − ) 2.448 2.060 2.351 κ (500 µ m) ( cm g − ) 1.198 1.062 1.164 κ (850 µ m) ( cm g − ) 0.4865 0.4202 0.4699 κ (1 .
28 mm) ( cm g − ) 0.2224 0.1866 0.2135 κ (2 .
10 mm) ( cm g − ) 0.1130 0.0914 0.1076 κ = absorption cross section per unit dust massfor the dust materials. For the silicates, DL07 assumed a density ρ = 3 . − , and for thecarbonaceous grains the density contributed by carbon alone was taken to be ρ C = 2 . − .With these adopted densities, the DL07 dust model that reproduces the extinction per H in thesolar neighborhood has M d /M H = 0 . M d /M H = 0 . ζ Oph, the observed depletionsimply M d /M H = 0 . ± . In the present paper we will assume solar neighborhood dust to have M d /M H = 0 . ∼ ρ sil reduced from 3 . − to 3 . − ,and the density of carbonaceous grain material reduced from 2 . − to 1 . − . Table 2gives A V / Σ M d for the renormalized DL07 model, as well as opacities at selected FIR wavelengths.With this renormalization, all dust masses and dust/H ratios in Aniano et al. (2012) should This value involves an assumption about the oscillator strength for C II ]2325˚A and a second assumption aboutthe depletion of oxygen – see Draine (2011b). be reduced by a factor f M = 0 . / . . The present study seeks to estimate dust masses in M31 by modeling the observed FIR andsubmillimeter emission. The dust is assumed to be heated by starlight, and the modeling has thefreedom to adjust the starlight intensities such that the grains are heated to temperatures suchthat their emission spectrum is consistent with the observed shape of the SED. Then, with thedust temperatures set by the starlight heating rates, the dust mass is proportional to the observedemission. Note that there is no single temperature characterizing the dust – even in a singleradiation field, grains of different size and composition have different temperatures, and the verysmall grains undergo temperature fluctuations due to single-photon heating. Additionally, a single“pixel” – which may be several hundred parsecs in transverse dimension – may include subregionswith different starlight intensities.For a given starlight intensity U , grain size a and composition, we solve for the temperatureprobability 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( p ν ) U,a, comp = 4 π (cid:90) dT ( dP/dT ) U,a, comp ( C abs ,ν ) a, comp B ν ( T ) , (1)where ( C abs ,ν ) a, comp is the absorption cross section at frequency ν . The emission per unit mass fora dust mixture exposed to starlight U is (cid:18) dL ν dM d (cid:19) U = (cid:80) comp (cid:82) da ( dn/da ) comp ( p ν ) U,a, comp (cid:80) comp (cid:82) da ( dn/da ) comp (4 π/ a ρ comp , (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: L ν,j = (cid:90) dU dM d ,j dU (cid:18) dL ν dM d (cid:19) U . (3)For a region j of solid angle Ω j , with dust mass M d ,j = Σ M d ,j Ω j D , the dust mass dM d ,j exposed tostarlight intensities in [ U, U + dU ] is assumed to be given by the simple parameterization proposedby DL07: dM d ,j dU = Σ M d ,j Ω j D (cid:34) (1 − γ j ) δ ( U − U min ,j ) + γ j ( α j − U − α j U − α j min ,j − U − α j max (cid:35) , (4)where Σ M d ,j is the total dust mass surface density in region j , δ is the Dirac δ -function, and α j > − γ 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 d /dU ∝ U − α j . For NGC 628 and NGC 6946, Aniano et al. (2012) found that the parameter U max couldbe 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: { Σ M d ,j , q PAH ,j , U min ,j , α j , γ j } . We require Σ M d ,j ≥
0, 0 ≤ γ j ≤
1, and 1 < α j <
3. We use oneadditional parameter to characterize the contribution of direct starlight (taken to have the spectrumof a 5000K 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 MIPS70and PACS70 cover essentially the same wavelengths, and likewise for MIPS160 and PACS160, theeffective 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 effectsof changes in the starlight spectrum when the bulge stars make a significant contribution to thedust heating. The model SED depends on the modeled distribution of dust temperatures and on the wavelength-dependence of the dust opacity κ ν . The wavelength-dependence of κ ν is tested by whether themodel can reproduce the observed shape of the SED, but if κ ν is in error by some constant factor A , the derived dust masses will be off by a factor 1 /A .The DL07 dust model uses adopted dielectric functions for the amorphous silicate and car-bonaceous grains which, with the assumption of spherical shape, allows absorption cross sectionsto be calculated. For solar-neighborhood abundances, the absorption cross section per H is givenin Table 3 for λ = 100 µ m, 250 µ m, and 500 µ m.Planck Collaboration et al. (2011) examined the FIR emission in the diffuse ISM, finding thatthe SED for low-velocity H I was consistent with κ ∝ ν β , with β ≈ .
8, and τ (250 µ m) /N H =(1 . ± . × − cm / H. As seen in Table 3, the
Planck estimate for τ /N H agrees well with theDL07 value at 100 µ m, but at 250 and 500 µ m the DL07 opacities are smaller, by a factor 1.3 at250 µ m, and a factor 1.5 at λ = 500 µ m. However, we will see below (Section 9) that the shape ofthe M31 SED seems more consistent with the β ≈ . β ≈ . Planck opacities are correct, then modeling the FIR emission using the DL07 modelwill tend to overestimate the dust mass. The present study does not use data longward of 500 µ m, Maps of the best-fit values of Σ M d ,j , q PAH ,j , U min ,j , α j , and γ j ∼ draine/m31dust ∼ µ m), but the possibility of a systematic error in the dust mass estimate should be kept in mind.Future studies employing the full Planck intensity data can be expected to shed further light ondust opacities (Planck Collaboration et al. 2013a,b,c).Table 3: Far-infrared Absorption/H τ ( λ ) /N H (10 − cm / H) λ ( µ m) DL07 Model Planck Observations a
100 5.0 5 . ± . b
250 0.75 1 . ± . . ± . ba Planck Collaboration et al. (2011). b Obtained from τ (250 µ m) /N H assuming β = 1 . Determination of the dust mass requires that the shape of the dust SED be well-determined, sothat the dust temperatures can be constrained. With the signal-to-noise properties of the presentobservations, we do not attempt to estimate the dust mass in a single 10 (cid:48)(cid:48) × (cid:48)(cid:48) pixel unless thedust luminosity per unit area (on the sky plane) Σ L d > Σ L d , min = 1 . × L (cid:12) kpc − (IR intensity I > . × − erg s − cm − sr − ). The irregular contour in Figure 1 and the other images bounds thecontiguous region with Σ L d > Σ L d , min (there are additional pixels with Σ L d > Σ L d , min outside thiscontour). For Σ L d (cid:46) Σ L d , min the single-pixel dust mass estimates are quite uncertain, especiallyat S350 resolution where the pixels are smaller, and SPIRE500 and MIPS160 cannot be used. AtM160 resolution the larger (18 (cid:48)(cid:48) × (cid:48)(cid:48) ) pixels and inclusion of SPIRE500 and MIPS160 help, buteven at M160 resolution the single-pixel dust mass estimates near Σ L d , min = 1 . × L (cid:12) kpc − havesubstantial 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 (cid:48)(cid:48) along the minor axis, and ∆ R = 46 . (cid:48)(cid:48) (169 pc) along the major axis. While these annuli are not resolved along the minor axis (theMIPS160 FWHM=39 (cid:48)(cid:48) ), they are resolved along the major axis. The boundary contour is established by modeling at M160 resolution (18 (cid:48)(cid:48) pixels), using all 13 bands to estimateΣ L d . When modeling at higher resolution (e.g., S350, with 10 (cid:48)(cid:48) pixels) we continue to use the contours established atM160 resolution to define the “galaxy mask”. If the center of a pixel falls within an annulus, we include the entire pixel in the sums for that annulus – we donot divide pixels that overlap annular boundaries. (a) Dust luminosity per area on the sky plane Σ L d at M160 resolution. (b) Σ L d as a function of radius R .Triangles: annular average of single-pixel modeling at S350 and M160 resolution. Circles: Σ L d based on modelingof SED of ∆ R = 677 pc annuli, using all photometric bands. Horizontal dot-dash line: surface brightness below whichwe do not attempt to model individual pixels. Dashed line: exponential fit (5) for R (cid:38)
10 kpc, with 3 . R that have Σ L d > Σ L d , min . The drop below 100% for R > . L d area in the SW. In order to study the dust at
R >
17 kpc, where a large fraction of the pixels have Σ L d < Σ L d , min (see Figure 2(c)), we use the background-subtracted images to extract the flux in each band fromlarger concentric annuli with widths 40 (cid:48)(cid:48) along the minor axis, and ∆ R = 188 (cid:48)(cid:48) (677 pc) alongthe major axis. These ∆ R = 677 pc annuli are fully resolved, even by MIPS160. The dust massand starlight intensity distribution in each annulus are then estimated by fitting the SED of theannulus. By integrating over the many pixels in each annulus, the random noise is substantiallysuppressed. We will see that reliable estimates of the mean dust mass surface density are possibleout to R ≈
25 kpc.
4. Dust Luminosity
As described in Aniano et al. (2012), the background-subtracted multiwavelength images arefirst used to estimate Σ L d , 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, andthen integrating over the model SED. While the model dust mass in individual pixels is unreliablefor Σ L d (cid:46) Σ L d , min , the pixel luminosity estimate is reliable down to considerably lower surfacebrightness. Figure 2(a) shows the dust luminosity surface density Σ L d at M160 resolution. From2Table 4: Global Quantities for M31 a Quantity
R <
17 kpc 17 −
20 kpc 20 −
25 kpc Total L d ( L (cid:12) ) 3 . × . × . × . × M d ( M (cid:12) ) 4 . × . × . × . × M H ( M (cid:12) ) 4 . × . × . × . × (cid:104) q PAH (cid:105)
MMP83 (cid:104) q PAH (cid:105) corr M d /M H a D = 744 kpc, inclination i = 77 . o the image it can be seen by eye that the dust luminosity is detectable down to Σ L d ≈ . L (cid:12) kpc − .The mean luminosity per area as a function of radius is found by averaging the single-pixelluminosities over the ∆ R = 177 pc annuli. Figure 2(b) shows Σ L d cos i , the dust luminosity surfacedensity on the disk plane, averaged over the ∆ R = 177 pc annuli, and plotted against the annularradius 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 resultsof modeling the SED of ∆ R = 677 pc annuli. The smooth decline in Σ L d with increasing R out to R ≈
23 kpc is an indication that background subtraction has been reasonably successful at leastdown to surface brightnesses Σ L d ≈ × L (cid:12) kpc − (Σ L d cos i ≈ × L (cid:12) kpc − ).Figure 2(c) shows the fraction f of the pixels at each R that have Σ L d > Σ L d , min . The fraction f = 1 for R <
15 kpc, but at R ≈ . R ≈
24 kpc the fractionhas fallen to f ≈ .
1; some fraction of these pixels may be raised above Σ L d , min by unresolvedbackground galaxies.The dust luminosity surface density Σ L d derived from the SEDs of ∆ R = 677 pc annuli, shownin Figure 2(b), is in good agreement with Σ L d obtained from pixel-based modeling. From theannular SEDs, the dust luminosity surface density for 17 kpc < R <
25 kpc isΣ L d cos i ≈ . × exp[ − ( R −
17 kpc) / . L (cid:12) kpc − . (5)For purposes of estimating the integrated dust properties from the pixel-based modeling, we sumover 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 . × L (cid:12) .Because νL ν peaks near 100 µ m, the global dust luminosity of M31 was already reasonablywell-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 ≈ . × L (cid:12) . Using MIPS All luminosities have been corrected to D = 744 kpc (a) Map of dust mass surface density Σ M d at M160 resolution. The irregular contour delineates thecontiguous region with Σ L d > Σ L d , min ; outside of this contour, the dust map is unreliable. (b) Radial profile of Σ M d projected onto the M31 disk plane, for cos i = 0 . M d from single-pixel modelingat S350 and M160 resolution. For each R , Σ M d is the average for pixels with Σ L d > Σ L d , min . Circles: Σ M d obtainedby modeling SED of ∆ R = 677 pc annuli, using all bands. photometry, Gordon et al. (2006) found L d ≈ . × L (cid:12) , consistent, within the uncertainties,with the present value of L d = 4 . × L (cid:12) .There is considerable structure in the dust luminosity image in Figure 2(a), with close similar-ities to the SPIRE 350 µ m image in Figure 1. While the structure is not purely “circular” (on theplane of the galaxy), the dust luminosity density extracted in (inclined) circular annuli centered onthe dynamical center, shown in Figure 2(b), exhibits a strong central peak (due to small amountsof relatively warm dust), two clear rings (at R = 5 . R = 11 . R ≈ . . IRAS imaging (Habinget al. 1984), and all three rings were noted in the
ISO µ m imaging by Haas et al. (1998).
5. Dust Mass, Dust-to-Gas Ratio, and Metallicity5.1. Dust Mass
Before the launch of
Herschel Space Observatory , the mass of dust in M31 had been estimatedby 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, and4Fig. 4.— (a) Dust luminosity interior to radius R , showing radii containing 1%, 10%, 50%, 90%, and 95% of thetotal dust luminosity within 25 kpc. (b) Dust mass interior to radius R , showing radii containing 1%, 10%, 50%,90%, and 95% of the total dust mass interior to 25 kpc. For R <
17 kpc L d ( r < R ) and M d ( r < R ) are based onsingle-pixel modeling at M160 resolution; for 17 kpc < R <
25 kpc the results are based on modeling the SEDs of∆ R = 677 pc annuli. to map the dust out to ∼
20 kpc. Figure 3(a) shows the dust mass surface density Σ M d on the planeof the sky, obtained by modeling at M160 resolution.In addition to constructing maps showing the dust properties in every pixel, we also averagethe 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 dustmasses in pixels with Σ L d > Σ L d , min , dividing by the total area of the annulus. For R < . L d > Σ L d , min , but this procedure will underestimate the actual mean surface densityfor 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 modelingat 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) agreeto within ∼
10% with dust mass estimates obtained at M160 resolution, using all cameras.Dust mass surface densities Σ M d 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 growingfraction of pixels has Σ L d < Σ L d , 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 in5Figure 3b:Σ M d cos i ≈ . × exp[ − ( R −
17 kpc) / . M (cid:12) kpc − for R >
17 kpc . (6)The dust mass surface density in Figure 3(b) shows two distinct peaks, corresponding to ringsat R = 5 . R = 11 . R = 15 . . I α ; Devereux et al. 1994),and a ∼
40% overdensity of stars with ages > ∼ I gasshows a similar deficiency in the same region.The DL07 dust model has A V = 0 . (cid:18) Σ M d M (cid:12) kpc − (cid:19) mag . (7)At R ≈
20 kpc, the dust mass surface density (projected onto the disk of M31) Σ M d cos i ≈ × M (cid:12) kpc − , corresponds to a visual extinction A V ≈ .
07 mag normal to the disk of M31, or A V ≈ . R = 11 . M d cos i ≈ × M (cid:12) kpc − , correspondingto A V ≈ . A V ≈ . M d ≈ . M (cid:12) kpc − , corresponding to A V ≈ Hubble SpaceTelescope observations of M31 (Dalcanton et al. 2012) may allow measurement of reddening towardmany individual stars in M31; it will be of great interest to compare the stellar reddening valueswith 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 Σ M d ,while for R >
17 kpc we use the annular SEDs to estimate the dust mass in each annulus. We finda total dust mass M d = 5 . × M (cid:12) within R = 25 kpc (see Table 4).It is difficult to estimate objectively the uncertainty in the estimate for M d . Calibrationuncertainties in the photometry itself are at least 10% for each of the MIPS, PACS, and SPIREbands. 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 assumptionsmade in the modeling about the physical properties of the dust, and simplified treatments of thestarlight intensity distribution. Overall, we adopt a tentative uncertainty estimate of 20% for theglobal dust mass given in Tables 4 and 5. The coefficient 0.74 differs from the value 0.67 in Aniano et al. (2012) because of the mass renormalizationdiscussed in Section 3.2. R . We find thathalf of the dust mass in M31 lies at R > . >
10% lies beyond 19 kpc.Table 5: Estimates for the Dust Mass in M31 M d (10 M (cid:12) ) a Data used Reference2 . ± . IRAS , extinction, H I Xu & Helou (1996)3 . ± . IRAS µ m Haas et al. (1998)4 IRAS , COBE-DIRBE,
ISO , MIPS Gordon et al. (2006)5 . ± .
45 PACS 100,160; SPIRE 250,350,500; annuli HELGA I Fritz et al. (2012)2 . . ± . ; M160 pixels and annuli This work a For D = 744 kpc.The HELGA survey (Fritz et al. 2012; Smith et al. 2012) obtained PACS and SPIRE imagingof M31 at a scan speed of 60 (cid:48)(cid:48) s − (three times faster than the scan speed for the PACS and SPIREimaging used in the present study). HELGA I (Fritz et al. 2012) measured the fluxes in a centralcircular aperture and five concentric annuli extending out to 20 kpc. Except for the central circle,the annular boundaries were ellipses. The λ ≥ µ m SED of each annulus was fit with a modifiedblackbody F ν ∝ ν β B ν ( T ), with fixed β = 1 .
8. The total dust mass M d = (5 . ± . × M (cid:12) found within 20 kpc by HELGA I is close to the present estimate M d ( R <
20 kpc) = (5 . ± . × M (cid:12) .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 MIPS70photometry 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 wereacceptable provided they did not exceed the MIPS70 photometry. In the regions included in theirstudy, HELGA II found a total dust mass M d = 2 . × M (cid:12) , but Smith et al. (2012) note thatthe pixels satisfying their 5 σ criterion accounted for only ∼
50% of the global 500 µ m emission fromM31. 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 usehere. Integrating this map out to R = 25 kpc yields an H I mass M (H I) = 6 . × M (cid:12) (88% ofthe total H I in their map). CO (1–0) has been mapped out to ∼
12 kpc, with the associated H Surface density Σ H = Σ(H I) + Σ(H ) from H I 21cm (Braun et al. 2009) and CO 1-0 (Nieten et al.2006, using X CO = 2 × cm − ( K km s − ) − ) at S350 resolution (left) and M160 resolution (right). Most ofthe structure seen at S350 resolution survives at M160 resolution. Note the deficiency of gas in the SW at radii15 kpc < R <
20 kpc. mass estimated to be M (H ) ≈ . × M (cid:12) (Nieten et al. 2006). We take the total H surfacedensity to be Σ H ≈ Σ(H I) + Σ(H ), where Σ(H ) is obtained from the observed CO 1-0 emissionassuming a constant X CO = 2 × H cm − ( K km s − ) − .We may be underestimating Σ H . H that is “CO dark” (i.e., not associated with CO 1-0emission) is assumed to be a small fraction of the total H mass. More importantly, we do notinclude H II gas in our estimate for Σ H . The ionized gas mass in bright H II regions is small, butthe mass in low-density diffuse H II may not be negligible. For our Galaxy, it is estimated thatdiffuse H II accounts for ∼
23% of the total ISM mass at
R <
20 kpc (Draine 2011b). The centerof M31 has Σ H ≈ . × M (cid:12) kpc − averaged over R < α ∼ × − erg cm − s − sr − (Devereux et al. 1994; Tabatabaei & Berkhuijsen2010) corresponds to an H II surface density 6 × ( cm − / (cid:104) n e (cid:105) ) T . M (cid:12) kpc − , where (cid:104) n e (cid:105) is theelectron density within the emitting regions. The density-sensitive [S II ]6716/[S II ]6731 line ratioindicates (cid:104) n e (cid:105) > . × cm − for R (cid:46) . ∼ Kplasma appears to be small compared to that in neutral gas. X-ray observations indicate that R (cid:46) T ≈ × K and a mean surface density ∼ × M (cid:12) kpc − (Bogd´an & Gilfanov 2008), only ∼
10% of the
R < I mass. Using a conversion factor X CO ( J = 1 →
0) = 2 × H cm − ( K km s − ) − (Bolatto et al. 2013). Dust to Gas mass ratio, areas with Ldust > 1.5e6 Mo/kpc2Dust to Gas mass ratio, areas with Ldust > 1.5e6 Mo/kpc2 M ( D u s t ) / M ( H ) Fig. 6.— (a) Dust-to-H mass ratio at M160 resolution, for pixels within the Σ L d = Σ L d , min contour. (b) Dust-to-Hmass ratio for each M160 resolution pixel with Σ L d > Σ L d , min . The central curve is the mean in each radial bin; theother curves are the mean ± σ , where σ is the estimated variance of the distribution. It is also possible that we may have overestimated Σ H in the center – Leroy et al. (2011) findthat X CO in the center ( R < R (cid:46) X CO .Figure 5 shows Σ H at S350 and M160 resolution. Note the deficiency of gas in the SW for R > . R = 16 and R ≈
20 kpc)is not known, but this part of the M31 disk has the appearance of having been affected by somerecent event. Block et al. (2006) argued that a nearly head-on collision with M32 ∼
210 Myr agocan account for the observed offset of the center of the 11 kpc ring, but such an encounter wouldnot seem likely to produce the observed deficiency of H I and dust in the SW at R ≈ −
20 kpc.It may be the result of an encounter with another member of the Local Group within the pastGyr. Lewis et al. (2013) discuss the Giant Stellar Stream, extending from the SW side of the disktoward 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 theGiant Stellar Stream and the H I filamentary structure is uncertain, but strongly suggestive of aclose passage of M33, possibly accounting for the deficiency of interstellar matter on the SW sideof the disk, ∼
18 kpc from the center of M31.Figure 6(a) is a map of the dust-to-H mass ratio at M160 resolution. Some azimuthal structure9Fig. 7.— (a) 100 × dust/H mass ratio in M31 as a function of galactocentric radius R . Triangles: present work,modeling at M160 resolution, using only pixels with Σ L d > Σ L d , min . Circles: present work, modeling the SED for∆ R = 677 pc annuli. Squares: HELGA I (Fritz et al. 2012). Dot-dash line: HELGA II (Smith et al. 2012). (b)Metallicity relative to solar, based on our dust modeling (triangles and circles, as in (a)), and H II region oxygenabundances (Zurita & Bresolin 2012); the line is shown as solid in the 8–16 kpc range where the oxygen abundancesare from the “direct” method. is evident, but the main feature is a conspicuous radial trend, with the dust-to-H ratio peaking atthe center and declining with R . Figure 6(b) shows results for every pixel (at M160 resolution)satisfying the condition Σ L d > Σ L d , min . Because of the low signal/noise ratio (S/N) in individualpixels, the scatter in the derived M d /M H is quite pronounced for R (cid:38)
17 kpc. However, the meandust/H ratio in each annulus exhibits a clear trend, decreasing with increasing R . The centraldust-to-H ratio is M d /M H ≈ . ∼ R ≈
20 kpc: a factor of ∼ R ≈
20 kpc.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 ringsfor
R <
23 kpc, including only pixels satisfying the criterion Σ L d > Σ L d , min , and calculating theratio of (total dust)/(total gas) for that ring. The dust/gas ratio therefore applies only to the pixelssatisfying the Σ L d > Σ L d , 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 theannulus do not satisfy the cut. This incompleteness is due in part to the general radial declinein surface brightness, but in part it reflects the conspicuous deficit of both gas and dust in theSW 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, and0the gas mass is the total mass of gas within the annulus. The dust-to-gas ratio is well-behavedout 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 thepixels with Σ L d > Σ L d , min . For R >
18 kpc the dust/gas ratios for the pixels with Σ L d > Σ L d , min is somewhat higher than the result from the annular photometry. This bias is attributable tothe fact that for Σ L d just below Σ L d , min , noise can raise the pixel above the threshold, leading tooverestimation of the dust mass. In this regime, dust mass estimates from the annular photometryand modeling should be more reliable.We find that the dust/H ratio for 0–25 kpc declines monotonically with increasing R , with M d M H ≈ . − R/ . R < . − R/
19 kpc) 8 kpc < R <
18 kpc0 . − R/ < R (cid:46)
25 kpc , (8)shown as a dashed line in Figure 7. Also shown in Figure 7(a) is the dust/gas ratio estimated byHELGA I (Fritz et al. 2012) for each of 5 radial zones. The HELGA I results are in fair agreementwith our findings interior to ∼ M dust /M H ≈ .
049 exp( − R/ . ∼
6. Metallicity of the ISM in M31
If depletions are similar to the local Milky Way, we expect M d M H ≈ . ZZ (cid:12) . (9)where Z is the mass fraction of elements other than H and He. This allows us to estimate themetallicity from our measured dust/H mass ratio (8): ZZ (cid:12) ≈ .
08 exp( − R/ . R < .
81 exp( − R/
19 kpc) 8 kpc < R <
18 kpc6 .
65 exp( − R/ < R (cid:46)
25 kpc . (10)This relation is plotted in Figure 7(b).Zurita & Bresolin (2012, hereafter ZB12) measured elemental abundances in M31 H II regionsusing “direct” methods between 8 and 16 kpc. When they allow for depletion of oxygen into grainsand a bias against H II regions with high oxygen abundances, they estimate that(O / H) ≈ . / H) (cid:12) exp( − R/
19 kpc) . (11)1Fig. 8.— (a) Starlight heating rate parameter U min in M31 at S350 resolution. (b) Radial profile of U min , estimatedat S350 resolution and M160 resolution (triangles) and using annular photometry (circles). Locations of dust masssurface density maxima are indicated. However, it is important to note that the H II region abundance determinations have appreciableuncertainties, and do not agree well with metallicity determinations in the atmospheres of B su-pergiants. Our result (Equation (10)) for the metallicity is in excellent agreement with the ZB12metallicity over the 8 kpc < R <
16 kpc range where the ZB12 H II region metallicities were basedon direct determinations of the gas temperature, and are therefore most reliable.According to Equation (10), the ISM in M31 has supersolar abundances for R (cid:46)
11 kpc. Thisis consistent with the high WC/WN stellar ratio observed in M31 (Neugent et al. 2012).
7. Dust Temperature and Starlight Properties
If the distribution of both stars and dust were known, the intensity and spectrum of thestarlight 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 thecomplex and correlated spatial distributions of both stars and dust.To model the infrared emission from dust, we take a much simpler approach, empirical approachto the starlight heating. Within a single “pixel” (which may include ∼ pc of disk area) thedust may be exposed to a wide range of starlight intensities, ranging from the general starlightbackground in the diffuse ISM, to high intensities found in star-forming regions. The DL07 model2Fig. 9.— Left: map of mean starlight heating rate parameter (cid:104) U (cid:105) = ¯ U in M31 at S350 resolution. Right: radialprofile of characteristic dust temperature. Locations of dust surface density maxima are indicated. adopts a parameterized distribution function for the starlight heating rate: the dust within a givenpixel 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 tothe data, we in effect use the dust grains as photometers to determine the intensity of starlight inthe regions where dust is present. The parameter U min is interpreted as being the starlight heatingrate in the diffuse ISM. The mean (weighted by dust mass) starlight heating rate within a pixel is (cid:104) U (cid:105) . 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 ≈ . 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).For starlight with the interstellar radiation field spectrum of Mathis et al. (1983, hereafterMMP83) the dust heating rate parameter U is U = u (cid:63) . × − erg cm − , (12)where u (cid:63) is the starlight energy density. For the DL07 model, the characteristic grain temperature(of the grains dominating the emission at λ > µ m) is related to the heating rate parameter U as T d , char ≈ U / K . (13)3Fig. 10.— (a) Radial profile of (cid:104) U (cid:105) in M31 from modeling the IR SED at S350 resolution and at M160 resolution(triangles), and using annular photometry (circles). The dashed curves show the estimated contributions to theheating from disk and bulge stars (see text). (b) Radial profile of (cid:104) U (cid:105) in the central 3 kpc. This is only a representative temperature – dust grains of different sizes and compositions illumi-nated by a single radiation field have different steady-state temperatures, and very small grainsundergo temperature fluctuations due to the quantized heating by stellar photons. Figure 9(b)shows T d , char as a function of radius.The mean starlight heating rate (cid:104) U (cid:105) is shown in Figure 9a. At S350 resolution, the centerhas (cid:104) U (cid:105) ≈
32 (see Figure 9), corresponding to T d , char ≈
33 K. At R = 15 kpc we find U ≈ .
3, or T d , char ≈
15 K, and at R = 20 kpc we find U ≈ .
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 Groveset al. (2012), who modeled the 100 − µ 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 higherthan the value 12 . R > (cid:104) U (cid:105) in Figure 9(a) is presumablyprimarily due to starlight from disk stars. In the star-forming parts of the disk, the radiationfield heating the dust is from a mixture of young and old stars, modified by dust attenuation; (cid:104) U (cid:105) 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 thethree-dimensional distributions of stars and dust. The observed (cid:104) U (cid:105) in Figure 10(a) shows an4approximately exponential decline with increasing R out to 20 kpc, which can be approximated by (cid:104) U (cid:105) disk ≈ . − R/
12 kpc) . (14)This is plotted in Figures 10; we will use it to estimate the contribution of disk starlight to (cid:104) U (cid:105) inthe central regions.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, isdiscussed by Groves et al. (2012), and elaborated further in Appendix B. From a three-dimensionalmodel for the stellar bulge, the energy density of bulge starlight is estimated to be u (0) (cid:63), bulge ( R ) = 5 . × − I ( R/r b ) erg cm − , (15)where I ( x ) is given by Equation (B4), r b = 0 .
58 kpc is the core radius of the bulge (see AppendixB) and the superscript (0) indicates that it is a theoretical estimate, with dust extinction ne-glected. For the starlight spectrum of the bulge, we estimate U (0)bulge = u (0) (cid:63), bulge . × − erg cm − = 30 I ( R/r b ) . (16)Equation (16) gives U (0)bulge ( R = 1 kpc) = 7 .
9, whereas our estimated heating rate from the IR SEDat R = 1 kpc is (cid:104) U (cid:105) ≈ . ∼ .
92. Thus at R = 1 kpc we might estimatethe heating contribution of the bulge stars to be ∼ .
5, 80% of our theoretical estimate U (0)bulge atthis radius. At smaller radii we find that U (0)bulge exceeds the inferred (cid:104) U (cid:105) by ∼ U (0)bulge and (cid:104) U (cid:105) estimated from modeling the IR SED: (1) Perhaps the bulge starlight hassimply been overestimated by ∼ (cid:104) U (cid:105) from the infrared observations, we have assumed the dust grains to have the properties of dust inthe DL07 model. The dust near the center of M31 may well differ from the solar-neighborhooddust on which the DL07 is based. If the dust grains at R (cid:46) lower ratio of (opticalabsorption cross section)/(FIR absorption cross section) than the DL07 model grains, a given radi-ation 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 Our estimate for u (cid:63), bulge (Equation (15)) is slightly below the value u (cid:63), bulge ( R ) = 5 . × − I ( R/r b ) erg cm − obtained by Groves et al. (2012). ∼ (cid:104) U (cid:105) , we will adopt the spatial profile expected forthe bulge starlight, but will scale down the heating rate by a factor 0 .
7. To this we add our empiricalestimate (Equation (14)) for the heating due to the disk stars. Thus: (cid:104) U (cid:105) = U bulge + (cid:104) U (cid:105) disk (17) U bulge = 0 . U (0)bulge = 21 I ( R/r b ) , r b = 0 .
58 kpc (18) (cid:104) U (cid:105) disk ≈ . − R/
12 kpc) . (19)This estimate of U is plotted in Figures 9 and 10. The agreement is good, except at R < . R < . R = 0 . ∼
15% lower than the T d , char values found here. Their lower dust temperature estimates may bedue to use of an earlier version of the data reduction pipeline and calibration factors, and possiblyalso to use of a modified blackbody rather than the multicomponent physical grain model usedhere.To summarize, we find relatively good agreement between the observed dust emission spectrumand that which would be expected for heating by a combination of the bulge starlight and adisk heating component. The observed dust temperatures are only slightly below what would beexpected. The inferred heating rate from the bulge stars is only ∼
30% lower than predicted for asimple 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 forsuch a discrepancy, this agreement is gratifying.The dust temperature T depends on the starlight heating rate parameter U and on the ratio (cid:104) C abs (cid:105) (cid:63) ( a ) / (cid:104) C abs ( a ) (cid:105) T , where (cid:104) C abs ( a ) (cid:105) (cid:63) is the dust absorption cross section averaged over thespectrum of the illuminating starlight for a grain of radius a , and (cid:104) C abs ( a ) (cid:105) T is the Planck-averagedabsorption cross section for grain temperature T . The fact that the observed dust temperature iswithin 5% of the predicted dust temperature indicates that the actual values of (cid:104) C abs (cid:105) (cid:63) / (cid:104) C abs (cid:105) T are close to the values in the DL07 grain model. This builds confidence in our grain model, andhence in the dust mass estimates, which are proportional to ρa/ (cid:104) C abs ( a ) (cid:105) T . Given the extremeenvironmental differences, it is remarkable that the dust properties near the center of M31 appearto be so similar to values inferred for dust in the solar neighborhood.6
8. PAH Abundance
The intensity in the IRAC bands includes both direct starlight and emission from dust pop-ulation. The observed intensities at λ ≥ . µ m are modeled as the sum of a stellar component(modeled as a 5000K blackbody) plus a nonstellar component ( F ν ) ns contributed primarily byPAHs. 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.8and 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 masscontributed by PAHs containing fewer than 10 C atoms. In the model fitting, q PAH is essentiallyproportional to the ratio of the power in the 6.2 and 7.7 µ m PAH emission features divided bythe total power radiated by the dust. When only IRAC photometry is available for λ < µ m,the q PAH parameter is, in practice, proportional to the nonstellar contribution to IRAC8.0 dividedby the aggregate 70–350 µ m luminosity. For a fixed starlight spectrum, F (8 µ m) is approximatelyproportional to q PAH and to the total infrared power, F TIR :( νF ν ) IRAC8 , ns ≈ A (cid:63) q PAH F TIR , (20)where ( νF ν ) IRAC8 , ns is the non-stellar contribution to the IRAC 8.0 µ m band, and the dimensionlesscoefficient A (cid:63) depends on the spectrum of the illuminating starlight. Draine & Li (2007) show that A MMP83 ≈ .
72 for the MMP83 spectrum of the starlight in the solar neighborhood.Varying the spectrum of the starlight can cause A (cid:63) to change, for two reasons. (1) The PAHabsorption cross section depends on wavelength differently from the absorption of overall grainmixture ( ∝ F TIR ), hence the fraction of the starlight power that is absorbed by PAHs will dependon the spectrum. (2) The PAH emission spectrum is the result of single-photon heating, hencethe fraction of energy absorbed by PAHs that is reradiated in the IRAC8.0 band depends on theilluminating spectrum. As an example of the dependence of A (cid:63) on the illuminating spectrum,Draine (2011a) showed that A ≈ . . A MMP83 for a 20kK blackbody cut off at 13 . A bulge ≈ . U = U bulge + U MMP83 , the effectivevalue is the dust luminosity-weighted mean A (cid:63) = A bulge U bulge + A MMP83 U MMP83 U bulge + U MMP83 = 1 . U bulge + 4 . U MMP83 U bulge + U MMP83 . (21)To correct the estimate of q PAH made assuming the MMP83 radiation field, we will take( q PAH ) corr = 4 . A (cid:63) × ( q PAH ) MMP83 (22)where ( q PAH ) MMP83 is the value of q PAH estimated assuming the dust to be heated by starlight withthe MMP83 spectrum, with U bulge and U MMP83 given by Equations (18) and (19).7Fig. 11.— (a) Map of PAH abundance parameter ( q PAH ) MMP83 in M31 at S350 resolution. The low values of q PAH in the NE may be due to problems with background subtraction in IRAC5.8 and IRAC8.0 (see text). (b) Radialprofiles: ( q PAH ) MMP83 (open symbols) and ( q PAH ) corr from Equation (22) (filled symbols). A map of ( q PAH ) MMP83 at S350 resolution is shown in Figure 11(a). Figure 11(b) shows theradial profile of ( q PAH ) MMP83 and ( q PAH ) corr . Interestingly, q PAH appears to peak in the 11 kpcring, attaining a value ( q PAH ) corr ≈ .
049 that is close to the value q PAH ≈ .
045 estimated for thediffuse ISM in the solar neighborhood.In Figure 11, ( q PAH ) MMP83 declines as one approaches the center, reaching a value ( q PAH ) MMP83 ≈ .
02 in the central regions. However, this decline is largely an artifact of assuming that the spectrumof the illuminating starlight is independent of radius, which is incorrect – in the central regionsthe starlight is dominated by light from the bulge stars, which is much redder than the MMP83spectrum.Figure 11(b) shows that when we allow for the starlight being increasingly dominated by anold stellar population as we move to the center, ( q PAH ) corr shows only limited variation as we movefrom R = 11 kpc to the central kpc. Evidently the balance between PAH formation and destructionin the ISM remains relatively constant from the central kpc out to R = 20 kpc. There does appearto be a systematic radial decline in q PAH for
R >
11 kpc, but this again might be an artifact of aradial gradient in the spectrum of the starlight as one moves from the 11.2 kpc ring – where starformation is active – to outer regions where there appears to be little contemporary star formation.At R (cid:38)
20 kpc q PAH estimated from the ∆ R = 677 pc annuli appears to rise. However, wesuspect this to be an artifact of imperfect background subtraction in the 5 . . µ m images.8Fig. 12.— For M160 resolution modeling: ratio of the model intensity divided by the observed intensity for theSPIRE 250 and 500 µ m bands. The DL07 model successfully fits the 500 µ m emission out to the edge of the “galaxymask”, at R ≈
20 kpc. In the inner regions, the DL07 model tends to underpredict SPIRE 250 µ m and overpredictSPIRE 500 µ m (see also Figure 13), but the deviations are only at the ∼
10% level.
Although background subtraction works well in other bands out to ∼
25 kpc, it appears less suc-cessful for the IRAC 5.8 µ m and 8 . µ m bands. The difficulty with background subtraction maybe 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 (Gordonet al. 2011) was able to minimize these problems by combining images taken with very different rollangles together with custom processing techniques, but we are simply using the M31 images fromBarmby et al. (2006). It is also possible that the Galactic foreground has structure arising fromvariations in PAH abundance or ionization state on ∼ . ◦ scales, which would not be identified bythe background subtraction procedures used here.
9. Dust Properties
In the present work, we attempt to reproduce the observed SED in each pixel using the DL07dust model and a parameterized distribution of starlight intensities. Above we have examinedthe values of the dust modeling parameters, such as the dust mass, q PAH , and properties of thestarlight intensity distribution. Here we compare the models with observations to see how wellthe model actually reproduces the data, and whether any systematic deviations are present thatindicate systematic problems with the modeling.9Fig. 13.—
Opacity spectral index β obs from Equation (23) obtained from modeling at M160 resolution (triangles)and using annular photometry (circles). The DL07 model overpredicts the SPIRE500/SPIRE250 band ratio by afactor > . R ≈ R ≈ β obs > β DL07 + 0 .
26. For R (cid:38)
10 kpc the DL07model agrees with the observations to within 10%, and the model appears to have nearly the correct value of β . Alsoshown is the HELGA II (Smith et al. 2012) result for β (see text). Aside from minor effects associated with variation in the PAH abundance when q PAH is allowedto vary, the composition (amorphous silicate and graphitic carbon) is assumed to be the sameeverywhere in the DL07 models used to fit the infrared emission. Further, the dust opacity isassumed to be independent of temperature. The DL07 dust model used here has an opacity with afixed dependence on frequency at long wavelengths ( λ (cid:38) µ m). Here we compare the model to theobserved emission from M31 to look for residuals that might be indicative of differences betweenthe 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 showsthe ratio of model to observation at SPIRE250 and SPIRE500. The first impression is that themodel is good: the ratio of model to observation is generally between 0.86 and 1.16, which seemsgood in view of noise in the observations, uncertainties in calibration, and the general uncertaintiesin the adopted dust opacities. However, systematic trends are evident: in the central regions ofM31 (excluding the center itself) the model tends to be low at SPIRE250, and high at SPIRE500.Let β ≡ ln [ κ (250 µ m) /κ (500 µ m)] / ln 2 be the effective power-law index of the dust opacitybetween 250 and 500 µ m. For the DL07 model (see Table 2) this ratio is β DL07 = 2 .
08. If the fitteddust temperatures were left unchanged, the 500/250 flux ratio could be brought into agreement0with observations if β were changed to β obs = β DL07 + ln ([ I ν (250 µ m) /I ν (500 µ m)] obs / [ I ν (250 µ m) /I ν (500 µ m)] model )ln(500 / . (23)Figure 13 shows β obs as a function of R for R <
22 kpc, using only M160 resolution pixels withΣ L d > Σ L d , min . For 1 < R <
10 kpc, β obs is larger than β DL07 , indicating that the opacity between250 and 500 µ m should fall more rapidly than ν . .The dust opacity ratio κ (250 µ m) /κ (500 µ m) could vary with location because the dust com-position is varying, or it could conceivably result from variations of the grain temperature, if thedust opacities are temperature-dependent. The DL07 model assumes the dust opacities to be in-dependent of temperature. However, some materials do exhibit temperature-dependent opacitiesin the laboratory (e.g., Mennella et al. 1998; Boudet et al. 2005; Coupeaud et al. 2011), and suchbehavior is expected in some models of amorphous solids (Meny et al. 2007). A number of stud-ies 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 fromboth observational noise and line-of-sight temperature variations (Shetty et al. 2009a,b; Kelly et al.2012). The “two level system” (TLS) model (Meny et al. 2007; Paradis et al. 2011) predicts thatincreasing dust temperature T should lead to a lower value of β . However, in M31 it appears thatthe 1–10 kpc regions where the dust is warmer than the outer disk have a larger value of β thanthe value in the outer disk, and the center – where the dust is hottest – has essentially the samevalue of β as the relatively cold dust at R ≈
20 kpc. The observed variations in β do not seem tobe consistent with what would be expected for the TLS model, unless one allows for substantialradial variations in the TLS model parameters themselves. Rather than attributing the apparentchanges in β to temperature, it would appear instead that the dust composition must be varyingwith radius in M31.The HELGA II collaboration (Smith et al. 2012) reported a sharp change in the dust propertiesat R ≈ β obtained from their modifiedblackbody fits. Interior to 3 kpc, the inferred dust temperature decreased from T = 27 . T = 16 . R = 3 . β 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 . 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 (cid:46) R (cid:38) µ m excess at large R – quite the contrary, the DL07 model tends to slightly over predict SPIRE500 in the outer regionsof M31.The HELGA II analysis was based on observations of M31 that were shallower than the dataused here, and also used an earlier calibration of the three SPIRE bands. It is possible that thedifferences between the HELGA II results and those of the present study may be due in part to1differences in signal/noise, calibration or data reduction pipeline. In the outer regions results arealso sensitive to background subtraction.If we were over- or under-subtracting the IR backgrounds, we would obtain very low or veryhigh 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 orCO 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 goodmatch 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.
10. Summary
The principal conclusions of this work are as follows:1. Consistent with previous observations, we find that the dust mass surface density in M31peaks in two rings, at R = 5 . R = 11 . R ≈ . M d = (5 . ± . × M (cid:12) within R = 25 kpc. 95% of this dustmass lies within R = 21 kpc, with the dust surface density peaking at R = 11 . R , from ∼ .
027 atthe center to ∼ . R = 25 kpc.4. The dust/H mass ratio parallels measurements of O/H in H II regions, consistent with aconstant fraction of the refractory elements Mg, Si, and Fe being in dust. Based on ourestimated dust/H mass ratio, we infer that the metallicity Z/Z (cid:12) varies from ∼ R = 0 to ∼ . R = 25 kpc – see Figure 7(b) and Equation (10).5. The starlight heating rate parameter (cid:104) U (cid:105) shows a nearly monotonic decline with galactocentricradius, from (cid:104) U (cid:105) ≈
50 at the center (at S350 resolution) to (cid:104) U (cid:105) ≈ . R ≈
20 kpc.6. We confirm the finding of Groves et al. (2012) that the dust heating in the central 2 kpc isdominated by light from the stellar bulge. We find that the starlight heating rates inferredfrom the observed IR emission are consistent with the heating rates expected from the bulgestarlight. It is remarkable that the dust properties near the center of M31 appear to be similarto the properties of dust in the solar neighborhood.7. After taking into account variation in the spectrum of the starlight, we find the PAH abun-dance q PAH to be approximately constant from the center of M31 out to ∼
20 kpc. There is2some 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 aresome systematic deviations. The dust at R ≈ − . < β < .
33 in the 250–500 µ m wavelength range, whereas the dust at R > β ≈ .
08, consistent with the dust in the DL07 model. We are in approximate agreementwith the radial variation of β found by HELGA II (Smith et al. 2012) in the central ∼ β values for R > R (cid:38)
10 kpc the DL07 model, with fixed opacity, is consistent with the observedphotometry, and returns dust masses that are consistent with the observations of the gas andmetallicity.We thank Edvige Corbelli, Stephen Eales, Jeremy Goodman, Matthew Smith, and the anony-mous 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.
A. PACS versus MIPS
MIPS and PACS have two wavelengths in common: 70 µ m and 160 µ m. While the filter responsefunctions are not identical, they are similar enough that we would expect MIPS and PACS tomeasure very similar flux densities for smooth SEDs. However, comparisons of MIPS and PACSimaging often shows discrepancies that are much larger than expected – see, e.g., the cases ofNGC 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.The PACS70/MIPS70 comparison is particularly striking: while there are a few regions nearthe center with PACS70/MIPS70 ∼ >
2. The fact thatthis occurs over extended regions makes it clear that this is not simply a consequence of randomnoise. The fact that it is non-uniform over the image makes it clear that it is not simply aresult of calibration. Background subtraction is of critical importance, but we have employedthe same background estimation procedures (Aniano et al. 2012) for both MIPS and PACS. It isnow thought that MIPS70 suffers from sub-linear behavior in high-surface brightness regions, butsuch high surface brightnesses are found only at the center of M31, whereas we see high values ofPACS70/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 com-parison, but there are still many regions – particularly along the 11 kpc ring – where PACS160/MIPS1603Fig. 14.—
Ratios of images convolved to the MIPS160 PSF, displayed within the “galaxy mask” where thesignal/noise is high enough that dust modeling is possible. Left: ratio of PACS70/MIPS70. Right: ratio ofPACS160/MIPS160. > .
4. This difference is much larger than the claimed uncertainties in the observed intensities atthese locations. MIPS160 is now thought to have sublinear response for I ν (cid:38)
50 MJy sr − (Paladiniet al. 2012), but most of M31 is below this value. This again emphasizes the value of being ableto include the MIPS160 data in the analysis, possible only if one degrades all the other imaging tothe MIPS160 PSF (FWHM 39 (cid:48)(cid:48) ). B. Heating by Bulge Starlight
Groves et al. (2012) have discussed the contribution of M31’s bulge stars to the heating ofdust. For the adopted bulge stellar luminosity density profile (Geehan et al. 2006) ρ (cid:63) ( R ) = L bulge πr b R/r b )(1 + R/r b ) , (B1)the starlight energy density due to bulge stars (neglecting possible absorption or scattering by dust)4is u (cid:63), bulge ( R ) = L bulge πr b c I ( R/r b ) (B2) I ( x ) ≡ x (cid:90) ∞ ln (cid:16) x + y | x − y | (cid:17) (1 + y ) dy (B3)= 1 x − x )2 x (1 + x ) − x ( x − . (B4) I ( x ) is logarithmically divergent for x →
0. Near x = 1, I ( x = 1 + (cid:15) ) = 12 + ln 28 − . x −
1) + O (( x − ) . (B5)The projected luminosity interior to radius R is L ( < R ) = (cid:90) R ρ (cid:63) ( r )4 πr dr + (cid:90) ∞ R ρ (cid:63) ( r ) (cid:104) − (cid:112) − ( R/r ) (cid:105) πr dr (B6)= L bulge (cid:34) − α (cid:90) π/ cos θdθ (1 + α sin θ ) (cid:35) α ≡ r b R (B7)= 0 . L bulge for α = 0 .
61 (B8)We obtain L bulge by integrating the average of the intrinsic and reddened spectrum obtainedby Groves et al. (2012) for L ( < . µ m , µ m]. Groveset al. (2012) assumed D = 780 kpc and r b = 0 .
61 kpc. Corrected to D = 744 kpc, we obtain L bulge (0 . − µ m) = 1 . × L (cid:12) , with r b = 0 .
58 kpc. Thus, the 0 . µ m < λ < µ m energydensity due to the bulge stars is u (cid:63), bulge ( R ) = 5 . × − I ( R/r b ) erg cm − . (B9)This energy density of bulge starlight gives a theoretical estimate for the dust heating rate, nor-malized to the heating rate for the MMP83 radiation field, U (0)bulge ( R ) = 30 I ( R/r b ) . (B10)For R = 1 kpc this gives U (0)bulge (1 kpc) = 7 .
9. This is our zero-th order estimate for the dust heatingrate in the central few kpc of M31. In Section 7 we show that this theoretical estimate for the dustheating rate agrees well with the dust heating rate inferred from the IR SED.
REFERENCES
Aniano, G., Draine, B. T., Calzetti, D., et al. 2012, ApJ, 756, 46Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. M. 2011, PASP, 123, 12185Barmby, P., Ashby, M. L. N., Bianchi, L., et al. 2006, ApJ, 650, L45Block, D. L., Bournaud, F., Combes, F., et al. 2006, Nature, 443, 832Bogd´an, ´A., & Gilfanov, M. 2008, MNRAS, 388, 56Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207Boudet, N., Mutschke, H., Nayral, C., et al. 2005, ApJ, 633, 272Braun, R., Thilker, D. A., Walterbos, R. A. M., & Corbelli, E. 2009, ApJ, 695, 937Chemin, L., Carignan, C., & Foster, T. 2009, ApJ, 705, 1395Ciardullo, R., Rubin, V. C., Jacoby, G. H., Ford, H. C., & Ford, W. K., Jr. 1988, AJ, 95, 438Corbelli, E., Lorenzoni, S., Walterbos, R., Braun, R., & Thilker, D. 2010, A&A, 511, A89Coupeaud, A., Demyk, K., Meny, C., et al. 2011, A&A, 535, A124Crane, P. C., Dickel, J. R., & Cowan, J. J. 1992, ApJ, 390, L9Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18de Vaucouleurs, G., de Vaucouleurs, A., Corwin, H. G., Jr., et al. 1991, Third Reference Catalogueof Bright Galaxies (New York: Springer)Devereux, N. A., Price, R., Wells, L. A., & Duric, N. 1994, AJ, 108, 1667Draine, B. T. 2003, ARA&A, 41, 241Draine, B. T. 2011a, in EAS Publications Series, Vol. 46, PAHs and the Universe, ed. C. Joblin &A. G. G. M. Tielens, 29Draine, B. T. 2011b, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: PrincetonUniv. Press)Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866Draine, B. T., & Li, A. 2001, ApJ, 551, 807Draine, B. T., & Li, A. 2007, ApJ, 657, 810Dupac, X., Bernard, J.-P., Boudet, N., et al. 2003, A&A, 404, L11Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10Fritz, J., Gentile, G., Smith, M. W. L., et al. 2012, A&A, 546, A34Geehan, J. J., Fardal, M. A., Babul, A., & Guhathakurta, P. 2006, MNRAS, 366, 9966Gordon, K. D., Bailin, J., Engelbracht, C. W., et al. 2006, ApJ, 638, L87Gordon, K. D., Meixner, M., Meade, M. R., et al. 2011, AJ, 142, 102Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3Griffin, M. J., North, C. E., Schulz, B., et al. 2013, MNRAS, 1306.1778Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892Haas, M., Lemke, D., Stickel, M., et al. 1998, A&A, 338, L33Habing, H. J., Miley, G., Young, E., et al. 1984, ApJ, 278, L59Hora, J. L., Fazio, G. G., Allen, L. E., et al. 2004, in Society of Photo-Optical InstrumentationEngineers (SPIE) Conference Series, Vol. 5487, 77Kelly, B. C., Shetty, R., Stutz, A. M., et al. 2012, ApJ, 752, 55Kennicutt, R. C., Calzetti, D., Aniano, G., et al. 2011, PASP, 123, 1347Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12Lewis, G. F., Braun, R., McConnachie, A. W., et al. 2013, ApJ, 763, 4Li, A., & Draine, B. T. 2001, ApJ, 554, 778Liang, Z., Fixsen, D. J., & Gold, B. 2012, submitted to MNRAS(arXiv:1201.0060)Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212Mennella, V., Brucato, J. R., Colangeli, L., et al. 1998, ApJ, 496, 1058Meny, C., Gromov, V., Boudet, N., et al. 2007, A&A, 468, 171Neugent, K. F., Massey, P., & Georgy, C. 2012, ApJ, 759, 11Nieten, C., Neininger, N., Gu´elin, M., et al. 2006, A&A, 453, 459Paladini, R., Linz, H., Altieri, B., & Ali, B. 2012, PACS ICC Document, PICC-NHSC-TR-034Paradis, D., Bernard, J. P., M´eny, C., & Gromov, V. 2011, A&A, 534, A118Paradis, D., Veneziani, M., Noriega-Crespo, A., et al. 2010, A&A, 520, L8Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013a, ArXiv e-prints, 1307.6815Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013b, arXiv:1303.5062, 1303.50627Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013c, ArXiv e-prints, 1303.5072Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, A&A, 536, A19Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2Popescu, C. C., & Tuffs, R. J. 2013, ArXiv e-prints, 1305.0232Reach, W. T., Megeath, S. T., Cohen, M., et al. 2005, PASP, 117, 978Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25Roussel, H. 2013, PASP, 125, 1126Shetty, R., Kauffmann, J., Schnee, S., & Goodman, A. A. 2009a, ApJ, 696, 676Shetty, R., Kauffmann, J., Schnee, S., Goodman, A. A., & Ercolano, B. 2009b, ApJ, 696, 2234Smith, M. W. L., Eales, S. A., Gomez, H. L., et al. 2012, ApJ, 756, 40Tabatabaei, F. S., & Berkhuijsen, E. M. 2010, A&A, 517, A77Vilardell, F., Ribas, I., Jordi, C., Fitzpatrick, E. L., & Guinan, E. F. 2010, A&A, 509, A70Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1Xu, C., & Helou, G. 1996, ApJ, 456, 163Zurita, A., & Bresolin, F. 2012, MNRAS, 427, 1463