EDP Sciences
Home All issues Volume 561 (January 2014) A&A, 561 (2014) A95 Full HTML
Free Access
Issue
A&A
Volume 561, January 2014
Article Number A95
Number of page(s) 14
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201321441
Published online 10 January 2014

© ESO, 2014

1. Introduction

Dust and gas are thoroughly mixed in the interstellar medium (ISM), with approximately half of the metals (atomic number Z > 2) in each. Unlike molecules, which radiate at specific frequencies, dust emission covers all frequencies and carries much more energy. Dust emission thus provides an alternative to spectral lines for studying the mass distribution in the ISM. Dust has the advantage over the commonly used CO molecule because it is not photo-dissociated by UV photons. However, dust emission depends strongly on the grain temperature, such that it is not straightforward to measure the dust mass. To estimate the grain temperature or range in temperatures, it is necessary to know how grain emission varies with wavelength.

For dust grains in local thermal equilibrium, the temperature is usually obtained using a modified black body (MBB) emission given by (1)where Sν is the intensity, Bν represents the Planck function, and τν denotes the optical depth of the dust, depending on frequency ν. In the optically thin limit, the above equation converts to (2)where Σdust is the dust mass surface density and κν = κ0   (ν/ν0)β is the dust emissivity with index β (κ0 is the grain cross-section per gram at frequency ν0).

There is some evidence supporting variations of the dust emission spectrum with environmental conditions from both Galactic (e.g. Paradis et al. 2009) and extra-galactic studies (e.g. Lisenfeld et al. 2000). The dust emissivity index β could change depending on grain properties such as structure, size distribution, or chemical composition and is still a matter of debate. These properties may be affected by different physical/environmental processes like shattering, sputtering, grain  −  grain collisions (mainly due to shocks, see e.g. Jones 2004, and references therein), condensation of molecules onto grains, and coagulation (e.g. Draine 2006). The extra-galactic observations presented here sample a beam of 160 pc with a depth of a few 100 pc (Combes et al. 2012). Along each line-of-sight, we therefore expect to sample regions of widely different physical properties such as temperatures and densities. The observed dust temperatures and β indices are therefore weighted averages of the local conditions along the lines of sight. This is discussed in more detail in Sect. 6. Neglecting the variation in the dust emissivity index could be misleading when determining the dust mass (e.g. Malinen et al. 2011). For instance, Reach et al. (1995) found a very massive cold-dust component in the Milky Way (MW), assuming a constant β, using the long-wavelength COBE data but this was not confirmed by Lagache et al. (1998).

Nearby galaxies provide ideal astrophysical laboratories (1) with different chemical abundances and (2) where  −  for face-on galaxies  −  the galactic disk is only crossed once by our line of sight, which makes it straightforward to distinguish spiral arms, associate gas (+dust) clouds with star-forming regions, and inter-arm regions.

The Herschel Space Observatory (Pilbratt et al. 2010) provides far-infrared (FIR) and sub-millimeter (submm) data at high angular resolution and sensitivity over the wavelengths required to investigate the physical properties of dust in galaxies. Recent studies have suggested that the emission of dust grains could differ from the standard MBB model (see Eq. (1)), which assumes a β = 2 emissivity law. This can lead to unphysical gas-to-dust mass ratios compared with those expected from the metallicity of the galaxies (Meixner et al. 2010; Galametz et al. 2010; Galliano et al. 2011; Israel & Maloney 2011). The emissivity indices reported from observations range from β = 1 to β = 2.5 (Chapin et al. 2011; Casey et al. 2011; Boselli et al. 2012). Theoretically, β < 1 is excluded at FIR and submm wavelengths, otherwise the Kramers-Kronig relation does not hold. The Mie theory implies that β = 2 for spherical grains of idealized dielectrics and metals, while β = 1 if the complex optical constant does not depend on frequency (Krügel 2003).

There is a systematic degeneracy between β and the dust temperature T given by Eq. (2), with an amplitude that depends on the signal-to-noise ratio (S/N), and the spectral coverage. Therefore, β and T cannot be reliably determined for a single position by fitting the FIR-submm dust emission, even for a homogeneous and single-temperature cloud (e.g. Kelly et al. 2012). This is because equally good fits (i.e. well within observational noise even for high S/N data, see e.g. Fig. 3 in Planck Collaboration 2011c) can be obtained for different (β, T) pairs. The degeneracy is such that β increases as T decreases, and vice versa. The scatter in the fits can produce an apparent anti-correlation between β and T (Shetty et al. 2009a,b).

The power of this degeneracy is illustrated in Fig. 1. For a modified black body of temperature 18 K and index β = 1.8, we added noise with a standard deviation of 1, 2, and 5% (left pannels from top to bottom) of the signal, corresponding to S/Ns of 100, 50, and 20. The mock data points were fit and the fits were plotted; this was repeated 100 times. The red triangles are the original no-noise data points at 160, 250, 350, and 500 microns. The right column shows the derived β and T values and the large open star indicates the (T = 18 K, β = 1.8) point. Despite the very high S/Ns, the scatter is broad and follows a banana-shaped anti-correlation pattern. The anti-correlation is generated by the noise despite the unrealistically high S/N levels used in these plots. For the S/Ns of 20, 50, and 100 used in Fig. 1, the uncertainties (standard deviation) are 3 K, 1.1 K, and 0.5 K in temperature and 0.4, 0.15, and 0.08 in β, respectively. For T = 18 K and β = 1.8, the data can be fit by a variety of values varying such that β ≈ 1.8−0.155(T − 18). In Sect. 4, we present a more detailed statistical (Monte Carlo) analysis showing how we constrain the dust temperature and β.

thumbnail Fig. 1

