EDP Sciences
Free Access
Issue
A&A
Volume 599, March 2017
Article Number A64
Number of page(s) 22
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201629251
Published online 01 March 2017

© ESO, 2017

1. Introduction

Starlight in galaxies is processed by dust grains through absorption and scattering. On average, about one third of starlight in normal spiral galaxies is attenuated by dust (Popescu et al. 2002; Skibba et al. 2011; Viaene et al. 2016). Dust resides around new stars in their birth clouds as well as in the diffuse interstellar medium (ISM). A Milky Way-like dust mixture is highly efficient in scattering and absorbing UV radiation, but also significant fractions of optical and sometimes even near-IR (NIR) light. All the energy that is absorbed by the dust is reprocessed and emitted at longer wavelengths.

In order to understand the processes governing the observed dust-starlight interplay, observations at all wavelengths are necessary. This allows the investigation of how the UV and optical energy from stars is converted into far-IR and sub-millimeter (FIR/submm) emission. Such knowledge is important since FIR/submm emission is a useful tool to study galaxies both near and far; the amount of dust obscuration roughly follows the star formation rate density in cosmic time (see e.g. Daddi et al. 2005; Gruppioni et al. 2013; Madau & Dickinson 2014). A correct understanding of dust heating allows the separation of dust into multiple temperature components. This in turn results in better dust mass estimates. Accurate dust masses can then be used to estimate the total ISM content in distant objects for which Hi content cannot be observed (see e.g. Corbelli et al. 2012; Eales et al. 2012; Scoville et al. 2014; Groves et al. 2015).

Nearby galaxies are good test laboratories for such dust heating studies. Their dust heating mechanisms have been studied through various techniques such as correlations with star formation rate (SFR) indicators and stellar mass tracers (Galametz et al. 2010; Boquien et al. 2011; Bendo et al. 2012; Foyle et al. 2013; Hughes et al. 2014; Bendo et al. 2015) or through panchromatic energy balance SED fitting (Groves et al. 2012; Aniano et al. 2012; Dale et al. 2012; Mentuch Cooper et al. 2012; Ciesla et al. 2014; Rémy-Ruyer et al. 2015; Boquien et al. 2016). While the situation differs from galaxy to galaxy, the current understanding is that the cold dust in most galaxies is heated by the evolved stellar populations which feed the general interstellar radiation field (ISRF). The warmer dust is then predominantly heated by the young and new stars, through strong UV radiation. In certain cases, even active galactic nuclei (AGN) can contribute to the dust heating (see e.g. Wu et al. 2007; Kirkpatrick et al. 2012; Schneider et al. 2015).

Usually, these investigations provide a qualitative view on the dominant heating source. In starbursts and galaxies with active star formation, dust is predominantly heated by young stellar populations. For early-type galaxies, little research has been conducted to estimate the dust heating sources. Based on their low star formation rates, one would expect that the evolved populations are dominant in these systems (see e.g. Smith et al. 2012b; Boquien et al. 2014). But other heating sources are possible, for example X-rays from hot halo gas (Natale et al. 2010). Most galaxies, however, fall in the ambiguous intermediate regime, where both evolved and young stellar populations contribute to the dust heating (Bendo et al. 2015). It is dangerous to infer properties like dust mass and SFR from FIR/submm data alone in these systems.

A powerful method to investigate dust heating processes is through dust continuum radiative transfer (RT) simulations, which mimic the transport of radiation through a dusty medium. For a review of this method, we refer the reader to Steinacker et al. (2013). Such simulations treat the dust-starlight interaction in 3D, allowing realistic geometries, multiple viewing angles and non-local heating. But it also comes with a large computational cost and the need for “simple”, well-behaved targets to model. Therefore, the most common application is in models of edge-on spiral galaxies (see e.g. Popescu et al. 2000; Misiriotis et al. 2001; Bianchi 2008; Baes et al. 2010; Popescu et al. 2011; MacLachlan et al. 2011; De Looze et al. 2012a; De Geyter et al. 2015) or early-type galaxies (De Looze et al. 2012b; Viaene et al. 2015).

Through consistent treatment of primary emission, scattering, absorption and thermal re-emission, dust heating mechanisms can be investigated. A general result of these studies is that starlight absorption by diffuse dust usually fails to heat this dust to sufficient levels to match the observed FIR/submm emission. Often, obscured star formation is necessary to balance the energy between UV/optical (absorption) and FIR/submm (emission).

In a recent attempt to quantify the dust heating mechanisms, De Looze et al. (2014) performed a detailed radiative transfer simulation of M 51. Based on observed images, they closely mimic the distribution of both stars and dust, and then perform an accurate treatment of the dust-starlight interaction. They find that the contribution of evolved stars to the thermal re-emission is significant and varies with location and observed IR wavelength. In the mid-IR (MIR), the evolved stellar populations contribute about 10% to the total dust heating from primary stellar emission, while young stars consequently contribute 90%. This changes at submm wavelengths, where the contribution of evolved stars is ~40%. There is also a significant enhancement of the contribution by young stars in the spiral arms (>80%) with respect to the inter-arm regions (<60%).

In this paper, we apply this novel technique to the Andromeda galaxy (M 31). It is the most massive object in our Local Group and close to our own Milky Way. In the optical, the galaxy is dominated by its bright bulge (Tempel et al. 2010). The smooth disk is intersected by several dark lanes where dust obscures the starlight. The dark dust patches coincide with bright emission in the FIR and submm (Fritz et al. 2012; Groves et al. 2012). On the other side of the spectrum, in the UV, a striking similarity with the FIR morphology is observed (Thilker et al. 2005). The morphological (anti-)correlations are a manifestation of intricate interactions between dust and starlight. It is also evident that the observed features point at a complex 3D structure. Looking at M 31 as it is projected on the sky inevitably leads to parameters which are summed or averaged along the line of sight. A 3D model of the sources and sinks for radiation in M 31 can help to understand our closest neighbour beyond the standard line-of-sight quantities. It grants us the opportunity to investigate the morphology from different viewing angles and its influence on the spectral energy distribution (SED), the attenuation law and the processes of scattering, absorption and re-emission by dust.

This endeavour sets new challenges to the radiative transfer modelling as M 31 is much larger and has a higher inclination angle than M 51 (in this work, we will use i = 77.5°, which is the average value for the disk of M 31 as derived from Hi data by Corbelli et al. 2010). On the other hand, M 31 is an obvious choice for studying dust heating mechanisms in great detail. Even for the hardest accessible wavelengths, the spatial resolution is still <140 pc along the major axis. Dust heating mechanisms in M 31 have been investigated previously using UV-FIR energy balance SED fitting (Montalto et al. 2009; Groves et al. 2012; Viaene et al. 2014), using FIR/submm SED fitting by Smith et al. (2012a) and Draine et al. (2014), and through FIR colours (Planck Collaboration Int. XXV 2015). All of these methods show that the emission from evolved stellar populations are the main contributor to the dust heating, especially at wavelengths beyond 160 μm. However, quantifying the contributions of the different heating sources and their wavelength dependency remains difficult. Our study approaches this problem from a new perspective (radiative transfer simulations), and will provide quantitative measures related to dust heating in this galaxy.

In Sect. 2 of this paper we describe the panchromatic dataset we use and core data products we derive from it. Our model and all its components are outlined in Sect. 3 and validated in Sect. 4. We present a 3D view of Andromeda in Sect. 5 and a dust heating analysis in Sect. 6. We present our main conclusions in Sect. 7.

2. A panchromatic dataset of M 31

The data used in this work is the same as for our previous study (Viaene et al. 2014, hereafter HELGA IV). An elaborate description of the data can be found in that paper, we give a brief overview here. The panchromatic dataset is a combination of dedicated observing campaigns using ground based and space telescopes. In the UV domain, M 31 was observed by GALEX (Martin et al. 2005) in the far-UV (FUV) and near-UV (NUV) bands (Thilker et al. 2005). Optical data were acquired from SDSS (York et al. 2000) and mosaicked by Tempel et al. (2011), taking special care of the background variations. The result is a large (2.5° × 8°) field of view in the ugriz bands. In the near and mid-infrared, we rely on observations from Spitzer (Werner et al. 2004) using the IRAC instrument (Fazio et al. 2004) and WISE (Wright et al. 2010). Andromeda was imaged in all four IRAC bands and is described in Barmby et al. (2006). High-quality mosaics of the WISE observations were created in all four bands by the WISE Nearby Galaxy Atlas team (Jarrett et al. 2013). In the FIR regime, we use observations from the MIPS (Rieke et al. 2004) instrument aboard Spitzer. Due to resolution and sensitivity restrictions, we only use the maps at 24 and 70 μm from Gordon et al. (2006).

Most recently, we completed our dataset with FIR/submm imaging using the Herschel Space Observatory (Pilbratt et al. 2010). The Herschel Exploitation of Local Galaxy Andromeda (HELGA) observed the galaxy with PACS (Poglitsch et al. 2010) at 100 and 160 μm and with SPIRE (Griffin et al. 2010) at 250, 350 and 500 μm. The PACS data were reduced using standard recipe with HIPE v12 and SCANAMORPHOS v24 (Roussel 2013). The SPIRE data using HIPE v12 and a custom method for the baseline subtraction (BriGAdE, Smith et al., in prep.). An elaborate account on the data acquisition and reduction can be found in Fritz et al. (2012, HELGA I) and Smith et al. (2012a, HELGA II). The Herschel data for M 31 is now publicly available online1 in its newly reduced form.

In HELGA IV, we have performed a series of data manipulations in order to bring this diverse dataset to a consistent set of images. First, foreground stars were masked based on UV and NIR colours. The masked areas were replaced by the local background, which often was the diffuse emission from M 31 itself. Second, all images were convolved to match the SPIRE 500 μm beam, which is the largest in our set. Finally, the frames were rebinned and regridded to match this beam size. The end products are a series of images with a pixel size of 36 arcsec (or 140 pc along the major axis of Andromeda) and with a beam of that same size. Throughout the paper, we use a distance to M 31 of 785 kpc (McConnachie et al. 2005).