Illustration of the β − T degeneracy: the panels on the left-hand side show fits to the 4 data points to which noise has been added, and the right column shows the values of β and T. The large triangles on the left show the original data points (before adding noise), which follow a modified black-body law with β = 1.8 and T = 18, typical values for normal spiral galaxies. The large open stars show this position in the right-hand column. From bottom to top, the noise added follows a Gaussian distribution with standard deviation of 5% (bottom panels), 2% (middle panels) and 1% (top panels). Despite this very low noise level, the scatter in β and T is significant and correlated. An enlargement of the spectral coverage would help reduce the degeneracy if no other processes such as emission from small or spinning dust grains play a role within the extended wavelength range.
Open with DEXTER

Several authors have found an anti-correlation between β and T (e.g. Dupac et al. 2003). While many authors (e.g. Paradis et al. 2009; Ysard et al. 2012) believe the anti-correlation is not entirely due to the β − T degeneracy, it remains difficult to isolate and identify the physical and degenerate parts of this relation (e.g. Galametz et al. 2012; Shetty et al. 2009a).

Laboratory experiments have also found an anti-correlation between β and T. However, in space, the anti-correlation is typical for dust temperatures of 10−20 K (the dominant temperatures of dust emitting in the submm), whereas laboratory experiments find the difference to be between 300 K dust (room temperature) and 10−20 K dust. The difference in dust spectra between 10 K and 20 K is very difficult to identify in laboratory experiments (e.g. Coupeaud et al. 2011). It is currently unclear whether laboratory measurements support a physical β − T anti-correlation for dust in the 10 − 20 K range.

Table 1

Positional data on M 33.

The nearest Scd galaxy, M 33 (NGC 598), at a distance of 840 kpc (1″ ≃  4 pc, Freedman et al. 1991), has been extensively studied at IR wavelengths. Its mild inclination (i = 56°, Regan & Vogel 1994, Table 1) makes it a suitable target for mapping a broad variety of astrophysical properties. Herschel observations of M 33 (HerM33es1, Kramer et al. 2010) enable us to study the dust spectral energy distribution (SED) at a resolution where complexes of giant molecular clouds and star-forming regions are resolved. Using the Herschel PACS and SPIRE data together with the Multiband Imaging Photometer Spitzer (MIPS) fluxes, Kramer et al. (2010) studied the dust SED in radial intervals of 2 kpc width (out to 8 kpc galactocentric radius) as well as the integrated SED of M 33. Assuming a constant β, they found that a two-component MBB fits the SEDs better than the isothermal models. They also showed that the global SED fits better with β = 1.5 than with β = 2. In a pixel-by-pixel approach, Xilouris et al. (2012) derived the temperature and luminosity density distribution of a two-component MBB for the best-fit β of 1.5, which was fixed across the galaxy. This paper investigates possible variations of β radially and in the disk of M 33, using different single- and two-component MBB approaches in which the degeneracy between β and T is treated properly. Studying both single- and two-component MBBs is important to assess the temperature mixing along the line of sight, which could in principle lead to a lower β value if a single-MBB model is assumed (Malinen et al. 2011). To obtain as much information as possible about the warm dust and to avoid contamination by nonthermalized grains (e.g. very small grains) as much as possible, we used the Herschel SPIRE/PACS and Spitzer MIPS maps at wavelengths λ ≥ 70   μm, probing primarily the dust emission from the big grains that are believed to be in thermal equilibrium with the interstellar radiation field (see e.g. Carey et al. 2009). In the diffuse Galactic ISM, these grains emit as an MBB with an equilibrium temperature of 17 − 18 K (e.g. Boulanger et al. 1996).

The paper is organized as follows. The data sets used in this study are described in Sect. 2. We present our methods and approaches as well as their resulting physical parameters in Sect. 3. Then, we investigate the systematics of the two-component MBB approaches by performing a Monte Carlo simulation (Sect. 4). Uncertainties of the physical parameters are presented in Sect. 5. We discuss and summarize the results in Sects. 6 and 7.

Table 2

Images of M 33.

thumbnail Fig. 2

Histogram of the reduced χ2 in the single-MBB fitting approach for selected β values between 1.4 and 1.8. Results are shown for fits to the high (20σ, left) and very high (60σ, right) S/N data.
Open with DEXTER

2. Data

Table 2 summarizes our data. M 33 was observed at 250 μm, 350 μm, and 500 μm using the SPIRE photometer (Griffin et al. 2010) onboard the Herschel space telescope. The data were reduced using the Herschel data reduction software HIPE 7.0 as detailed in Xilouris et al. (2012). We also used the 160 μm data taken with the Herschel PACS instrument (Poglitsch et al. 2010). The PACS data were reduced using the scanamorphos algorithm (Roussel 2013) as discussed in detail in Boquien et al. (2011). At shorter wavelengths, we used the 70 μm Spitzer MIPS data in the two-component MBB approach2.

The largest calibration uncertainties (i.e., for diffuse emission) are 15% for the MIPS (Carey et al. 2009), 20% for the PACS (Poglitsch et al. 2010), and 15% for the SPIRE bands (Griffin et al. 2010).

We also used Hα, HI, and CO(2−1) line emission data to investigate the connection between the dust spectrum and tracers of the ionized and neutral gas phases. The CO(2−1) line emission was observed with the IRAM-30 m telescope and is detailed in Gratier et al. (2010). The HI-21 cm line was mapped with the VLA (Gratier et al. 2010). The Hα data are taken from the Kitt Peak National Observatory (KPNO; Hoopes & Walterbos 2000). The maps used in our analysis were all convolved to the resolution of the 500 μm SPIRE image (~40″) by using the dedicated convolution kernels provided by Aniano et al. (2011), which were projected on the same grid of 10″ pixel and center position.

3. Methodology and results

We independently constrained β and T by using (a) single-component MBB fits for which we assumed that the grains are described by a fixed β but a varying temperature, and all physically reasonable values of β were tested to determine which β results in the lowest residuals and (b) two-component MBB models to which we applied an iterative method to solve systems of independent equations using the Newton-Raphson method. Approach (a) aims at deriving β and T as a function of the galactocentric radius and searches for trends with FIR luminosity. We derived β, T, and mass surface densities pixel by pixel across the galaxy in approach (b).

3.1. Single-component approach

In the single-component MBB model (Eq. (1)), it is assumed that the grain properties do not change over the region or structure studied, therefore we fit only the temperature for each pixel (and the optical depth, which acts as a general scaling). The temperature was fit for a fixed β but we repeated the fitting process for many different values of β. For a given β, the overall quality of the fit was determined by minimizing the reduced χ2 as in Xilouris et al. (2012) over the whole set of pixels3. While making the hypothesis that β, which represents the grain properties, is unlikely to undergo real variations at the scale of the regions studied, this process avoids the β − T degeneracy while allowing us to determine the best β.

For each pixel, fits were first made for the 160 − 500 μm bands where cold dust dominates (Compiègne et al. 2011). At 100 μm, the emission comes from a mixture of cold and warm dust. The observed 100 μm emission thus represents an upper limit to the emission from the cool component, which we fit from 160 − 500 μm. When the 100 μm flux derived from the fit was higher than the observed flux, we refitted using the 100 μm flux as well, which resulted in a more realistic dust temperature.

thumbnail Fig. 3

Radial variation of the best-β calculated by minimizing the reduced χ2 for each of the annuli in the single modified black-body (MBB) fitting approach (top). The different symbols show a coherent decrease in β with radius and represent values obtained for different flux cutoffs as indicated. Bottom: dust temperature calculated for each flux cutoff at each radius. The error bars show the 1 sigma dispersion.
Open with DEXTER

Exploring the best β for a selected set of β = 1, 1.3, 1.5, 1.7, and 2, Xilouris et al. (2012) showed that the data from all wavelengths above the 3σ rms noise level are best fitted by β = 1.5. Similarly, we explored how the best β value varies for progressively higher S/N data. While choosing only high-quality data is desirable as long as enough pixels are available (which is the case), it introduces a selection of regions with more and more star formation and probably also a higher fraction of molecular gas than in the lower S/N regions. Compared with the S/N = 3 case (Xilouris et al. 2012, see Fig. 3 in), where β = 1.5 yields the lowest residuals, we found that the more actively star-forming high S/N regions have steeper β, for instance, β = 1.6 for S/N = 20 and β = 1.7 for S/N = 60 (see Fig. 2).

To better investigate the change in β, we divided M 33 into radial bins of 0.5 kpc and made the same calculations for each β between 1 and 2.4 with steps of 0.1 and determined the lowest residuals. Figure 3 shows that the best value for β decreases with radius from 1.8 to about 1.2. Though differently, the dust temperature also decreases with radius. The higher flux-cutoff (10 MJy/sr) has a higher temperature for two reasons. First, dust tends to be warmer at higher fluxes. Second, since the β calculated at large radii are lower for the 10 MJy/sr cutoff, the temperatures obtained are higher. These results are very similar to those of Fig. 1 in Tabatabaei et al. (2012), which uses entirely different numerical techniques.

3.2. Two-component approaches

Here we introduce a second dust component heated mainly by young massive stars to a higher temperature (warm dust) in addition to the dust heated by the general interstellar radiation field (ISRF, cold dust). The aim of this approach is to differentiate the dust properties in arm/interarm and inner/outer disk regions, including diffuse emission. Hence, we performed this analysis pixel by pixel using the most sensitive Herschel and Spitzer maps available in the wavelength range of 70 μm to 500 μm. The specific intensity of a two-component MBB dust emission is given by (3)where the indices c and w stand for the cold and warm dust, respectively. The Newton-Raphson iteration method from the Numerical Recipes was used to solve the five equations corresponding to the five wavelengths at 70, 160, 250, 350, and 500 μm to derive Tc, Tw, Σdust,c, Σdust,w, and β.

To decrease the probability of falling into a local minimum of the merit function, the initial guess for the physical parameters was determined by synthesizing the observed intensities. We first generated uniform samples of 5000 random sets of variables corresponding to each physical parameter: cold and warm dust-temperatures, mass surface densities, and β. Each random set of physical parameters was then used to synthesize the intensities given by Eq. (3) at each wavelength. To define the initial values for the Newton-Raphson method, 20 parameter sets leading to the synthesized intensities closest to the observed values were selected by maximizing the likelihood. Thus, the Newton-Raphson method was applied 20 times for each pixel of the maps, resulting in 20 sets of solutions out of which the set that minimized the merit function ( ≡ , where σobs is the corresponding calibration uncertainty) was taken as the final set of solutions. The initialization scheme is completely independent for every two pixels. Hence, there is no correlated error due to this procedure.

Considering the cold and warm dust, different treatments for β were used:

  • I)

    Both components of dust are assumed to have the same variableemissivity index (βc = βw = β).

  • II)

    Because most of the dust was in the cool phase, the variable β is attributed to this major component (βc = β), while the warm dust was assumed to emit with a fixed emissivity index across the disk. Although we had a fixed value of βw (note that the number of unknowns must be equal to the number of equations), we repeated the analysis for different possibilities of βw = 1, 1.5, and 2.

These calculations were all performed for every pixel of the maps with intensities higher than the 3σ noise level at each wavelength. The pixel-by-pixel analysis furthermore helps to reduce dust-temperature mixing because the range of temperatures is limited in a given pixel. Using a Monte Carlo simulation, the two-component MBB approaches were then synthesized to determine the systematics and reliability of each method. The temperature and surface density of the cold dust do not depend on the choice of βw, but the physical parameters of the warm dust are best reproduced using the βw = 2 model (see Sect. 4). This agrees with Casey (2012), who derived a mid-infrared power-law slope of  ≃ 2, indicative of optically thin dust with a shallow radial density profile from clumpy, hot regions. Hence, we focused on the results of the two-component MBB with βw = 2, although the general distributions of the physical parameters across the galaxy are similar to those of other two-component MBB models. The probability density function (PDF) of residuals between the observed and modeled intensities (Fig. 4) also supports the preference for the βw = 2 case for the warm dust component because it results in a smaller uncertainty U4 at 70 μm than the other two-component MBB models (see the U values in Fig. 4, upper panel).

thumbnail Fig. 4