We work within a limited field of view to ensure good coverage at all wavelengths. The studied region is limited by an ellipse with centre at α = 00:43:06.28 and δ = +41:21:12.22, semi-major axis of 1.425°, semi-minor axis of 0.400° and a position angle of 38.1°. The corresponding physical scale is about 19.5 kpc along the major axis. This corresponds to the B-band D25 radius (de Vaucouleurs et al. 1991) and encloses the star-forming ring and the brightest outer dust features. A few faint substructures beyond 25 kpc are not included in this field-of-view. This should not affect the dust heating analysis, which is limited to the inner 20 kpc of the galaxy.

Although foreground stars and background galaxies are masked in the images by HELGA IV, M 32 and NGC 205 are not. Both dwarf galaxies are satellites of Andromeda and M 32 in particular lies in front of our area of interest. This object is thought to have played an important role in the recent history of Andromeda due to a head-on collision (see e.g. Block et al. 2006; Gordon et al. 2006; Dierickx et al. 2014). M 32 is currently not part of the disk of Andromeda and its emission will influence our measurements of the main galaxy itself. We therefore mask out this dwarf companion where it is detected (from the FUV to the 24 μm band). To do that, we fit a 2nd order 2D polynomial to the pixels surrounding M 32 and replace the emission from M 32 by the interpolated local background.

Finally, we use integrated fluxes from IRAS (Neugebauer et al. 1984) and Planck (Planck Collaboration I 2011) as additional constraints on the global SED of M 31. IRAS fluxes at 12, 25, 60 and 100 μm were taken from Miville-Deschênes & Lagache (2005). We use the Planck 350, 550 and 850 μm fluxes derived by Planck Collaboration Int. XXV (2015). Most of the images described above will serve as observational constraints on the model. However, we use a handful of images as an input for our model construction. They are assumed to trace different components of M 31 and are discussed in Sect. 3.2.

3. Constructing a 3D model

3.1. Radiative transfer simulations

We make use of SKIRT2 (Baes et al. 2011; Camps & Baes 2015), an advanced dust continuum RT code that uses the Monte Carlo approach. The code allows panchromatic RT simulations using efficient dust grids (Camps et al. 2013; Saftly et al. 2014), a vast suite of possible geometries and geometry decorators (Baes & Camps 2015) and was recently updated with a user-friendly interface. It supports several dust models, and a full treatment of stochastically heated dust grains (Camps et al. 2015).

A key feature for this investigation is the input of a 2D FITS image as a possible geometry, first demonstrated by De Looze et al. (2014). It uses the relative surface brightness in the image pixels to set up a 2D density distribution in the virtual 3D space. This density distribution can then be deprojected using a simple cos(i) factor, with i the inclination angle. To give the deprojected 2D geometry a 3D nature, the density is smeared out in the vertical direction according to an exponential profile with a vertical scale height. This powerful technique allows more realistic stellar geometries such as asymmetric features or clumpy and filamentary structures. For M 31 in particular, with its highly disturbed and clumpy rings (HELGA I), this is a big step forward. Of course, the deprojection of a highly inclined galaxy such as M 31 (i = 77.5°) is always degenerate. Brighter spots in the disk will be smeared out in the direction of deprojection and the assumption of an exponential vertical profile will again smear the bright spots in the vertical direction. This technique does, by construction, conserve flux during the conversion from 2D to 3D.

3.2. Model set-up

The radiative transfer model of Andromeda contains several stellar components and two dust components. The details of our model setup are outlined in Appendix A. We give a brief overview here. The components are summarised in Table 1.

Table 1

Overview of the different stellar populations (SP) and dust component in our model.

To keep the balance between a realistic model and computational feasibility, the number of components has to be limited. We split the stars in M 31 into four components: an evolved bulge population (12 Gyr), an evolved disk population (8 Gyr), a young disk population (100 Myr), and an population of ionising stars (10 Myr). The dust in M 31 is split in an interstellar component, and a small fraction of dust around the ionising stellar component. All stellar components have two aspects: geometrical distribution in 3D space, and a panchromatic SED.

The geometrical distributions are based on observed images. For the evolved populations, we take an analytical 3D Einasto (1965) profile for the bulge. This bulge is subtracted from the IRAC 3.6 μm image to obtain the disk as seen on the sky. To make the disk 3-dimensional, it is deprojected, and a vertical dimension is added in the form of an exponential surface brightness profile. The young stellar component is based on the FUV image of M 31, with a flux correction for UV emission from evolved stars. The ionising stellar component is based on the MIPS 24 μm image, also corrected for flux originating from evolved stars. The latter two observed images are also deprojected and a vertical exponential profile is added to obtain a 3D geometry. However, the vertical scale height of the young and ionising component is set to be much smaller than for the evolved stellar disk. SEDs are assigned to the distributions, and normalised to the flux in a single waveband, see Table 1.

The dust around stars in their birth clouds is embedded as a subgrid recipe in the SED template for the ionising stellar population. It thus follows the distribution of this stellar population. The distribution for the interstellar dust is derived from the dust mass map from Viaene et al. (2014). This map is again deprojected and a vertical exponential profile is added to create a 3D distribution. The optical properties of the dust are set according to the THEMIS dust mix model (Jones et al. 2013).

The dust grid in our simulations is constructed using a hierarchical binary tree. The cells are built by splitting the parent cell in turns in the x, y and z direction. Cells are not further split when their dust mass fraction is below 10-6. For more information about this dust grid, we refer the reader to Saftly et al. (2014). The grid is automatically constructed and has 1.45 million dust cells for the full M 31 model. Note that these cells are not equal in size. Their volume ranges between (305000 pc)3, with a median value of (89 pc)3.

3.3. Model optimization

Table 2 gives an overview of the main parameters that have to be set in the simulation. For a detailed account on their physical meaning and the initial guess values, we refer the reader to Appendix A.

In order to find a set of parameters that provides a good working model, we fix all but three parameters: the intrinsic FUV luminosities of the young and ionising stellar components, and the total dust mass. For computational reasons, we can only explore a small phase space. We therefore make our initial guess values as accurate as possible (see Appendix A). For the dust mass, we search the parameter space within a ± 25% range around our initial guess. The intrinsic FUV radiation of the young component ranges from 0.75−2.0 times the initial guess value. The intrinsic FUV luminosities for the ionising component ranges from 0.4 to 2.6 × 108L.

Table 2

Overview of the parameters in our model.

Using this parameter grid, we run a panchromatic RT simulation for each of the 1050 possible combinations of the free parameters. We compute the stellar emission and the effects of dust at 128 wavelengths between 0.05 and 1000μm. The wavelength grid is logarithmically spaced over this range, but has a finer spacing in the MIR range (130 μm) to capture the aromatic features. The best fitting model SED is determined by the lowest χ2 between observed and model fluxes. The model fluxes were obtained by multiplying the model SED with the filter response curves and integrating the flux. The MIR dominates the χ2 when using a uniform weighing because there are many data points in this regime. However, to study the dust heating in M 31, it is more important to adequately fit the UV attenuation and the dust emission peak. To ensure this, we lower the weights of the MIR data points to give priority to a good fit in the FUV and FIR/submm.

thumbnail Fig. 1

Panchromatic SED of Andromeda. The black line is the best fitting RT model, run at high resolution (512 wavelengths). The green points are the observed integrated luminosities for M 31. The red, blue and cyan lines represent the SEDs for simulations with only one stellar component: evolved, young and ionising, respectively. The interstellar dust component is still present in these simulations. The bottom panel shows the residuals between the observations and the best fitting (black) model.

Open with DEXTER

The best fitting model is re-simulated at higher resolution to extract images of higher signal-to-noise (S/N). The wavelength sampling was increased by a factor 4 to 512 wavelengths between 0.05 and 1000 μm. Additionally, we use 20 times more photon packages (20 × 106 per wavelength in total) to accurately sample extinction, scattering, and emission. The SED of the high-resolution model is shown in Fig. 1, together with the separate stellar components. To this end, we re-run the simulation each time using only one of the three stellar components (and still including the dust). Note that dust emission depends in a non-linear way on the absorbed energy, so the individual SEDs do not add up to the total SED. However, we can already learn some interesting trends on the energy balance of the galaxy.

The total SED of M 31 is clearly dominated by the evolved stellar populations. They produce nearly all of the optical and NIR light, and are the strongest source in the MIR. They are also indirectly responsible for the bulk of the dust emission. We come back to this in Sect. 6. Only in the UV domain, the contribution of the young (non-ionising) stellar component is higher than that of the evolved stellar populations. The emission of the ionising stellar populations is comparable to the evolved stellar populations in this domain. However, this is still about 5 times lower than the young, non-ionising component.

The SED of the ionising component has two emission peaks in the FIR/submm. This points to two very distinct dust components. The emission peak round 40 μm is associated with warmer dust, located inside the molecular clouds and heated by the ionising stars. The second peak in the SED of the ionising component occurs around 250 μm, indicating colder dust. This is emission from dust in the diffuse ISM, heated by radiation escaping the molecular clouds.

The emission from dust in the diffuse ISM peaks at different wavelengths in the SEDs of the different components. If only evolved stellar populations are included, the peak occurs at a much shorter wavelength (~150 μm) than the peaks for the young (~210 μm) and ionising (~250 μm) component. This indicates that the observed average dust temperature of Andromeda is mostly influenced by its evolved stellar populations. We perform a more quantitative analysis of the dust heating in Sect. 6.2.

4. Model validation

4.1. Global SED

Before investigating the 3D structure of M 31 and the dust heating mechanisms in this galaxy, we perform a series of quality checks on the best fitting model. It has a dust mass of Mdust = 4.7 × 107M. The intrinsic FUV luminosity of the young and ionising component is and , respectively. This model SED is shown in the top panel of Fig. 1, together with the observed fluxes. While the correspondence is generally good, there is still room for higher MIR/FIR luminosities to fit the observations in this area. However, that requires more FUV emission and/or more dust mass, which will affect the UV part of the SED. Since we tried these options in our modelling routine and did not get a lower χ2, we assume this situation is not preferable.