Probability density function (PDF) of residuals between observations and model for the two-component MBB approaches at 70, 160, 250, 350, and 500 μm (from top to bottom). The uncertainty of each approach (U) is also given for comparison.
Open with DEXTER

thumbnail Fig. 5

Distribution of the cold dust-mass surface density in M 33. A sketch of the optical spiral arms (Sandage & Humphreys 1980) is overlaid on the cold dust map. The wedge gives the surface density units in μg cm-2.

Open with DEXTER

3.2.1. Distribution of cold and warm dust

Figures 5 and 65 show the resulting distributions of the mass surface densities of the cold and warm dust for the βw = 2 case. The cold dust has been found not only in the main spiral arms (arm I in the north and south), but also in the weaker arms, flocculent structures, and dust lanes. Dense patches of cold dust with surface densities of Σdust,c > 45 μg cm-2 occur along arm I and some parts of arms II and III (both in the north and south). Moreover, filaments of dense dust have been detected in the central 2.5 kpc. The densest region corresponds to the high concentration of dust in the giant HII/star-forming complex NGC 604.

The warm dust is mostly concentrated in the central 2.5 kpc and along the main spiral arm I and shows clumpy structures. A good match between clumps of warm dust and star-forming/HII regions is found by a comparison with the Hα emission (Fig. 6, right). Apart from the central region and the main arms, the rest of the galaxy is covered by many spotty warm dust features, most of which can be found in the Hα map as well. Therefore not only the bright HII regions (such as NGC 604, IC 133, NGC 595, NGC 588, and B690), but also the weaker ones can be traced in the warm dust map.

thumbnail Fig. 6

Left: distribution of the warm dust-mass surface density in M 33. Overlaid are contours of the Hα emission, shown in black or white for better contrast against the background. The wedge gives the surface density units in μg cm-2. Right: the Hα map of M 33.
Open with DEXTER

In addition to the tight spatial correlation between Hα emission and the warm dust in Fig. 6, Fig. 7 shows that the warm dust column and the Hα emission are proportional, with IHα ∝ (0.92 ±    0.02)Σdust,w. The correlation is best (highest Pearson correlation coefficient, see Sect. 6) when βw = 2, and this is one of the reasons for preferring the βw = 2 model over βw = 1 and 1.5. The Monte Carlo simulation presented in Sect. 4 provides additional support for βw = 2.

The average mass surface density of the warm dust in star-forming complexes in the central 4 kpc is Σdust,w = 0.8 ± 0.3   μg cm-2. The cold-to-warm dust-mass ratio has a Gaussian distribution (in log) that peaks at  ~ 100, typical of spiral galaxies (Misiriotis et al. 2006; Tabatabaei et al. 2013), and a standard deviation of one order of magnitude. We stress that this result was derived without any pre-assumption about the filling factor of the cold or warm dust.

3.2.2. Distribution of the dust emissivity index β

Figure 8 shows that β varies between 1.2 and 2 across the galaxy, in agreement with Fig. 3. The mean value of the dust emissivity index over the disk is β = 1.5 ±    0.2. This agrees with β derived globally by Kramer et al. (2010) using the integrated SED and Xilouris et al. (2012) by exploring the best-fit β described in Sect. 3.1. As in the inner disk, a steep dust spectrum with β ≃ 2 is found in NGC 604. Out to NGC 604, the main spiral arms are visible in the β map, because they have higher index (β > 1.5) than their adjacent inter-arm regions with lower β. This confirms the higher β for the high S/N cutoffs derived for these radii in the single-component MBB approach (Sect. 3.1). There are small-scale fluctuations with variations of  ≲ 0.2 all over the β map which could be due to noise or to the systematics of the approach (Sects. 4 and 5).

3.2.3. Cold and warm dust temperatures

The temperatures of both cold and warm dust components are in general lower in the outer disk than in the inner disk of M 33. In the central 4 kpc, the cold dust temperature Tc varies mainly between 20 K and 24 K with an average of 21 ± 2 K (Figs. 9 and 3). In the outer parts, Tc ranges between 14 K and 19 K with a mean value of 17 ± 2 K. The higher cold dust temperature in the inner disk with respect to the outer disk is consistent with reports by Kramer et al. (2010), Braine et al. (2010), and Xilouris et al. (2012). A similar trend was also found in other disk galaxies based on Herschel data (e.g. Pohlen et al. 2010; Engelbracht et al. 2010; Fritz et al. 2012). The corresponding Tc values in HII regions agree with those derived by Relano et al. (2013) using the Draine & Li (2007) models and following the prescription of Galametz et al. (2012). Warm dust temperatures of Tw ≳ 50 K are found in the spiral arms and star-forming regions with a maximum of 60 ± 10 K in NGC 604, NGC 595, and IC 133. In the outer parts, the warm dust has a temperature of Tw = 35 ± 10 K on the average (the errors given are discussed in Sect. 5).

Based on the physical parameters retrieved, the modeled SEDs are shown for two selected pixels with different β in Fig. 10.

4. Monte Carlo simulations with the Newton-Raphson method

Using a Monte Carlo simulation, we assessed the systematics and reliability of the two-component approaches. This is particularly useful in determining a possible artificial correlation between β and T. The Monte Carlo simulations were performed assuming a uniform probability density function. We first generated random numbers for each of the physical parameters in a selected range given in Table 3. For surface densities, we used a uniform distribution in logarithmic space. The extreme values in the selected ranges may seem to be very different from reality. Compared with more moderate values, however, this selection ensures that we covered the whole physical parameter space. Each set of physical parameters generated (input) leads to synthesized intensities using Eq. (3), as described in Sect. 3.2. The synthesized intensities were then treated as observed data in the Newton-Raphson method, resulting in new sets of physical parameters (output). In the ideal case, that is, when there are no systematics caused by a specific model or approach, the output parameters are similar to the input parameters. This case is indicated by a one-to-one (equality) relation between the input and output parameters in Figs. 11 and 12.

thumbnail Fig. 7

Linear correlation between the warm dust-mass surface density Σdust,w and the Hα emission measure. The scatter in the low surface-density tail coincides with the systematic of the approach for this range of surface density (see Sect. 4).
Open with DEXTER

In all the two-component MBB models, the simulations show that Tc can be best reproduced between 10 K and 20 K and the scatter around the equality relation is symmetric for the βw = 2 case (other models tend to overestimate Tc). In this case, Tc is reproduced with an accuracy better than 3 K in the range 10 − 20 K. This case also shows a symmetric scatter for Tw for 35 K  < Tw <  50 K, where it is best reproduced (with a dispersion of 11 K), while other models tend to underestimate Tw in this range. The mass surface density of the cold dust generally shows a small scatter for all models and is much better constrained than that of the warm dust. The accuracy of reproducing the dust mass surface density is about 5.5 times better for the cold than for the warm component in logarithmic scale for the βw = 2 case. This case is again more successful in reproducing Σdust,w than other cases, especially for mass surface densities Σdust,w < 10-6 g cm-2. Interestingly, we found a coincidence between the low-density scatter in the Hα − Σdust,w correlation plot (Fig. 7) and the low-density deviations in the simulated Σdust,w (Fig. 12).Where Tc < 10 K and Σdust,c > 10-5 g cm-2, some points follow a shallower slope than the equality relations. These points are most likely local minima in the merit algorithm. However, synthesizing the observed data reduces the chance of falling into a local minimum, as explained in Sect. 3.2. Figure 11 shows that β is not as well constrained as Tc. To give an example, for an input value of β = 1.5 (with βw = 2), the output value of β varies between 1.1 and 2 with a dispersion of 0.3. The method tends to overestimate β for β < 1.5 and underestimates it for β > 2.

To study the sensitivity of the different MBB approaches to noise, we added Gaussian noise with a width of 5σ rms of the observed maps to the synthesized intensities at various wavelengths and repeated the simulations6. Adding the noise increased the dispersion in β by  ~ 50% in the two-component MBB with βw = 2. It is important to note that none of the model approaches led to an artificial correlation between temperature, mass surface density, and β for either simulated or real results (see Fig. 13 and Sect. 5.1).

thumbnail Fig. 8

Distribution of the dust emissivity index β in the disk of M 33 obtained using the two-component approach with βw = 2.
Open with DEXTER

thumbnail Fig. 9

Distribution of the cold (left) and warm (right) dust temperatures in M 33. The wedge shows the temperature values in Kelvin.
Open with DEXTER

5. Sources of uncertainties

In addition to the uncertainties due to the method, our results might be affected by statistical and systematic uncertainties on measured fluxes. To study the effect of the statistical uncertainties, we first measured the stochastic pixel noise by calculating the dispersion in intensities in background regions. In the reduced images, background dispersion might in principle be due to the propagation of original pixel noise and artifacts into the background pixels as well as possible imperfect foreground subtraction. We therefore simulated the observed intensities by adding random noise with a normal distribution of a 3σ background noise level measured at each wavelength. Using this statistical experiment, we created 1000 mock datasets used as the data points (fluxes) in our two-component approach with βw = 2, which resulted in a set of 1000 values for each physical parameter. The statistical uncertainties (δstat) were then determined from the standard deviation of the results (Table 4).

thumbnail Fig. 10

Modeled dust SED for the two-component modified black body for 2 different pixels with coordinates RA = 1h 33m 31.35s, Dec = 30:30:21.18 (top) and RA = 1h 33m 17.31s, Dec = 30:21:47.99 (bottom). The dashed, dot-dashed, and solid curves represent the SEDs of the warm, cold, and total dust emission, respectively.

Open with DEXTER

The flux uncertainty might also be due to systematic uncertainties such as calibration errors and artifacts caused by observing cameras (e.g. Aniano et al. 2012), which do not have the same properties as the statistical errors. In most cases, systematics are addressed via calibration uncertainties, because the systematics due to cameras are not well-defined. For NGC 6946 and NGC 0628, Aniano et al. (2012) investigated the systematics due to the MIPS and PACS cameras by comparing their data at a same wavelength (after convolving them to the same resolution and pixel size). For M 33, this comparison is possible only at 160 μm, because at other wavelengths the data are taken with only one camera. We considered the relative difference in the 160 μm intensities obtained with MIPS and PACS as an estimate for the systematic error if it is larger than the calibration uncertainty at each pixel. In the region of interest R < 6 kpc (inside the optical radius R25 ~ 7.5 kpc), the relative difference between the MIPS and PACS data is  ≲ 20%, which is similar to the calibration uncertainty. Hence, the systematics at 160 μm are dominated by the calibration uncertainty. We assumed a similar situation for the data at other wavelengths. This is supported by using the largest calibration uncertainties at each wavelength7 (Sect. 2). The uncertainties in the physical parameters due to the systematics were then determined by adding the calibration uncertainties (see Sect. 2) to the synthesized intensities retrieved in the Monte Carlo experiment in Sect. 4. These intensities were then used to derive the physical parameters using the Newton-Raphson method. This way the systematics due to the method are also taken into account. Table 4 shows the median of the absolute residuals8 ( | output -input | ), δsys.

Table 3

Setup of the numerical experiment in the two-component MBB approaches.

thumbnail Fig. 11

Monte Carlo simulation of the two-component MBB model. From left to right: Tc, Tw, and βc( = β). From top to bottom: βw = 2, βw = 1.5, βw = 1, and βw = βc. The red line shows the equality between the output and input parameters.
Open with DEXTER

thumbnail Fig. 12

Monte Carlo simulation of the two-component MBB model. Left: Σdust,c and right: Σdust,w. From top to bottom: βw = 2, βw = 1.5, βw = 1, and βw = βc. The red line shows the equality between the output and input parameters.
Open with DEXTER

The error on each physical parameter is then given by . As expected, δ is dominated by the systematics and not by the statistical noise.