The dust mass is on the high end of previous dust mass estimates (HELGA I; HELGA II; HELGA IV; Draine et al. 2014), but still consistent. One must note that a different dust model is assumed here. Despite comparable optical properties (tuned to fit the Milky Way extinction curve), the composition and grain density can differ significantly. This can lead to a difference of a factor 23 (Jones et al. 2013; Rémy-Ruyer et al. 2015; Dalcanton et al. 2015; Planck Collaboration Int. I 2016).

Our other two free parameters will have a significant influence on the SFR (both from UV or FIR measurements). We therefore “observe” our model SED and apply several recipes from Kennicutt & Evans (2012) to convert flux to a SFR. The resulting values for the SFR are summarised in Table 3. For the FUV and NUV estimates, the intrinsic UV fluxes were used to avoid the problem of dust attenuation. The SFR values derived from single band fluxes lie close to each other. They are well in line with previous estimates of the SFR in M 31, which lie between 0.2 and 0.4 M yr-1 (see e.g. Tabatabaei & Berkhuijsen 2010; HELGA IV; Rahmani et al. 2016; or Ford et al. 2013 (HELGA III)).

A caveat to using only 24 μm or 70 μm flux is that it assumes all UV radiation from star formation is absorbed by dust. This is not the case for M 31, so these bands naturally underestimate the SFR. We therefore also apply the hybrid SFR estimator from Leroy et al. (2008), which combines the observed FUV with the observed 24 μm flux. This recipe captures both the obscured and unobscured star formation. We find a value of 0.31 M yr-1, which is significantly higher than the single band FIR tracers. However, it is still of the same order of magnitude and consistent with literature values for the SFR in Andromeda.

Offset from these values is the SFR derived from the total infrared light (TIR). This calibration gives a SFR rate which is roughly three times higher than the other estimates. This is most likely because the dust heating in M 31 is not dominated by the young stellar populations, as we will show in Sect. 6. TIR is therefore not a reliable tracer of SFR in Andromeda. In fact, even for galaxies with more star formation than M 31, TIR can lead to an overestimated SFR (Boquien et al. 2014).

Table 3

Star formation rates for M 31, derived from synthetic broad band fluxes from our model SED.

To estimate the uncertainty on the free parameters, we construct probability distribution functions (PDFs). For each of the tested models, we compute the probability as proportional to exp(−χ2/ 2). For each parameter value we then add the probabilities and plot the PDFs in Fig. 2. We find that dust mass and are fairly well constrained (the PDFs resemble a normal distribution). The , on the other hand, has a more or less flat probability distribution and thus is poorly constrained. This is not so surprising since this parameter corresponds to constraining of the cyan curve in Fig. 1. This normalization is fitted in a parameter range far below the luminosity of the old and young stellar components, making it difficult to find the optimal solution. The median (50th percentile) values that come from these PDFs are , and . Here the upper and lower errors correspond to the 84th and 16th percentile of the PDFs, respectively.

thumbnail Fig. 2

PDFs of the three free parameters in our model optimization: the total dust mass Mdust (left), the intrinsic FUV luminosity from the young stellar component (middle) and the intrinsic FUV luminosity from the ionising stellar component (right). Black dashed lines are the parameter values for the best fitting model. Red dashed lines are the median (50th percentile) values from the PDFs.

Open with DEXTER

4.2. Attenuation law

We can accurately reconstruct the attenuation law in our model since we know the sources and sinks of the starlight. The attenuation curve for the best fit model is given in Fig. 3, normalised to the V-band attenuation. We find a steadily increasing attenuation with decreasing wavelength, with on top of that a broad bump that peaks around 0.22 μm. This curve is a combination of attenuation by diffuse dust as modelled by the dust mass map described in Sect. A.4, and attenuation by dust in star-forming clouds. The latter is incorporated in the MAPPINGS III SED templates (Groves et al. 2008), which were used for the ionising stellar component (see also Sect. A.3). The attenuation laws for both components are also shown in Fig. 3. Note that the MAPPINGS attenuation law is directly provided by Groves et al. (2008), while the diffuse dust attenuation curve is a combination of the THEMIS extinction law and the relative star-dust geometry.

The dust around new, ionising stars contributes less than 0.1% of the absorbed energy in the V band. However, this ratio increases to 5% in the UV bands. The total attenuation curve is thus mainly shaped by the diffuse dust. The effect of the MAPPINGS III attenuation curve is visible in a slightly higher NUV bump, and in a sharper increase at the shortest wavelengths.