To evaluate the role of the systematics on the decrease of β with distance from the center shown in Sect. 3.2.2, yet another experiment was performed. After adding the systematic uncertainties to the observed intensities, a vertical cut passing through the center of M 33 was solved 100 times using the Newton-Raphson method. This resulted in 100 sets of solutions for each pixel along the vertical cut9. The uncertainties in the physical parameters were then determined from the standard deviations of the parameter distributions. Figure 14 (left) shows the distribution (profile) of β and its uncertainty δβ along the cut. From the center to R ≃    6 kpc (indicated by the blue lines), δβ increases from  ≃ 0.15 to  ≃ 0.3. This is equivalent to a relative change in β from  ≃ 6% to  ≃ 20% as shown in Fig. 14 (right). Hence, the magnitude of the derived radial decrease of beta is significant, that is, larger than three times the uncertainty.

thumbnail Fig. 13

Top: scatter plot of the cold dust emissivity index β vs. temperature, colour-coded for three radial intervals. Bottom: scatter plot of the cold dust temperature vs. the warm dust temperature showing no correlation.
Open with DEXTER

6. Discussion

6.1. The β − T correlation

One of the most important results of this work is that we clearly show that not only the temperature but also β decreases with galactocentric distance in M 33. This is a real effect, shown by different physical and numerical approaches, and additionally, the effect is strong enough to be seen against the well-known numerical β − T anti-correlation (Shetty et al. 2009b).

The techniques adopted in this work – solving a system of equations or deriving the best-β leaving only the temperature as a variable – are robust despite the strength of the β − T degeneracy. However, for any given pixel there is a tendency for β and T to be anti-correlated simply because dust temperatures derived with lower β are higher, whatever the method. When averaged over radial annuli, it becomes clear that both β and T decrease with radius, though their decreasing trends are different (e.g. Fig. 3).

In M 31, with similar data and thus similar linear resolution as for M 33, Smith et al. (2012) found that T and β show an opposite behavior as a function of radius (their Fig. 8), such that the temperature increases with radius and β decreases. This may be linked to the β − T degeneracy because the dust temperature decreases with radius in M 31, as shown by Fritz et al. (2012) using the same data, which was also reported in other studies (Groves et al. 2012; Tabatabaei & Berkhuijsen 2010).

It is known that the dust temperature decreases with radius in M 33 due to the radial decrease in the interstellar radiation field at all wavelengths (Verley et al. 2009; Kramer et al. 2010; Braine et al. 2010; Boquien et al. 2011; Xilouris et al. 2012). A trend of decreasing β could only be expected if β is also linked to a factor that varies (on average) with radius. In addition to the ISRF, many properties of spiral galaxies vary with distance from the center: the metallicity decreases with radius, as does the H2/HI ratio, the stellar surface density, and the star formation rate. Thus, the gradient in β we found may be due to any of these factors, among which the metallicity, the H2/HI ratio, and the star formation rate are strongly linked. This motivated us to search for a correlation between β and different environmental conditions imposed by star formation activity and different phases of the ISM.

6.2. Connection to star formation and ISM tracers

Within the inner disk, the spiral arms and NGC 604 are visible as high-β regions (Fig. 8). In Figs. 3 and 8, a clear change in β is apparent at 4 kpc radius, where β changes from approximately Galactic (β ≃ 1.8, Planck Collaboration 2011a) to lower than that of the LMC (β ≃ 1.5, Planck Collaboration 2011b). This is also the radius at which Gardan et al. (2007) noticed a break in the H2 formation efficiency in M 33 (see their Fig. 13). Starting at the same radius, the [CII]/FIR ratio is also found to increase (Kramer et al. 2013).

In Fig. 15, we investigate possible connections between the observed β and star formation via the Hα brightness, the CO(2 − 1) and the HI line emission. After excluding pixels with a S/N lower than 3, the maps of the Hα, CO(2 − 1), and HI integrated intensities10 were used to cross-correlate with β in the same way as in Sect. 5.1. Figure 15 shows that β is higher in regions with star formation as traced by Hα emission (top panel). The link between β and the molecular gas (middle panel), which is the fuel for star formation, is slightly weaker and agrees with the observation that most but not all molecular clouds in M 33 host star formation (Gratier et al. 2010). The atomic gas surface density is only weakly linked to the current star formation and there is no clear correlation between HI surface density and β (lower panel). The fact that, unlike Hα emission and molecular gas, the atomic gas is not concentrated in the central disk and along the spiral arms could explain the weak β–HI correlation obtained.

Because both Hα and molecular gas trace star formation, this indicates that the dust spectrum may differ in star-forming and non-star-forming regions. Star formation could in principle influence the grain composition and size distribution (via shattering of dust grains by strong shocks or dust coagulation), which could modify the dust emissivity index β. Although such a difference has not been proven, it is possible that the properties of the dust in the atomic and molecular gas clouds differ because of the temperature and density differences.

Star formation replenishes the ISM with metals, for instance, through stellar winds of young massive stars and supernovae. Concerning the dust grains, this is likely to increase the abundance ratio of silicates to carbonaceous grains11 (Woosley et al. 2002). According to the standard dust models, β = 2 is consistent with silicate grains, while β = 1 for carbonaceous grains (e.g. Desert et al. 1990). Observations of the Milky Way show that silicates dominate the ISM (e.g. Molster & Waters 2003; Molster & Kemper 2005, crystalline silicates are not found often in the diffuse ISM, however) which is consistent with the observed dust emissivity index of close to 2 (β ≃ 1.8, Planck Collaboration 2011a). Galaxies with lower metallicity than the Milky Way, however, are found to have flatter SEDs (e.g. Galliano et al. 2005; Planck Collaboration 2011b), indicating the general role of metal abundance on the dust emissivity variation. Similarly, M 33, which has a metallicity lower than solar by about a factor of 2 (Magrini et al. 2010), shows a flatter average value of β (~1.5), higher in the central 4 kpc and lower in the outer parts. The fact that most of the young massive stars are concentrated in the central 4 kpc of M 33 while the evolved carbon stars are smoothly distributed across the disk (e.g. Verley et al. 2009) could provide a distinct difference in the dust composition and hence the emissivity index.

Table 4

Uncertainties of the obtained physical parameters.

The dust grains in the ISM are subject to a variety of destruction and modification processes such as vaporization, shattering, and coagulation, with the latter two processes being more important in the ISM (e.g. Tielens et al. 1994; Jones et al. 1996). Shattering and fragmentation of dust grains occurs in grain-grain collisions caused by shocks, which can shift the peak in the grain size distribution to lower values (Jones et al. 1996) and may change crystalline grains into amorphous grains (Vollmer 2009), leading to a flatter SED (Seki & Yamamoto 1980). The small grains are more dominant in the diffuse and low-density ISM, considering that they have been accelerated by shock waves traversing the ISM, than in the spiral arms with sites of grain growth in dense molecular clouds.

thumbnail Fig. 14

Distribution of the emissivity index β (right) and its relative uncertainty δβ/β (left) along a vertical cut through the center. The red and blue lines indicate the center of M 33 and R ≃ 6 kpc, respectively.
Open with DEXTER

thumbnail Fig. 15

Scatter plots of the dust emissivity index β versus the Hα, CO(1 − 0), and HI emission from M 33. The y-axis presents the logarithm of the intensities in units of cm-6 pc for the Hα emission and K km s-1 for the CO(1 − 0) and HI emssions. Also shown are Pearson correlation coefficients (r), the ordinary least-squares (OLS) fits, and the bisector fits.
Open with DEXTER

We note the possibility of systematic deviations of the calculated β values from the true β. As shown in Sect. 4, the dynamical range in the true (intrinsic) β could be broader than that obtained in the two-component MBB approach. The β values lower than  ~ 1.5 could be overestimated and those higher than  ~ 2 underestimated. The true variation in β may be stronger than that shown in Fig. 8.

In our observations, as probably in all astrophysical situations, the dust within the telescope beam is not at a single temperature. The emission spectrum of dust at similar but different temperatures is necessarily broader than single-temperature dust and, if fit by a single MBB, the resulting β will be lower than for the real grains. This depends strongly on the density, however (Malinen et al. 2011).

Could the temperature mixing be particularly dominant for the outer disk of M 33, leading to the radial decrease in β that we find? We have several arguments that suggest this is not the case.

  • Firstly, it would appear somewhat counter-intuitive to have thebroader mix of dust temperatures in regions with fewer dustgrains (the dust surface density is considerably lower in the outerdisk), particularly with respect to the central regions. This agreeswith Malinen et al. (2011), whoshowed that the observed β is closer to the true β in regions with lowerdust densities. A weaker temperature variation in the outer thanthe inner parts is visible in Fig. 10 bycomparing Tc and Tw in each of those regions.

  • Secondly, several flux cutoffs were used for the single-component model (Fig. 3) and, irrespective of the cutoff, a clear gradient in β is observed. By selecting pixels with strong fluxes, we select outer disk regions with a level of star formation similar to what is observed in the inner disk, and indeed, the dust is fairly warm. However, these pixels have a low β. Similarly, using the two-component model, the two hot giant HII regions in the outer parts, IC 132 and IC 133, appear as weak sources in the β map (Fig. 8), unlike the HII regions in the inner disk.

  • Thirdly, not only the single-component, but also the two-component MBB model, which is designed to resolve the temperature mixing at first order of approximation, leads to a radial decrease in β. Comparing β from the two methods, we deduced that the single-component method indeed leads to a lower β, particularly in the center. However, the difference in β is not larger in the outer than in the inner disk (see Fig. 16, which shows the difference in β obtained from the two methods vs. radius).

thumbnail Fig. 16

Difference in the dust emissivity index β derived using the two-component MBB and the single-component MBB models vs. radius.
Open with DEXTER

7. Summary

Using the Herschel SPIRE/PACS and Spitzer MIPS maps at λ ≥ 70 μm, we investigated the physical properties of dust using single- and two-component modified black-body models and studied the variation of the dust emissivity index β across M 33 at a linear resolution of 160 pc (pixels of  ~50 pc). Although the single- and two-component models used modified black-bodies, the procedures were very different and both have been designed to avoid the β − T degeneracy. The single-component model assumed that β does not vary rapidly from pixel to pixel at the scale observed in M 33 and, by testing many values of β, determined the value leading to the lowest residuals. The two-component models solved exactly for the parameters (T, β, surface density) for a combination of warm and cold dust for data covering 70 μm to 500 μm, whereas the single-component models used the data from 160 μm to 500 μm.

Both analyses found that (i) β is higher in actively star-forming regions; (ii) β is higher in the inner disk and decreases sharply beyond a galactocentric distance of 4 kpc; and of course (iii) that the dust temperature decreases with radius. Both β and T decrease with radius in M 33, contrary to reported anti-correlations.

In addition, the two-component model showed that β is higher along the spiral arms within the inner disk. The warm dust follows the Hα emission, indicating that the separation into warm and cold dust is appropriate. Proper behavior of the two-component model was also verified by extensive Monte Carlo simulations.

In an attempt to identify the reason for the variation of β, we found that β is correlated both with Hα intensity and with CO(2 − 1) line strength, but not significantly with the HI column density. Because CO and Hα emission decrease radially, we were unable to distinguish their effects on β. The radial decrease in β appears not to be due to an increasingly broad mixture of dust temperatures.

1

Herschel M 33 extended survey open time key project, http://www.iram.es/IRAMES/hermesWiki

2

The PACS data at 100 μm were not used because their sensitivity is lower than the MIPS data (see e.g. Aniano et al. 2012; Hinz et al. 2012). In the single-component MBB approach, we used the PACS data at 100 μm to reject SED fits resulting in intensities higher than the observed intensities at 100 μm.

3

We minimized over all bands; for four bands fitting two parameters (temperature and optical depth), this is equivalent to a reduced χ2.

4

U is defined as the mean of relative residuals over all pixels.

5

The maps of the physical parameters shown in Figs. 5, 6, 8, and 9 are slightly smoothed to improve legibility.

6

Note that Figs. 11 and 12 show the simulations before the noise was added.