In Fig. 3, we also compare against several literature attenuation curves: (a) the average Milky Way extinction curve from Fitzpatrick & Massa (2007), converted to an attenuation curve using RV = 3.001 (see also their Eq. (10); (b) the SMC bar region from Gordon et al. (2003), again this is an extinction curve, converted to an idealised attenuation curve; (c) the average attenuation law for 10 000 local (z ≲ 0.1) star-forming galaxies from Battisti et al. (2016); (d) the average attenuation measurements from Dong et al. (2014) for the central 1 of M 31.

thumbnail Fig. 3

Attenuation laws from our model, normalised to the V-band attenuation. The black solid line is the global law. The dashed line reflects the attenuation law for the diffuse dust and the dotted black line for the dust in star-forming regions. Several literature measurements are shown; red curve: MW curve for R(V) = 3 from Fitzpatrick & Massa (2007), blue dashed line: SMC bar region (Gordon et al. 2003), and green dash-dotted line: attenuation curve from Battisti et al. (2016).

Open with DEXTER

thumbnail Fig. 4

Comparison of the model with observations in selected wavebands. First column show the observed images, second column the model images. The 3rd column shows the residuals derived as in Eq. (1). From top to bottom row: FUV, g, 3.6 μm, 8 μm, 24 μm, 100 μm and 250 μm.

Open with DEXTER

At NIR and optical wavelengths (λ> 0.4 μm) our attenuation curve follows the others almost perfectly. At shorter wavelengths, the curves start to diverge. The UV bump in our model curve is much broader than that of the MW. The broadness of the bump of this dust model is discussed in Jones et al. (2013). It relates to the grain size distributions and their band gap properties. The discrepancy is amplified because the MW curve is an idealised attenuation curve, converted from an extinction law. This does not fully incorporate the large-scale spatial distribution of stars and dust. An additional complication in comparing attenuation curves is that their shape can depend on the resolution of the study (Boquien et al. 2015).

On the long wavelength side of the bump, our model is consistent with the empirical measurements of Dong et al. (2014). Their UV bump does seem to be less pronounced. Unfortunately, they have no measurements shortward of 0.193 μm. The broad UV bump can explain why our model has difficulties to reproduce the observed NUV flux (see Fig. 1).

Beyond the UV bump, our model curve steadily increases and falls somewhere between the MW and SMC idealised attenuation laws. In general, our model produces an attenuation law that is broadly in line with observations, indicating that our treatment of extinction and scattering yields realistic results.

4.3. Images vs. observations

We now compare model and observations at a spatially resolved scale. Figure 4 shows a selection of important bands across the spectrum. We only show pixels with an observational S/N above 2. The easiest way to visualise the difference between observation (left column) and model (middle column) is through the residual image (right column): (1)A negative (positive) fraction means the model overestimates (underestimates) the observed flux.

The median (absolute value) deviation between model and observations across all bands is 22%. Most of the model pixels fall within 50% of their observed counterpart. Some bands (g, 3.6 μm and 250 μm) have much smoother residuals, while others (FUV, 8 μm, 24 μm, 100 μm) display clear structure. The duality appears to be related to the appearance of M 31 across wavelengths. In the NIR, Andromeda is a smooth galaxy without sharp features. The residuals in this regime are particularly regular. Bands in the MIR regime, with contributions of evolved stars, hot dust and aromatic features, show a very clumpy and ring-like Andromeda. The residuals in these bands exhibit strong features.

In general, the model tends to underestimate the flux in the rings of M 31 and overestimate the flux in the diffuse inter-ring regions by about 20% on average. A first explanation for this can be the vertical dimension that was added to the model. Upon deprojecting the images, bright regions tend to be smeared out in the direction of deprojection. After that, the light is again smeared out in the vertical direction. In Appendix C, the effect of alternative scale heights for the disk components is investigated. However, we find that the effect is minor and cannot explain the bulk of the residual features.

A second aspect that could be responsible for at least a fraction of the discrepancy between observations and model is the limited resolution, and in particular the fact that we are not capable of modelling the star-forming regions in detail. Due to their intrinsic spherical shell geometry, these models emit isotropically and have a higher average attenuation per unit dust mass than models with a more asymmetric or clumpy distribution (Witt & Gordon 1996, 2000; Városi & Dwek 1999; Indebetouw et al. 2006; Whelan et al. 2011). This effect is particularly relevant in the UV, where the ionising stellar populations are relatively luminous and attenuation effects are pronounced.

As we discuss in the construction of our model (see Appendix A.2), we probably underestimate the contribution of young stellar populations in the dustiest regions of the galaxy. Vice versa, we overestimate their contribution in the least attenuated regions. This could be another cause for the discrepancy observed in the residual maps in Fig. 4. In the model, the radiation field in the dustiest regions therefore does not hold enough UV energy. As a result, the model underestimates the observed flux in these regions, and overestimates the flux in the more dust-free areas.

There may be a fourth effect contributing to the discrepancy. De Looze et al. (2014) also observed residual features that could be associated with the spiral arms in their M 51 model. They argue that there is a difference in both stellar populations and dust properties between arm and inter-arm regions. In HELGA II it seemed that such variations caused differences in the dust emissivity index β in their pixel-by-pixel modified black body fits of M 31. However, a variable β appears to correlate with 3.6 μm emission (Kirkpatrick et al. 2014). This suggests that variations in beta may not actually represent variations in dust emissivity but instead may represent a variation in dust heating sources. On the other hand, Galametz et al. (2012) found that allowing β to vary can lead to unphysical dust temperatures in pixel-by-pixel energy balance SED fitting.

Nevertheless, it is not unlikely that spiral arms host different dust mixtures than inter-arm regions. While we allow for 3 stellar components in our model, we take the first order assumption of just one type of dust mixture for the entire galaxy. Grain size distributions may vary significantly depending on radiation field and ISM phase (see e.g., Ysard et al. 2013, 2015). However, it would lead us too far to investigate the changing size distributions and grain properties in the different regions of M 31. Furthermore, this would multiply the number of parameters in our model. We therefore choose to stay at the current level of complexity, while being aware of the caveats.

4.4. FIR colours

As a final check for our model, we compare subsequent FIR/submm colours with observed ones in Fig. 5. Note that we define colours as logarithmic flux ratios. We only considered pixels with S/N> 3 as the PACS maps can have large uncertainties at low surface brightness. There is generally a good correspondence between model and observation. This was already evident from the global SED in Fig. 1, where the cold dust peak is adequately reproduced. The colour maps, however, show that this is also the case on the resolved scale.

The observed 70/100 and 100/160 colours are more “clumpy” than in the model. This is mostly due to the low S/N in the MIPS and PACS observations. On average, however, the correspondence is satisfactory, although the observed maps are slightly redder. We chose not to use these colours for further analysis because of the higher uncertainties that come with it.

Of particular interest are the 160/250 and 250/350 surface brightness ratios, which cover the wavelength range where the transition of dust heating sources occurs. Consequently, they have been used in the past to identify the dominant dust heating sources (see e.g. Bendo et al. 2010; Boquien et al. 2011; Boselli et al. 2010, 2012; Bendo et al. 2012, 2015; Planck Collaboration Int. XXV 2015, and references therein). There is a reasonable agreement in the morphology, with blue colours in the centre, and redder colours in the outskirts. The model maps are again smoother than the observed colour maps, but the difference is less pronounced than at shorter wavelengths. In Fig. 1, we find that our model falls below the observation at 160 μm, while it agrees with the 250 μm flux. The model thus underestimates the 160/250 colour on the global scale. The resolved colour map indicates that this underestimation comes from the star-forming ring. The colours are bluer in the observed image. This is partially compensated by the centre and outskirts, where the model produces slightly redder colours than observed. The 250/350 colour is also underestimated and the offset is more general. The observed colour is bluer than the model in most of the pixels.

For completeness, we also present the 350/500 colour. This is often neglected on a local scale because of the low angular resolution that comes with including the SPIRE 500 μm band. Additionally, both the 350 μm and 500 μm band trace the Rayleigh-Jeans tail of the cold dust emission. It is therefore less informative than colours at shorter wavelengths. Again, observations and model are well in line with each other. The main deviations are the slightly redder outer regions in the observed colour map, but the effect is minor.

Taken together, the observed colours are fairly well reproduced by our model. This points at an adequate treatment of the dust physics, which is reassuring.

thumbnail Fig. 5

Observed (left) and model (right) colours tracing the FIR dust emission peak. A colour here corresponds to the flux ratio: log (flux1/ flux2), that is, blue colours indicate a high flux ratio, red colours a low ratio.

Open with DEXTER

5. A face-on view of M 31

The 3D RT simulations of Andromeda allow us to view the galaxy from different angles. Most informative is the face-on view shown in Fig. 6 in six bands. The images highlight the diversity and complexity of the galaxy across the spectrum. In the FUV band, the clumpy star formation regions are structured in several rings and arcs. The bulge is more compact than in the optical, but still the brightest region. Most of the bright spots are elongated in the direction of the deprojection. This effect is inherent to the method and cannot be corrected for.

At 3.6 μm the galaxy has a smooth disk with only the ring at 10 kpc as the reminder of the complex ring-arc morphology visible at 0.153 μm. The bulge is large and dominates the NIR emission. A small gap can be observed in the middle of the image, along the vertical direction. This feature is an artefact of the bulge subtraction from the disk. Because we approximate the observed boxy bulge with an Einasto profile, some residual sidelobes remain in the disk. These sidelobes are smeared out in the direction of deprojection. This creates the two brighter features on the upper right and lower left of the vertical “gap” in the model 3.6 μm image.

The bulge is still the brightest feature at 8 μm. Through the FIR and submm (24, 100 and 250 μm images) the look of Andromeda gradually evolves. The bulge becomes less dominant and the inner and outer rings turn brighter. The rings also exhibit a clear broadening going toward the submm as the extent of the diffuse dust becomes apparent at these wavelengths. While the bulge is still visible at 100 μm, there is no sign of it at 250 μm and M 31 becomes one swirl of clumpy rings.

thumbnail Fig. 6

Face-on model images of M 31 in the FUV, 3.6, 8.0, 24, 100 and 250 μm. Colour scale is dif ferent for each image for clarity. The green circles on the 250 μm image (bottom right) correspond to radii of 10 kpc and 15 kpc from the centre.

Open with DEXTER

6. Dust heating mechanisms

We now aim to identify the dust heating sources and quantify their respective contribution. We will first look at the intrinsic 3D radiation field and link that to the heating of dust (Sect. 6.1). Then, we define and compute the dust heating parameters for the projected model on the sky (Sect. 6.2). This way, we can link them to observational properties (Sect. 6.3).

6.1. 3D analysis

We can consider each dust cell in the model as a voxel (i.e. a 3D pixel) that processes stellar radiation. During the simulation, the absorbed luminosity is stored in every cell and at every wavelength. Integrating over wavelength gives the total absorbed luminosity Labs in each cell.

Such an approach does require a good sampling of the radiation field in each dust cell during the simulation. If only a few photon packages pass through a dust cell, the Monte Carlo noise will be high. This issue is not relevant when looking at the line-of-sight projected pixels, since many photon packages travel through the line-of-sight and provide sufficient sampling. To investigate the absorption for the dust cells, we run a separate set of simulations with more photon packages (108 per wavelength). The number of wavelengths is limited to 25, sampled logarithmically between 0.1 and 10 μm, to compensate the computational cost that comes with shooting more photon packages. This wavelength range is adequate since the absorbed radiation mainly originates from stellar sources. We integrate the absorbed luminosity along the wavelength axis to compute Labs. The coarser sampling of the wavelength grid is compensated by the gain in S/N due to the integration. This method is similar to the formalism outlined in Natale et al. (2015), Eqs. (3)(5). They start from the stellar radiation field to compute the total absorbed luminosity, whereas we start from the absorbed energy per cell. In a well-sampled simulation the two formalisms are equivalent.

We perform this simulation 4 times: including all stellar components yields . Including only evolved, young or ionising stellar populations yields , and , respectively. is the luminosity that originates from the MAPPINGS template SED, and is absorbed by the diffuse dust component. There is also an internally absorbed luminosity in the subgrid implementation of the MAPPINGS model for each dust cell. This can be seen as the absorbed energy from ionising stars by dust in star-forming clouds (as opposed to energy absorption due to diffuse ISM dust, ).

To test the sampling of the radiation field in our simulations, we plot the sum of , and against the total absorption by the diffuse dust in Fig. 7. It is not necessary to count the internal absorption in the star-forming clouds here because the energy balance is internally secured in the subgrid model. The 1:1 line is plotted in red and shows a good correspondence between the three quantities. Only at low absorption luminosities, the variance starts to increase. This indicates that the radiation field in each dust cell is well sampled and we can separate the contribution from evolved, young and ionising stellar populations to the total absorbed luminosity. We then define the contribution of young and ionising stellar populations to the absorbed energy per dust cell as (2)where we consider the contributions of the young and ionising stellar populations together under the umbrella of “unevolved stellar populations”. Consequently, the absorbed luminosity fraction from the evolved stellar populations is .

thumbnail Fig. 7

Left: 2D histogram comparing the integrated absorbed luminosity per dust cell including all stellar components with the same quantity, but as the sum of separate simulations including only evolved, young and ionising stellar populations. Middle: distribution of the bolometric absorption fraction of unevolved stellar emission per dust cell: . The histogram is weighted by mass fraction of each dust cell and normalised. The red line indicates the median . Right: 2D histogram showing the radial distribution of . The bins are weighted by the mass fraction of the dust cells. Red indicates a high number of data points, blue a low number. The spread at a fixed radius is a combination of dust cells at the same radius and with different vertical location in the disk. The white line is the weighted mean heating fraction per radius.

Open with DEXTER

First, we take a global measure of the dust heating by summing, for each stellar component, the absorbed luminosities over the dust cells. Equation (2) can then be applied to these global absorption luminosities to yield the global absorption fraction for the galaxy. We find that on a global scale for M 31. In other words, 91% of the energy absorbed by dust originates from the evolved stellar populations. This is a particularly high number (De Looze et al. 2014 found 37% for M 51 using a comparable technique) and most likely dominated by strong radiation field of the bulge.

Second, we take a closer look at the dust heating of the individual dust cells in the middle and right panel of Fig. 7. The middle panel shows the distribution of for the dust cells. We weighted each dust cell with its mass fraction to account for the difference in size and dust content between the individual dust cells. The distribution peaks near the median value of 0.10 and is skewed towards higher fractions. However, the contribution of unevolved stellar populations to the absorbed energy hardly exceeds 50%. This indicates that, even at the level of individual dust cells, the heating by evolved stellar populations dominates.

The radial variation of the dust heating fraction can be investigated by plotting as a function of galactocentric radius for each dust cell. In the right panel of Fig. 7, this is shown for each dust cell, irrespective of its vertical coordinate. The 2D histogram is weighted by the mass fraction, and the weighted mean for each 1 kpc bin is shown as a white line. The contribution of the unevolved stellar populations is only a few percent in the inner few kpc, where the bulge is dominant. There is a gradual increase when moving to larger galactocentric radii. However, there is a large variation in for a fixed radius, especially at larger radii. Each radius contains dust cells from all around the galaxy, and from different vertical coordinates. The rings in M 31 are also not concentric. The variance per radius therefore points towards quite different environments for these dust cells, even though they reside at the same galactocentric distance.

Most interesting are the prominent spikes in the radial distribution. They point towards radii with an enhanced contribution of unevolved stellar populations to the absorbed energy (see also the circles in the bottom right panel of Fig. 6). Two moderate spikes are visible at roughly 6 and 8 kpc. They can be associated with the filamentary structures in the inner disk of M 31. A broad peak occurs between 10 and 15 kpc which corresponds to the dust cells in the main star-forming ring. Just outside the 15 kpc radius, a second sharp peak occurs, which is linked to the second large ring of Andromeda. Beyond this radius, there is an overall increase in , but without any distinct features. At these radii, the overall flux from the galaxy becomes very faint, and the dust grid contains fewer cells. The Monte Carlo noise starts to dominate, which makes it too speculative to interpret any structure in the outskirts.

It is remarkable that the dust heating fraction by unevolved stellar populations stays so low up to 10 kpc (with only a small number of dust cells as an exception). The radiation field is strong in the centre, where relatively little dust resides. Photons from evolved stellar populations may thus travel through the galaxy without being absorbed. Once they reach the dusty ring at 10 kpc, their chance of being absorbed increases strongly. This explains why the dust heating in the star-forming rings of M 31 still has an important contribution from the evolved stellar populations. This effect of non-local heating is a crucial element in energy balance studies and difficult to capture in 2D models of galaxies, or models without geometrical assumptions like energy balance SED fitting.

As a caveat to these results, we remind the reader of the residual maps in Fig. 4. There, the FUV flux is underestimated in the dustiest regions. The effect is not so large for NIR emission. We argued that this could be caused by non-uniform FUV attenuation across the disk (see Sect. 4.3). In the dustiest regions of our model, there is not enough energy coming from unevolved populations. As a result, we underestimate their contribution to the dust heating fraction. Unfortunately, it is difficult to estimate the magnitude of this effect at the present time.

6.2. Projected dust heating in M 31

We now look at the dust heating projected on the sky. While the intrinsic model and radiative transfer simulations are still in 3D, we use the 2D maps projected on the sky, recreating the observed line-of-sight for each pixel. This exercise is useful to connect the 3D model to 2D observations.

As in the previous section, we treat the contributions of the “unevolved” components as one single component. This includes absorbed energy from young and ionising stellar populations, plus the internally absorbed energy from the subgrid implementation of star-forming regions. We define the total dust heating fraction as the fraction of absorbed energy coming from a particular stellar population. For the unevolved populations for example, we denote this as Funev.. Naturally, Fevolved and Funev. add up to unity.

It is not straightforward to split these quantities into separate wavelengths because dust heating is non-local in wavelength. Dust emission at a particular wavelength is caused by the total radiation field and not by the stellar radiation at a single wavelength. Moreover, it depends in a non-linear way on the total absorbed energy. Therefore, while Funev. + Fevolved = 1 on a global scale, this is not true for individual wavelengths. Yet again it is useful to investigate whether warm dust (emitting in the MIR) is heated differently than cold dust (emitting in the submm).

We follow the definition of De Looze et al. (2014) for the wavelength dependent dust heating fraction, which is based on the absorbed energy per resolution element in the RT simulations. An alternative definition, based on the total radiation field, was first proposed by Popescu et al. (2000) and described in detail in Natale et al. (2015). It is possible that this equivalent method is quantitatively offset with our definition. However, both methods are viable to qualitatively explore the wavelength dependency of the dust heating fraction.

De Looze et al. (2014) proposed a way to approximate Fλ,evolved and Fλ,unev., while maintaining energy conservation within a single waveband. This involves running the full radiative transfer simulation with different combinations of stellar populations (see again Fig. 1, bottom panel), and recording the FIR/submm flux at every wavelength, Sλ. De Looze et al. (2014) then define: (3)and (4)where Sλ,tot is the dust emission heated by evolved, young and ionising stellar populations. Sλ,evolved is the dust emission when only the evolved stellar components are included in the simulations, and Sλ,unev. when only the young and ionising components are included. Note that there is no need to consider any missing subgrid contribution from the star-forming clouds here. Since we are dealing with emission, this is already included in the total SED of the ionising stellar component (see Sect. A.3). In the above definitions, and are the approximate heating fractions by unevolved and evolved stellar populations. They now add up to unity by definition.

This should be compared with the case where the dust heating fractions are not adjusted, that is (5)and (6)We make this comparison to estimate the uncertainty on the dust heating fraction, when looking on a per wavelength basis. Both quantities are shown for the unevolved stellar populations in Fig. 8. The amount of absorbed energy defines the shape and strength of the dust emission spectrum. The area between the curves indicates that the effect of this non-local, non-linear dust heating is most prominent in the submm regime. It can cause a discrepancy of over 5% in the dust heating fraction. Interestingly, the effect is negligible at shorter wavelengths (<50 μm).

thumbnail Fig. 8

Dust heating fraction as a function of wavelength for the unevolved stellar populations. The solid black line corresponds to the definition in Eq. (3), the dashed line to Eq. (5). The shaded area between the curves indicates the magnitude of the effect of non-local heating in wavelength space.

Open with DEXTER

The general shape of the curves is qualitatively the same, however. At almost all wavelengths the unevolved stellar populations are not the main heating source in M 31. From 1050 μm, is mostly constant at 40−50%, with several peaks. These peaks correspond to stochastially heated small grains, sensitive to the high-energy photons coming from the unevolved populations. Beyond 50 μm, the contribution of the unevolved populations drops significantly, with a minimal contribution of 15% at around 100 μm. Interestingly, rises again when it comes to heating the colder dust at submm wavelengths. This does not conform to the classical view of dust heating, where unevolved stellar populations have their peak contribution at short wavelengths (10–100 μm), but fall off after that, leaving the evolved stellar populations to heat the diffuse dust that emits at submm wavelengths (see e.g. De Looze et al. 2014; Natale et al. 2015; Bendo et al. 2015).

It seems counter-intuitive that the contribution of the unevolved stellar populations starts to increase again beyond 100 μm. A possible explanation is that the centre of the galaxy is heated by evolved stars and is warmer (or has a more intense ISRF) than the outer parts of the galaxy. In HELGA II and Draine et al. (2014) the radiation field in the centre of M 31 was found strong enough to heat the diffuse dust there to temperatures above 30 K, so dust emission peaks at shorter wavelengths. Consequently, the emission from the centre of the galaxy contributes more to the global 100 μm emission than to the global 250 μm emission. Hence, when integrating all of the emission over the galaxy, the 100 μm emission appears to have a lower in Fig. 8 than dust at 250 μm.

When we look at spatial heating maps at different wavelengths, shown in Fig. 9, we find the classical morphological segregation between disk and bulge, with the influence of the bulge almost stretching out to the 10 kpc ring. In the rings the unevolved stellar populations are dominant at short wavelengths (up to 90% at 24 μm). Their contribution slowly declines at longer wavelengths. In the large bulge of M 31, the contribution of evolved stellar populations is always high, but peaks at 100 μm. There are large spatial variations in the heating fraction at shorter wavelengths. In the submm, these variations are almost gone and the heating fraction is more uniform across the disk.

thumbnail Fig. 9

Maps of the dust heating fraction by unevolved stellar populations for a selection of wavebands. Top row (left to right): 12, 24 and 70 μm. Bottom row (left to right): 100, 250 and 500 μm.

Open with DEXTER

Going back to the global, aggregate heating fractions (not the approximate, wavelength-dependent quantities), we find that Funev. = 21%. In other words, almost 80% of the total dust emission in Andromeda is due to heating by the evolved stellar populations. This is somewhat lower than the 91% found in the previous section. This underlines an important difference between the 3D and 2D determination of the dust heating sources. In the projected view, it is possible that the full magnitude of the bulge radiation field is not captured adequately. This results in a lower contribution of the evolved stellar populations to the dust heating. Nevertheless, the evolved stellar populations dominate the global dust heating in both scenarios. Our findings are broadly in line with the dust heating analysis based on Planck colours (Planck Collaboration Int. XXV 2015). These authors report that M 31 cold dust is mainly heated by the evolved stellar populations as traced by the IRAC 3.6 μm emission. They do find that their 70/100 μm colour correlates better with the unevolved stellar populations (traced by the 24 μm emission). While we do see a rise in at wavelengths <100 μm, we do not find the unevolved stellar populations to become the dominant energy source. We come back to the point of FIR colours as dust heating tracers in the next section.

6.3. Tracers of dust heating

The analysis we have made so far, involving a full, self-consistent radiative transfer model, is a suitable way to explore the details of dust heating. On the other side, it is very demanding – both in terms of time and data3 – to construct a full RT model for each galaxy one would like to study. Hence, we have looked for possible correlations, in M 31, between the total heating fraction and some other observed properties. Here we use the global heating fraction per pixel (i.e. derived by integrating the dust SED) Funev.. For this analysis, only pixels within a galactocentric radius of 18 kpc were used. This is to avoid low S/N areas to obscure the main trends. For each correlation, we compute Kendall’s rank coefficient τ (Kendall 1938; Kendall & Gibbons 1990). This is a non-parametric correlation coefficient where τ = ± 1 denotes a perfect (anti-)correlation, and τ = 0 means there is no correlation.

In Fig. 10, we plot Funev. against NUVr colour, commonly used in the context of galaxy evolution as an observable tracer of the specific SFR (see e.g. da Cunha et al. 2010; Cortese et al. 2012; HELGA IV; Viaene et al. 2016). There is a clear positive trend (τ = 0.73), where “blue” pixels correspond to higher contributions of unevolved stellar populations to the dust heating than “red” pixels.

These results are in line with the recent insights by Boquien et al. (2016), based on panchromatic energy balance SED fitting of sub-galactic regions. They find that sSFR directly reflects the relative contributions to the total radiation field. Vice versa, they find that FIR bands like 70 μm and 100 μm are more sensitive to the absolute radiation field. This is because, FIR bands are much more influenced by dust temperature (so total energy), and less by the relative contributions to the total energy.

We therefore do not observe convincing trends between Funev. and two commonly used FIR colours to trace dust heating: the 160/250 and 250/350 flux ratios (middle and right panel of Fig. 10). The dust heating can take any fraction between 0 and 50% for most values of 160/250 or 250/350 colour. In the former, there is a significant scatter towards red 160/250 colours (or low dust temperatures). These pixels are located in the outer regions of the disk. They have cold dust, but the influence of unevolved stellar populations is relatively higher than at shorter galactocentric radii. The 250/350 colour range is particularly narrow, which makes it difficult to use as a predictor.

thumbnail Fig. 10

Density plots of the dust heating fraction from unevolved stellar populations Funev. vs. NUVr and two commonly used FIR colours. Only pixels within a galactocentric radius of 18 kpc are shown. Kendall’s correlation coefficient τ is given in each panel. A red colour indicates a small number of data points in that area of the plot, yellow at a large number.

Open with DEXTER

thumbnail Fig. 11

Density plots of the projected dust absorption luminosity from unevolved stellar populations Lunev. vs. NUVr and two commonly used FIR colours. Only pixels within a galactocentric radius of 18 kpc are shown. Kendall’s correlation coefficient τ is given in each panel. A red colour indicates a small number of data points in that area of the plot, yellow at a large number.

Open with DEXTER

Bendo et al. (2015) also found that the fraction of dust heating is hard to predict using FIR colours. On the other hand, they were able to link FIR colours to the absolute energy output from evolved and unevolved stellar populations, traced by 3.6 μm and Hα emission, respectively. The same results were confirmed for M 31 by Planck Collaboration Int. XXV (2015). In our model framework, we can push this one step further and look at the total absorbed energy per pixel for each stellar population. In Fig. 11, we plot this value for the unevolved stellar populations (Lunev.). This quantity is similar to the Labs parameters used in Sect. 6.1, but now projected to the 2D pixels on the sky.

However, the trends we observe for Lunev. are less convincing then for Funev.. NUVr is still the strongest correlation (τ = 0.21), but still weak and clearly bimodal. For the 250/350 colour, the trend is weak (τ = 0.19) and the colour range is narrow. Again, the 160/250 colour (τ = 0.07) has a large variance and is also not suitable as a predictor for Lunev.. This is rather surprising because in M 31, the 160/250 flux ratio should be most sensitive to the shape of the cold dust emission spectrum (see for example Fig. 1). The peak of the dust distribution lies at 160 μm, while 250 μm is part of the Rayleigh-Jeans tail of the cold dust emission. One would naively expect a link to the total absorbed energy. It appears that the variation in dust temperature on a local, pixel-by-pixel scale is too large to capture in a single FIR colour.

We now look for possible correlations between the dust heating and physical properties. Our previous investigation of M 31 (HELGA IV) offers a range of physical parameters for each pixel. In Fig. 12, the most commonly used physical parameters are plotted against Funev.. There is a highly scattered positive trend with Mdust (τ = 0.36). This could be interpreted in the context of star formation, where dusty regions are often associated with star-forming regions. Consequently, these regions are dominated by the radiation field of unevolved stellar populations. This is in accordance with a positive trend with SFR (τ = 0.26). However, the scatter in this correlation is large, suggesting other elements are in play as well. Interestingly, no clear trend with the total dust luminosity, Ldust, was found either (τ = −0.12). Infrared luminous galaxies often point at high levels of star formation, and so one would naively expect a higher contribution of unevolved stellar populations to the heating fraction. However, for the individual pixel-regions in M 31, this is clearly not the case.

We find a negative correlation (τ = −0.34) with stellar mass (M). If one assumes that the bulk of stellar mass is in evolved stars, then it is not surprising that the pixels with high stellar mass correspond to low Funev.. But here as well, the trend is ambiguous and suggests the influence of a third parameter. We further investigate this by looking at the specific dust mass (Mdust/M). This quantity gives an idea of the dust content, relative to the stellar mass. High Mdust/M values point at dusty regions and those can indeed be associated with high Funev. (τ = 0.64). The heating fraction of unevolved stellar populations approaches zero for the lowest specific dust masses. This suggests that at low dust content, there is also no significant radiation field from young or ionising stellar populations.

Another significant correlation is with sSFR (τ = 0.51). This correlation is likely to be intrinsic as the trends with SFR or M separately show much more scatter. This trend is not completely unexpected, since we already found a positive link between Funev. and NUVr colour, and this colour correlates with sSFR in M 31 (HELGA IV). sSFR traces the hardness of the UV radiation field (Ciesla et al. 2014) and is linked to several FIR/submm colours in galaxies (e.g. Boselli et al. 2010, 2012; Boquien et al. 2011; Cortese et al. 2014). Moreover, De Looze et al. (2014) found a tight link between the dust heating fraction of unevolved stellar populations and sSFR for M 51. We overplot their data points and best fit curve in the central, bottom-row panel of Fig. 12. In M 51, higher fractions of Funev. can be achieved compared to Andromeda. The regions of M 31 follow the average relation for the M 51 regions closely for all but the lower sSFR values. In fact, both datasets seem to lie in the continuation of one another. This is quite intriguing, since the physical size along the major axis of the M 31 pixels is about 3.5 times smaller than for the M 51 pixels. On top of that is the difference in inclination of both galaxies, where a far greater distance is probed along the line of sight in M 31.

thumbnail Fig. 12

Density plots of the dust heating fraction from unevolved stellar populations Funev. vs. a number of physical parameters, obtained from HELGA IV. A red colour indicates a small number of data points in that area of the plot, yellow at a large number. Kendall’s correlation coefficient τ is given in each panel. First row: dust mass, stellar mass, dust-to-stellar mass ratio. Second row: star formation rate, specific SFR, temperature of the cold dust in the ISM. The black line is a best fit linear model, the green line is the relation found by De Looze et al. (2014) for M 51. Parameters are derived according to the MAGPHYS SED model.

Open with DEXTER

The bottom line is that extensive galaxy parameters such as Mdust, Ldust, Mstar or SFR are not suitable to disentangle the dominant dust heating sources. A better option is to resort to intensive parameters that somehow express the ratio of evolved and unevolved stellar populations, for example Mdust/M or sSFR. While there is still some variance in the correlation with sSFR, it is small enough to identify the dominant heating source. Of course, a more thorough investigation involving multiple galaxies is desirable to quantify this relation.

7. Conclusions

We have constructed a highly detailed model for the Andromeda galaxy to investigate the dust heating mechanisms in this galaxy. Our model is based on observed morphologies and uses 3D panchromatic radiative transfer to simulate the dust-starlight interactions in a realistic setting. The main points of this work are:

  • The integrated SED of M 31 is fitted well and theresulting attenuation curve is consistent with observations, butwith a broader UV bump. We are able to constrain 2 of our 3 freeparameters: the dust mass and the luminosity of the young stellarcomponent. The intrinsic luminosity of the ionising stellarpopulation is more difficult to constrain.

  • The model is able to reproduce the observed morphologies fairly well from FUV to submm wavelengths. The median (absolute value) deviation between model and observations across all bands is 22%. The flux in the rings is generally underestimated, and the flux in the inter-ring regions is overestimated. Lowering the vertical scale height of the galaxy somewhat mitigates this effect, but cannot resolve it. The discrepancies are a combination of deprojection effects, variations in the nature and size distribution of the dust grains, and the subgrid treatment of the star-forming regions. With an inclination of 77.5°, M 31 is about the limiting case for these kind of deprojected radiative transfer models.

  • The dust in Andromeda is mainly heated by the evolved stellar populations. From a 3D analysis of the radiation field, we find that 91% of absorbed stellar radiation originates in evolved stellar populations. This high value is mainly due to the bright bulge, which dominates the radiation field out to the main star-forming ring at 10 kpc. Inside and beyond the star-forming ring, the contribution of unevolved stellar populations (i.e. ionising and non-ionising young stellar populations) to the radiation field increases, but usually remains in the 10−30% range.

  • The sSFR (and its observational counterpart, NUVr colour) is a promising tracer for the total dust heating fraction and thus for the relative contributions to the total IR emission. We find that regions in M 31 match the best fit relation derived from M 51 pixels. In fact, the two datasets make a rather smooth and continuous sequence. More research is required to assess whether sSFR (and NUVr) is a general tracer of dust heating fractions in galaxies at a local scale.

Our study has shown that the heating of dust by stellar populations is a complex problem with a large influence of geometry in three dimensions. Effects like non-local heating make it difficult to draw conclusions from 2D “on-sky” analysis. The contribution of evolved stellar populations is dominant in M 31, and can also be significant in other galaxies. One should therefore be careful to directly link dust emission to the properties of unevolved stellar populations.

As a final remark, we want to underline that this is one of the first attempts to construct a detailed geometrical model of stars and dust of a resolved and moderately inclined galaxy, based purely on 3D radiative transfer simulations. Because 3D radiative transfer models are computationally demanding, we were forced to make a number of simplifying assumptions on the geometry, and limit the number of free parameters in the model. In the future, we will refine our method on several points and apply it to galaxies of different morphologies and inclination.


2

SKIRT is publicly available at http://skirt.ugent.be

3

For 512 wavelengths and 2 × 107 photon packages per wavelength, running one model takes about 4000 CPU hours and 15 GB of RAM.

4

In this study, we consider stellar populations to be evolved when their UV radiation is not significant any more and they mainly produce optical/NIR emission.

Acknowledgments

This project has received funding from DustPedia, a European Union’s Seventh Framework Programme for research, technological development and demonstration under grant agreement No. 606874. S.V. gratefully acknowledges the support of the Flemish Fund for Scientific Research (FWO-Vlaanderen). I.D.L. is a post-doctoral researcher of the FWO-Vlaanderen (Belgium). E.T. and A.T. acknowledge the support from the ESF grants IUT40-2, IUT26-2. M. Boquien acknowledges funding by the FIC-R Fund, allocated to the project 30321072. This work has been realised in the frame of the CHARM framework (Contemporary physical challenges in Heliospheric and AstRophysical Models), a phase VII Interuniversity Attraction Pole (IAP) programme organised by BELSPO, the BELgian federal Science Policy Office. We thank all the people involved in the construction and the launch of Herschel. SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC and UKSA (UK); and NASA (USA). HIPE is a joint development (are joint developments) by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia.

References

Appendix A: Model components

thumbnail Fig. A.1

Effect of galaxy thickness on residuals in the FUV band (top row) and at 250 μm (bottom row). Thefirst column shows residuals for the original thickness, middle and right columns for 2/3 and 1/2 of that thickness, respectively.

Open with DEXTER

For our aim to investigate dust heating properties, we simplify the galaxy to a handful of “average” galaxy components. They are listed in Table 1 and we discuss them in more detail below.

Appendix A.1: Evolved stellar populations

The 3D structure of the evolved stellar populations (>100 Myr 4) is a necessary input for our model. We simplify this distribution to two components: a disk and a bulge. For M 31, we derive the distribution of the evolved stellar populations from the IRAC 3.6 μm image. The derived geometries must be deprojected to construct their 3D equivalent. Deprojecting an image assumes all features reside in the same (flat) disk. The bulge, however, has a large vertical dimension, so it would be smeared out into a bar-like structure in the deprojected image. This issue is avoided by subtracting the bulge from the IRAC 3.6 μm image before deprojection. An accurate representation of the bulge of M 31 is needed in 2D to subtract from the 3.6 μm image but also in 3D, to implement it in our RT model. We thus require an analytical model that can be fitted to the 2D image, but can be deprojected in a realistic way. The big issue here is the boxiness of Andromeda’s bulge (Beaton et al. 2007).

Several analytical bulge-disk decompositions are available for M 31 (e.g., Courteau et al. 2011; Tamm et al. 2012), but they all produce elliptical isophotes. Using N-body simulations, Athanassoula & Beaton (2006) managed to construct a 3D model that approximately reproduced the boxy bulge of M 31. They argue that a bar is responsible for the observed boxiness. Unfortunately, they do not provide an analytical prescription. The fitting of 2D analytical boxy bulges is possible using for example GALFIT (Peng et al. 2010). However, to our knowledge there is no 3D analytical prescription to mimic a boxy bulge. Deriving such a prescription falls beyond the scope of this paper. We therefore resort to the next best 3D representation of Andromeda’s bulge: the flattened Einasto profile derived by Tamm et al. (2012; see their Eq. (B.4)). It describes the surface brightness S(r) as a function of radius r: (A.1)This prescription has two free parameters: the concentration radius of ac = 1.155 kpc which contains half of the total brightness, and a structure parameter N = 2.7. The numerical constant dN is only dependent on N and ensures that half of the brightness is containted in ac. An intrinsic flattening in the vertical direction of q = 0.72 was imposed in the sense that (A.2)We find a difference in position angle of 13.8° between the major axis of the bulge and the major axis of the disk. The flattened Sérsic (1963) bulge derived by Courteau et al. (2011) was also tested, but the resulting disk (subtracting the bulge from the IRAC 3.6 μm image) showed stronger artefacts.

Tamm et al. (2012) found a bulge-to-disk ratio (B/D) of 0.54 from SDSS-i band imaging and we find that this ratio also works well at 3.6 μm. Upon subtraction of the bulge model from the 3.6 μm image, some pixels show a slight over- or under-subtraction caused by the small deviation of the bulge in M 31 from the smooth Einasto profile. To ensure no negative emission in the centre of the disk, we replace the inner circle with a radius of 5 pixels (685 pc) by the local disk background just like we have masked M 32 (see Sect. 2). This correction is unlikely to affect the radiation field in the nucleus of M 31 as the bulge is dominant here.

The corrected, bulge-subtracted 3.6 μm image was used as input geometry for the evolved stellar disk. Following De Looze et al. (2014), we assume an exponential vertical profile for the disk. In practice, this introduces an extra parameter: the vertical scale height. We rely on the modelling efforts of Tempel et al. (2010), who derive radial scale lengths and flattening parameters for Einasto profile fits to M 31. They take the 3D nature of the galaxy into account by concurrently modelling the extinction by dust. Their results are equivalent to a vertical scale height of 538 pc for the evolved stellar disk. This value also matches the average scale height found by De Geyter et al. (2014) for a sample of 12 edge-on galaxies.

A panchromatic approach requires each model component to have a luminosity at every UV to submm wavelength. In practice this is achieved by selecting a suitable template SED for each component. The luminosity is then scaled to match a certain value at a representative wavelength. For the evolved stars we make use of the Bruzual & Charlot (2003) simple stellar populations (SSPs). In this paper, we make use of their low resolution version of the Padova 1994 (Bressan et al. 2001) model, which uses the Chabrier (2003) initial mass function. The templates are characterised by metallicity and age.

For the bulge, we assume an average age of 12 Gyr and Solar metallicity (Z = 0.02). In the disk, we assume the same metallicity and an average age of 8 Gyr for the stars. The estimates for age and metallicity are taken from Saglia et al. (2010), who analysed several long-slit spectra across the bulge and inner disk of M 31. We normalise this SED for both bulge and disk to the total IRAC 3.6 μm luminosity and use B/D = 0.54. This is the only stellar SED we can unambiguously normalise as the IRAC 3.6 μm band is a “pure” band, that is, it traces only one component (evolved stars) and suffers least from attenuation by dust.

Appendix A.2: Young stellar populations

Following the general picture of large galaxies, we assume that the young stellar populations reside in a disk that is thinner than the evolved stellar disk. These stars are bright in the UV, but can not ionise hydrogen molecules. They should have had the time to migrate from or dissolve their dusty birth clouds, which is typically after 106−107 yr (Hartmann et al. 2001; Murray 2011; Bailey & Basu 2014). Their emission is therefore only obscured by dust in the diffuse ISM.

Evolved stars can contribute significantly to the FUV flux (e.g. Kennicutt et al. 2009). The bulge of M 31 is known to have old stars which contribute to the UV flux (Saglia et al. 2010; Rosenfield et al. 2012). We correct for this using the formula derived by Leroy et al. (2008) to obtain , the surface brightness from young stellar populations in the FUV: (A.3)In HELGA III this was applied to M 31 and αFUV = 8.0 × 10-4 was found, based on the same FUV and IRAC 3.6 μm observations of M 31 as used in this work. Note that this correction only affects the geometrical distribution of the young stellar population. In practice, this distributes most of these stars in the disk, and very little in the bulge. The corrected FUV image is used as the 2D input geometry for the young stellar populations. Similar to the disk of evolved stellar populations in the previous section, we add an exponential profile in the third dimension. We set the vertical scale height to be 190 pc, following Tempel et al. (2010).

We assume that the relative surface brightness in the FUV band gives a good approximation of the distribution of these stars, even though the intrinsic flux may differ from the observed one. This is a first order approximation as dust extinction is not uniform across the disk. Regions with visible dust lanes are likely to absorb a larger fraction of the FUV radiation from young stars. This will decrease their relative importance in the map. Vice versa, the distribution of young stars will be boosted in the inter-ring regions.

Ideally, one would independently vary the amplitude of the FUV emissivity (for the fixed assumed vertical distribution) for each resolution element on the images. It is clear that such an approach is currently intractable, as it would require far too many RT calculations to explore all permutations of the FUV brightness distribution.

As SED template for this component, a Bruzual & Charlot (2003) template was again used. There is a known metallicity gradient in M 31 (e.g., Mattsson et al. 2014, HELGA V), but it is shallow within the optical disk. As we are working with average stellar components, we keep the metallicity at the Solar value. We set a mean age of 100 Myr, following De Looze et al. (2014). Using a different age would mainly affect the UV flux of this component, and change the normalization of the SED. However, we found that this has only little effect on the FIR/submm SED. The SED is normalised to an intrinsic FUV luminosity.

The total intrinsic FUV in M 31 can of course not be measured directly, so we leave this parameter free. We can, however, make an educated guess. Using Eq. (A.3)we find that the contribution of young stellar populations to the total FUV flux is about 84%. The extinction in the FUV band of the global SED model of HELGA IV is 45% in flux. When we correct the observed FUV flux with these two factors, we obtain an initial guess of the intrinsic FUV flux for the young stellar component.

Appendix A.3: Ionising stellar populations

Embedded star formation is a crucial element in the energy balance between UV/optical (absorption) and FIR/submm (emission). This was first recognised by Silva et al. (1998) and Popescu et al. (2000), and confirmed by consequent work (Misiriotis et al. 2001; Baes et al. 2010; Popescu et al. 2011; MacLachlan et al. 2011; De Looze et al. 2012b,a, 2014). The ionising stars (10 Myr) are still embedded in their birth clouds and produce a hard, ionising radiation. Due to the dust in these clouds it is near impossible to trace their emission directly. Ionising stars ionise their surrounding gas. The ionised H ii gas can be traced through its Hα emission. Devereux et al. (1994) mapped a large part of Andromeda in Hα. Unfortunately their map does not cover the entire disk and there are several artefacts present after continuum subtraction. Another indirect tracer is the warm dust surrounding the stellar birth clouds. The emission of warm dust was already derived for M 31 in HELGA III, using the MIPS 24 μm image. This waveband can hold a significant contribution of warm dust residing in the atmospheres of evolved stars, Leroy et al. (2008) therefore derive the following correction to obtain , the surface brightness of the ionising stellar component in the FUV: (A.4)This was also applied to M 31 in HELGA III. They found that α24 = 0.1. Despite removing this contribution, the map may still be contaminated by emission from diffuse dust, heated by the general ISRF (see also Kennicutt et al. 2009). One should keep in mind that this seriously hampers the use of the 24 μm band as quantitative tracer of the total SFR. However, as we are only interested in the spatial distribution of the obscured SFR, the effect of diffuse dust emission is less important. In the third (vertical) dimension, we again use an exponential profile, as for the young and evolved stellar disks. The vertical scale height is set to 190 pc, similar as for the young stellar disk.

For the SED of the ionising stars we follow the approach of De Looze et al. (2014) and use the standard MAPPINGS III template SEDs for obscured star formation from Groves et al. (2008). These SEDs are isotropic representations of the H ii regions around massive clusters of young stars. We adopt the same values for metallicity (Z = 0.02), compactness (), pressure of the surrounding ISM (P0 = 1012 K/m3) and cloud covering factor (f = 0.2) as in De Looze et al. (2014) and refer the reader to this work for more specifics about the use of these templates in SKIRT.

As for the young stellar component, we also normalise this SED to the total intrinsic FUV emission of the ionising stars. For the same reason (extinction by dust) we leave this intrinsic luminosity as a free parameter. We can again make an initial guess for this parameter: the Groves et al. (2008) templates are normalised to a SFR of 1 M yr-1. We take the template luminosity in the FUV band and scale this to 0.2 M yr-1, which is a good approximation of the SFR in M 31 (Tabatabaei & Berkhuijsen 2010; HELGA III; HELGA IV; Rahmani et al. 2016). We use this value as an initial guess for the intrinsic FUV luminosity of the ionising stars.

Appendix A.4: Interstellar dust

To map the dust in a galaxy one can resort to one of the FIR/submm bands. However, as dust emission peaks at different wavelengths depending on its temperature, choosing a single band is not straightforward.

Alternatively a map of the UV attenuation AFUV can be created (Montalto et al. 2009; De Looze et al. 2014). This method takes into account multiple FIR/submm wavebands and links this to the FUV morphology, where dust can be seen in extinction. Consequently, it should produce a tighter constraint on the dust morphology. There are many recipes available for AFUV, based on different calibrations (Cortese et al. 2008; Hao et al. 2011; Boquien et al. 2012; Viaene et al. 2016). The drawback of these recipes is that they indirectly rely on simplified geometries to convert optical depth to dust attenuation.

For M 31 in particular, a high-resolution dust attenuation map was constructed by Dalcanton et al. (2015) from extinction of individual stars. Unfortunately, their map only covers one third of our field of view. Draine et al. (2014) constructed a dust mass map based on pixel-by-pixel SED fits from MIR to submm wavelengths. Independently, in HELGA IV, we performed a panchromatic pixel-by-pixel fit of M 31. Their dust map is also consistent within their model with the attenuation of UV and optical light. They found that their map was consistent with that of Draine et al. (2014), within the model uncertainties.

We will use the HELGA IV dust mass map as an input for the relative dust mass distribution in our model, but do not fix the total dust mass. To this 2D geometry we again add an exponential profile in the third dimension. Tempel et al. (2010) find a dust scale height of 238 pc, which yields a dust-to-stellar scale height 0.5. This is consistent with RT models of edge-on galaxies (Xilouris et al. 1999; Bianchi 2007; De Geyter et al. 2014).

We use the THEMIS dust model compiled by Jones et al. (2013). The model consists of amorphous hydrocarbons and silicates and we use their parametrization for diffuse dust. THEMIS is a framework of dust properties, based of a vast set of laboratory data, and characterised by a series of free parameters. We adopt the parametrization that best suits Milky Way dust (Jones et al. 2013; Ysard et al. 2015). This model naturally produces the spectral shape in the MIR/FIR/submm, and is consistent in explaining both extinction and emission by the same dust mixture for diffuse dust in the Milky Way (Ysard et al. 2015).

The dust component in our RT model is normalised by the total dust mass. The amount of diffuse dust was determined from the dust mass map from HELGA IV. However, a correction for dust in local star-forming clouds is necessary. Indeed, the SED of the obscured ionising stars holds a certain amount of dust (see Sect. A.3). The total observed dust mass in M 31 is the sum of this star-forming cloud dust, and the diffuse dust component. This introduces no extra free parameters, but couples the normalization of the total dust mass to the normalization of the SED of the ionising stars.

We leave the total dust mass as a free parameter. As an initial guess for this value we take the total dust mass given from the MAGPHYS dust map, which is normalised using the Dunne et al. (2000) dust extinction coefficient κ350 μm = 0.454 m2/kg. In the Jones et al. (2013) model, the average κ350 μm becomes 0.305 m2/kg. This enhances the total dust mass by a factor of 1.5. However, we note that this average κ350 μm should in general not be used to convert dust emission into dust masses. The value for the extinction coefficient varies strongly with dust mixture and is still an average factor, with large variance. Still, for our initial guess, it is sufficient.

Appendix B: Observed fluxes

Table B.1

Integrated fluxes for M 31 adopted in this paper, listed by increasing central wavelength.

Table B.1 lists the integrated fluxes for M 31 used in this paper. These are the integrated fluxes from the region we consider in our model. This an ellipse with centre at α = 00:43:06.28 and δ = +41:21:12.22, semi-major axis of 1.425°, semi-minor axis of 0.400° and a position angle of 38.1°. The values may differ slightly from other literature measurements, based on different apertures.

Appendix C: Alternative scale heights

The best-fitting model for M 31 tends to underestimate the flux in the rings, and overestimate the flux in the inter-ring regions. We test whether this could be induced by the addition of a vertical dimension.

We therefore ran the best-fitting model with smaller scale heights for stars and dust. We show the effect on the FUV and 250 μm band in Fig. A.1. If the galaxy model is 1/3 thinner, we observe little effect on the residuals. The median deviation from observations decreases from 34% to 30% in the FUV and from 23% to 22% at 250 μm. Running the model at only 1/2 of the original thickness slightly decreases these values again. The median deviation from observation now becomes 29% in the FUV and 21% at 250 μm.

So, even with a lower scale length, the model still overestimates the flux in the inter-ring regions and underestimates the flux in the rings. We therefore choose not to change our initial assumptions about the scale height of the individual components as these are compatible with the empirical optical depth study of Tamm et al. (2012).

All Tables

Table 1

Overview of the different stellar populations (SP) and dust component in our model.

Table 2

Overview of the parameters in our model.

Table 3

Star formation rates for M 31, derived from synthetic broad band fluxes from our model SED.

Table B.1

Integrated fluxes for M 31 adopted in this paper, listed by increasing central wavelength.

All Figures

thumbnail Fig. 1

Panchromatic SED of Andromeda. The black line is the best fitting RT model, run at high resolution (512 wavelengths). The green points are the observed integrated luminosities for M 31. The red, blue and cyan lines represent the SEDs for simulations with only one stellar component: evolved, young and ionising, respectively. The interstellar dust component is still present in these simulations. The bottom panel shows the residuals between the observations and the best fitting (black) model.

Open with DEXTER
In the text
thumbnail Fig. 2

PDFs of the three free parameters in our model optimization: the total dust mass Mdust (left), the intrinsic FUV luminosity from the young stellar component (middle) and the intrinsic FUV luminosity from the ionising stellar component (right). Black dashed lines are the parameter values for the best fitting model. Red dashed lines are the median (50th percentile) values from the PDFs.

Open with DEXTER
In the text
thumbnail Fig. 3

Attenuation laws from our model, normalised to the V-band attenuation. The black solid line is the global law. The dashed line reflects the attenuation law for the diffuse dust and the dotted black line for the dust in star-forming regions. Several literature measurements are shown; red curve: MW curve for R(V) = 3 from Fitzpatrick & Massa (2007), blue dashed line: SMC bar region (Gordon et al. 2003), and green dash-dotted line: attenuation curve from Battisti et al. (2016).

Open with DEXTER
In the text
thumbnail Fig. 4

Comparison of the model with observations in selected wavebands. First column show the observed images, second column the model images. The 3rd column shows the residuals derived as in Eq. (1). From top to bottom row: FUV, g, 3.6 μm, 8 μm, 24 μm, 100 μm and 250 μm.

Open with DEXTER
In the text
thumbnail Fig. 5

Observed (left) and model (right) colours tracing the FIR dust emission peak. A colour here corresponds to the flux ratio: log (flux1/ flux2), that is, blue colours indicate a high flux ratio, red colours a low ratio.

Open with DEXTER
In the text
thumbnail Fig. 6

Face-on model images of M 31 in the FUV, 3.6, 8.0, 24, 100 and 250 μm. Colour scale is dif ferent for each image for clarity. The green circles on the 250 μm image (bottom right) correspond to radii of 10 kpc and 15 kpc from the centre.

Open with DEXTER
In the text
thumbnail Fig. 7

Left: 2D histogram comparing the integrated absorbed luminosity per dust cell including all stellar components with the same quantity, but as the sum of separate simulations including only evolved, young and ionising stellar populations. Middle: distribution of the bolometric absorption fraction of unevolved stellar emission per dust cell: . The histogram is weighted by mass fraction of each dust cell and normalised. The red line indicates the median . Right: 2D histogram showing the radial distribution of . The bins are weighted by the mass fraction of the dust cells. Red indicates a high number of data points, blue a low number. The spread at a fixed radius is a combination of dust cells at the same radius and with different vertical location in the disk. The white line is the weighted mean heating fraction per radius.

Open with DEXTER
In the text
thumbnail Fig. 8

Dust heating fraction as a function of wavelength for the unevolved stellar populations. The solid black line corresponds to the definition in Eq. (3), the dashed line to Eq. (5). The shaded area between the curves indicates the magnitude of the effect of non-local heating in wavelength space.

Open with DEXTER
In the text
thumbnail Fig. 9

Maps of the dust heating fraction by unevolved stellar populations for a selection of wavebands. Top row (left to right): 12, 24 and 70 μm. Bottom row (left to right): 100, 250 and 500 μm.

Open with DEXTER
In the text
thumbnail Fig. 10

Density plots of the dust heating fraction from unevolved stellar populations Funev. vs. NUVr and two commonly used FIR colours. Only pixels within a galactocentric radius of 18 kpc are shown. Kendall’s correlation coefficient τ is given in each panel. A red colour indicates a small number of data points in that area of the plot, yellow at a large number.

Open with DEXTER
In the text
thumbnail Fig. 11

Density plots of the projected dust absorption luminosity from unevolved stellar populations Lunev. vs. NUVr and two commonly used FIR colours. Only pixels within a galactocentric radius of 18 kpc are shown. Kendall’s correlation coefficient τ is given in each panel. A red colour indicates a small number of data points in that area of the plot, yellow at a large number.

Open with DEXTER
In the text
thumbnail Fig. 12

Density plots of the dust heating fraction from unevolved stellar populations Funev. vs. a number of physical parameters, obtained from HELGA IV. A red colour indicates a small number of data points in that area of the plot, yellow at a large number. Kendall’s correlation coefficient τ is given in each panel. First row: dust mass, stellar mass, dust-to-stellar mass ratio. Second row: star formation rate, specific SFR, temperature of the cold dust in the ISM. The black line is a best fit linear model, the green line is the relation found by De Looze et al. (2014) for M 51. Parameters are derived according to the MAGPHYS SED model.

Open with DEXTER
In the text
thumbnail Fig. A.1

Effect of galaxy thickness on residuals in the FUV band (top row) and at 250 μm (bottom row). Thefirst column shows residuals for the original thickness, middle and right columns for 2/3 and 1/2 of that thickness, respectively.

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.