7

The calibration uncertainties used are larger than those of point sources by more than a factor of two. They are also larger than those assumed by Aniano et al. (2012).

8

The residuals have a Gaussian distribution.

9

This cut includes  ~280 pixels for which Iλ > 3σ rms noise.

10

For the sake of consistency, regions not covered by the CO(2 − 1) observations (Gratier et al. 2010) were masked out in the Hα and HI maps before the analysis.

11

The low-mass, evolved carbon (AGB) stars are the main source of the carbonaceous grains in the ISM, particularly in low-metallicity environments (Lattanzio 1989).

Acknowledgments

We thank Eva Schinnerer for support and helpful discussions. We also thank the anonymous referee, Brent Groves, Hendrik Linz, Svitlana Zhukovska, and Amy Stutz for discussions and useful comments. F.S.T. acknowledges the support by the DFG via the grant TA 801/1-1.

References

All Tables

Table 1

Positional data on M 33.

In the text
Table 2

Images of M 33.

In the text
Table 3

Setup of the numerical experiment in the two-component MBB approaches.

In the text
Table 4

Uncertainties of the obtained physical parameters.

In the text

All Figures

thumbnail Fig. 1

Illustration of the β − T degeneracy: the panels on the left-hand side show fits to the 4 data points to which noise has been added, and the right column shows the values of β and T. The large triangles on the left show the original data points (before adding noise), which follow a modified black-body law with β = 1.8 and T = 18, typical values for normal spiral galaxies. The large open stars show this position in the right-hand column. From bottom to top, the noise added follows a Gaussian distribution with standard deviation of 5% (bottom panels), 2% (middle panels) and 1% (top panels). Despite this very low noise level, the scatter in β and T is significant and correlated. An enlargement of the spectral coverage would help reduce the degeneracy if no other processes such as emission from small or spinning dust grains play a role within the extended wavelength range.
Open with DEXTER
In the text
thumbnail Fig. 2

Histogram of the reduced χ2 in the single-MBB fitting approach for selected β values between 1.4 and 1.8. Results are shown for fits to the high (20σ, left) and very high (60σ, right) S/N data.
Open with DEXTER
In the text
thumbnail Fig. 3

Radial variation of the best-β calculated by minimizing the reduced χ2 for each of the annuli in the single modified black-body (MBB) fitting approach (top). The different symbols show a coherent decrease in β with radius and represent values obtained for different flux cutoffs as indicated. Bottom: dust temperature calculated for each flux cutoff at each radius. The error bars show the 1 sigma dispersion.
Open with DEXTER
In the text
thumbnail Fig. 4

Probability density function (PDF) of residuals between observations and model for the two-component MBB approaches at 70, 160, 250, 350, and 500 μm (from top to bottom). The uncertainty of each approach (U) is also given for comparison.
Open with DEXTER
In the text
thumbnail Fig. 5

Distribution of the cold dust-mass surface density in M 33. A sketch of the optical spiral arms (Sandage & Humphreys 1980) is overlaid on the cold dust map. The wedge gives the surface density units in μg cm-2.

Open with DEXTER
In the text
thumbnail Fig. 6

Left: distribution of the warm dust-mass surface density in M 33. Overlaid are contours of the Hα emission, shown in black or white for better contrast against the background. The wedge gives the surface density units in μg cm-2. Right: the Hα map of M 33.
Open with DEXTER
In the text
thumbnail Fig. 7

Linear correlation between the warm dust-mass surface density Σdust,w and the Hα emission measure. The scatter in the low surface-density tail coincides with the systematic of the approach for this range of surface density (see Sect. 4).
Open with DEXTER
In the text
thumbnail Fig. 8

Distribution of the dust emissivity index β in the disk of M 33 obtained using the two-component approach with βw = 2.
Open with DEXTER
In the text
thumbnail Fig. 9

Distribution of the cold (left) and warm (right) dust temperatures in M 33. The wedge shows the temperature values in Kelvin.
Open with DEXTER
In the text
thumbnail Fig. 10

Modeled dust SED for the two-component modified black body for 2 different pixels with coordinates RA = 1h 33m 31.35s, Dec = 30:30:21.18 (top) and RA = 1h 33m 17.31s, Dec = 30:21:47.99 (bottom). The dashed, dot-dashed, and solid curves represent the SEDs of the warm, cold, and total dust emission, respectively.

Open with DEXTER
In the text
thumbnail Fig. 11

Monte Carlo simulation of the two-component MBB model. From left to right: Tc, Tw, and βc( = β). From top to bottom: βw = 2, βw = 1.5, βw = 1, and βw = βc. The red line shows the equality between the output and input parameters.
Open with DEXTER
In the text
thumbnail Fig. 12

Monte Carlo simulation of the two-component MBB model. Left: Σdust,c and right: Σdust,w. From top to bottom: βw = 2, βw = 1.5, βw = 1, and βw = βc. The red line shows the equality between the output and input parameters.
Open with DEXTER
In the text
thumbnail Fig. 13

Top: scatter plot of the cold dust emissivity index β vs. temperature, colour-coded for three radial intervals. Bottom: scatter plot of the cold dust temperature vs. the warm dust temperature showing no correlation.
Open with DEXTER
In the text
thumbnail Fig. 14

Distribution of the emissivity index β (right) and its relative uncertainty δβ/β (left) along a vertical cut through the center. The red and blue lines indicate the center of M 33 and R ≃ 6 kpc, respectively.
Open with DEXTER
In the text
thumbnail Fig. 15

Scatter plots of the dust emissivity index β versus the Hα, CO(1 − 0), and HI emission from M 33. The y-axis presents the logarithm of the intensities in units of cm-6 pc for the Hα emission and K km s-1 for the CO(1 − 0) and HI emssions. Also shown are Pearson correlation coefficients (r), the ordinary least-squares (OLS) fits, and the bisector fits.
Open with DEXTER
In the text
thumbnail Fig. 16

Difference in the dust emissivity index β derived using the two-component MBB and the single-component MBB models vs. radius.
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.