# The Herschel Exploitation of Local Galaxy Andromeda (HELGA)

VII. A SKIRT radiative transfer model and insights on dust heating*

S. Viaene<fnref n="1" />, M. Baes<fnref n="1" />, A. Tamm<fnref n="2" />, E. Tempel<fnref n="2" />, G. Bendo<fnref n="3" />, J. A. D. L. Blommaert<fnref n="4" />,<fnref n="5" />, M. Boquien<fnref n="6" />, A. Boselli<fnref n="7" />,
P. Camps<fnref n="1" />, A. Cooray<fnref n="8" />, I. De Looze<fnref n="1" />,<fnref n="9" />, P. De Vis<fnref n="1" />, J. A. Fernández-Ontiveros<fnref n="10" />, J. Fritz<fnref n="11" />, M. Galametz<fnref n="12" />, G. Gentile<fnref n="1" />,
S. Madden<fnref n="13" />, M. W. L. Smith<fnref n="14" />, L. Spinoglio<fnref n="10" />, and S. Verstocken<fnref n="1" />

Received 5 July 2016 / Accepted 26 September 2016

## ABSTRACT

The radiation from stars heats dust grains in the diffuse interstellar medium and in star-forming regions in galaxies. Modelling this interaction provides information on dust in galaxies, a vital ingredient for their evolution. It is not straightforward to identify the stellar populations heating the dust, and to link attenuation to emission on a sub-galactic scale. Radiative transfer models are able to simulate this dust-starlight interaction in a realistic, three-dimensional setting. We investigate the dust heating mechanisms on a local and global galactic scale, using the Andromeda galaxy (M31) as our laboratory.

We have performed a series of panchromatic radiative transfer simulations of Andromeda with our code SKIRT. The high inclination angle of M31 complicates the 3D modelling and causes projection effects. However, the observed morphology and flux density are reproduced fairly well from UV to sub-millimeter wavelengths. Our model reveals a realistic attenuation curve, compatible with previous, observational estimates. We find that the dust in M31 is mainly (91% of the absorbed luminosity) heated by the evolved stellar populations. The bright bulge produces a strong radiation field and induces non-local heating up to the main star-forming ring at 10 kpc.

The relative contribution of unevolved stellar populations to the dust heating varies strongly with wavelength and with galactocentric distance. The dust heating fraction of unevolved stellar populations correlates strongly with *NUV*-*r* colour and specific star formation rate. These two related parameters are promising probes for the dust heating sources at a local scale.

**Key words.** galaxies: individual: M 31 – galaxies: ISM – infrared: ISM – galaxies: fundamental parameters – dust, extinction – radiative transfer

## 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

<footnote n="1">¹ Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281, 9000 Gent, Belgium
e-mail: sebastien.viaene@ugent.be</footnote>

<footnote n="2">² Tartu Observatory, Observatooriumi 1, 61602 Tõravere, Estonia</footnote>

<footnote n="3">³ UK ALMA Regional Centre Node, Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="4">⁴ KU Leuven, Instituut voor Sterrenkunde, Celestijnenlaan 200D 2401, 3001 Leuven, Belgium</footnote>

<footnote n="5">⁵ Vrije Universiteit Brussel, Department of Physics and Astrophysics, Pleinlaan 2, 1050 Brussels, Belgium</footnote>

<footnote n="6">⁶ Unidad de Astronomía, Fac. Cs. Básicas, Universidad de Antofagasta, Avda. U. de Antofagasta 02800, Antofagasta, Chile</footnote>

<footnote n="7">⁷ Laboratoire d'Astrophysique de Marseille - LAM, Université d'Aix-Marseille & CNRS, UMR 7326, 38 rue F. Joliot-Curie,
13388 Marseille Cedex 13, France</footnote>

<footnote n="8">⁸ Department of Physics and Astronomy, University of California, Irvine CA 92697, USA</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="10">¹⁰ Istituto di Astrofisica e Planetologia Spaziali (INAF-IAPS), via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>

<footnote n="11">¹¹ Instituto de Radioastronomía y Astrofísica, IRAF, UNAM, Campus Morelia, A.P. 3-72, C.P. 58089, Mexico</footnote>

<footnote n="12">¹² European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching-bei-München, Germany</footnote>

<footnote n="13">¹³ Laboratoire AIM, CEA, Université Paris Diderot, IRFU/Service d'Astrophysique, Bat. 709, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="14">¹⁴ School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote>* Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>
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 M51. 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 (M31). 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 M31 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 M31 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 M31 is much larger and has a higher inclination angle than M51 (in this work, we will use $i = 77.5^{\circ}$, which is the average value for the disk of M31 as derived from HI data by Corbelli et al. 2010). On the other hand, M31 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 M31 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 M31

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, M31 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 M31 is now publicly available online<fnref n="1" /> 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 M31 itself. Second, all images were convolved to match the SPIRE 500 µm beam, which is the largest in our set. Finally, the frames were re-binned 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 M31 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 $D_{25}$ 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, M32 and NGC 205 are not. Both dwarf galaxies are satellites of Andromeda and M32 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). M32 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 M32 and replace the emission from M32 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 M31. 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 M31 and are discussed in Sect. 3.2.

# 3. Constructing a 3D model

## 3.1. Radiative transfer simulations

We make use of SKIRT<fnref n="2" /> (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 M31 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 M31 ($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.

To keep the balance between a realistic model and computational feasibility, the number of components has to be limited. We split the stars in M31 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 M31 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

<footnote n="1">¹ http://irsa.ipac.caltech.edu/data/Herschel/HELGA</footnote>

<footnote n="2">² SKIRT is publicly available at http://skirt.ugent.be</footnote>
<table><caption>**Table 1.** Overview of the different stellar populations (SP) and dust component in our model.</caption><thead><tr><th>Component</th><th>2D geometry</th><th>Vertical dimension</th><th>SED template</th><th>Normalization</th></tr></thead><tbody><tr><td><em>Bulge</em></td><td></td><td></td><td></td><td></td></tr><tr><td>Evolved SP (12 Gyr)</td><td>3D Einasto profile</td><td></td><td>Bruzual &amp; Charlot (2003)</td><td>3.6 µm</td></tr><tr><td><em>Disk</em></td><td></td><td></td><td></td><td></td></tr><tr><td>Evolved SP (8 Gyr)</td><td>Corrected<sup>a</sup> 3.6 µm image</td><td>Exponential profile</td><td>Bruzual &amp; Charlot (2003)</td><td>3.6 µm</td></tr><tr><td>Young SP (100 Myr)</td><td>Corrected<sup>b</sup> FUV image</td><td>Exponential profile</td><td>Bruzual &amp; Charlot (2003)</td><td>FUV</td></tr><tr><td>Ionising SP (10 Myr)</td><td>Corrected<sup>b</sup> 24 µm image</td><td>Exponential profile</td><td>MAPPINGS III<sup>c</sup></td><td>FUV</td></tr><tr><td><em>Dust</em></td><td></td><td></td><td></td><td></td></tr><tr><td>Interstellar dust</td><td>MAGPHYS dust mass map<sup>d</sup></td><td>Exponential profile</td><td>THEMIS<sup>e</sup> dust mix</td><td>Total dust mass</td></tr></tbody></table>

**Notes.** We used an inclination angle of 77.5° to deproject the input images. (a) The Einasto bulge was first subtracted from the image. (b) The contribution from the evolved stars was first subtracted from the image. (c) The SED templates for obscured star formation from Groves et al. (2008). (d) The dust mass map derived from the pixel-by-pixel MAGPHYS energy balance SED models from Viaene et al. (2014). (e) Jones et al. (2013).

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 M31, 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 M31 model. Note that these cells are not equal in size. Their volume ranges between $(30-5000 \text{ pc})^3$, with a median value of $(89 \text{ 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 \times 10^8 L_\odot$.

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 (1–30 µm) to capture the aromatic features. The best fitting model SED is determined by the lowest $\chi^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 $\chi^2$ when using a uniform weighing because there are many data points in this regime. However, to study the dust heating in M31, 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.

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 \times 10^6$ 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 M31 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
<table><caption>**Table 2.** Overview of the parameters in our model.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Initial guess</th><th>Range</th><th>Best fit</th></tr></thead><tbody><tr><td><strong>Evolved stars bulge (&gt;100 Myr)</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Total luminosity at 3.6 µm</td><td>λL<sub>λ</sub><sup>bulge</sup> [L<sub>☉</sub>]</td><td>1.53 × 10<sup>9</sup></td><td>fixed</td><td></td></tr><tr><td>Central radius of the Einasto profile</td><td>a<sub>c</sub> [kpc]</td><td>1.155</td><td>fixed</td><td></td></tr><tr><td>Einasto structure parameter</td><td>N</td><td>2.7</td><td>fixed</td><td></td></tr><tr><td>Flattening factor</td><td>q</td><td>0.72</td><td>fixed</td><td></td></tr><tr><td><strong>Evolved stars disk (&gt;100 Myr)</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Total luminosity at 3.6 µm</td><td>λL<sub>λ</sub><sup>evolved</sup> [L<sub>☉</sub>]</td><td>2.86 × 10<sup>9</sup></td><td>fixed</td><td></td></tr><tr><td>2D geometry of the evolved stars</td><td>2D image</td><td>IRAC 3.6 µm<sup>a</sup></td><td>fixed</td><td></td></tr><tr><td>Vertical scale height</td><td>h<sub>z</sub> [pc]</td><td>538</td><td>fixed</td><td></td></tr><tr><td><strong>Young stars (100 Myr)</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Total intrinsic FUV luminosity</td><td>λL<sub>λ</sub><sup>young</sup> [L<sub>☉</sub>]</td><td>0.93 × 10<sup>9</sup></td><td>(0.7–1.9) × 10<sup>9</sup></td><td>0.96 × 10<sup>9</sup></td></tr><tr><td>2D geometry of the young stars</td><td>2D image</td><td>GALEX FUV<sup>b</sup></td><td>fixed</td><td></td></tr><tr><td>Vertical scale height</td><td>h<sub>z</sub> [pc]</td><td>190</td><td>fixed</td><td></td></tr><tr><td><strong>Ionising stars (10 Myr)</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Total intrinsic FUV luminosity</td><td>λL<sub>λ</sub><sup>ion.</sup> [L<sub>☉</sub>]</td><td>0.4 × 10<sup>8</sup></td><td>(0.4–2.6) × 10<sup>8</sup></td><td>2.0 × 10<sup>8</sup></td></tr><tr><td>2D geometry of the ionising stars</td><td>2D image</td><td>MIPS 24 µm<sup>c</sup></td><td>fixed</td><td></td></tr><tr><td>Vertical scale height</td><td>h<sub>z</sub> [pc]</td><td>190</td><td>fixed</td><td></td></tr><tr><td><strong>Interstellar Dust</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Total dust mass</td><td>M<sub>dust</sub> [M<sub>☉</sub>]</td><td>4.3 × 10<sup>7</sup></td><td>(3.2–5.4) × 10<sup>7</sup></td><td>4.7 × 10<sup>7</sup></td></tr><tr><td>2D geometry for the interstellar dust</td><td>2D image</td><td>MAGPHYS map<sup>d</sup></td><td>fixed</td><td></td></tr><tr><td>Vertical scale height</td><td>h<sub>z</sub> [pc]</td><td>238</td><td>fixed</td><td></td></tr></tbody></table>

**Notes.** We used an inclination angle of 77.5° to deproject the input images. Most parameters remain fixed during the SED modelling. For our 3 variable parameters we give the range of our parameter grid. A detailed description of all parameters and the initial guess values can be found in Appendix A. Luminosities are in units of bolometric solar luminosity $L_{☉}$. (a) The Einasto bulge was first subtracted from the image, see Sect. A.1. (b) The contribution from the evolved stars was first subtracted from the image, see Sect. A.2. (c) The contribution from the evolved stars was first subtracted from the image, see Sect. A.3. (d) The dust mass map derived from the pixel-by-pixel MAGPHYS energy balance SED models from Viaene et al. (2014).

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 M31 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 $M_{dust} = 4.7 \times 10^7 M_{☉}$. The intrinsic FUV luminosity of the young and ionising component is $\lambda L_{\lambda}^{\text{young}} = 9.6 \times 10^8 L_{☉}$ and

$\lambda L_{\lambda}^{\text{ion.}} = 2.0 \times 10^8 L_{☉}$, 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 $\chi^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 2–3 (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 M31, which lie between 0.2 and 0.4 $M_{☉} \text{yr}^{-1}$ (see e.g. Tabatabaei & Berkhuijsen 2010; HELGA IV; Rahmani et al. 2016; or Ford et al. 2013 (HELGA III)).
<figure><img src="image_1.png" /><figcaption>**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.</figcaption></figure>

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_{\odot}$ 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).

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(-\chi^2/2)$. For each parameter value we then add the probabilities and plot the PDFs in Fig. 2. We find that dust mass and $\lambda L_{\lambda}^{\text{young}}$ are fairly well constrained (the PDFs resemble a normal distribution). The $\lambda L_{\lambda}^{\text{ion}}$, 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

<table><caption>**Table 3.** Star formation rates for M 31, derived from synthetic broad band fluxes from our model SED.</caption><thead><tr><td>Band</td><td>SFR[$M_{\odot}$ yr$^{-1}$]</td><td>Recipe</td></tr></thead><tbody><tr><td>FUV</td><td>0.24</td><td>a, b</td></tr><tr><td>NUV</td><td>0.31</td><td>a, b</td></tr><tr><td>24 µm</td><td>0.23</td><td>c</td></tr><tr><td>70 µm</td><td>0.19</td><td>d</td></tr><tr><td>TIR (3–1000 µm)</td><td>0.92</td><td>a, b</td></tr><tr><td>FUV + 24 µm</td><td>0.31</td><td>e</td></tr></tbody></table>

**Notes.** For the *FUV* and *NUV* estimates, the intrinsic UV fluxes were used to avoid the problem of dust extinction. Conversions are based on the Kroupa & Weidner (2003) initial mass function. (a) Murphy et al. (2011); (b) Hao et al. (2011); (c) Rieke et al. (2009); (d) Calzetti et al. (2010); (e) Leroy et al. (2008).

optimal solution. The median (50th percentile) values that come from these PDFs are $M_{\text{dust}} = 4.72_{-0.39}^{+0.36} \times 10^7 M_{\odot}$, $\lambda L_{\lambda}^{\text{young}} = 10.5_{-1.6}^{+1.7} \times 10^8 L_{\odot}$ and $\lambda L_{\lambda}^{\text{ion}} = 1.4_{-0.6}^{+2.2} \times 10^8 L_{\odot}$. Here the upper and lower errors correspond to the 84th and 16th percentile of the PDFs, respectively.

## 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
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** PDFs of the three free parameters in our model optimization: the total dust mass $M_{\text{dust}}$ (left), the intrinsic FUV luminosity from the young stellar component $\lambda L_{\lambda}^{\text{young}}$ (middle) and the intrinsic FUV luminosity from the ionising stellar component $\lambda L_{\lambda}^{\text{ion.}}$ (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.</figcaption></figure>

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 $R_V = 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 \le 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.

At NIR and optical wavelengths ($\lambda > 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

<figure><img src="image_3.png" /><figcaption>**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).</figcaption></figure>

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
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

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):

$$ \text{residual} = \frac{\text{observation} - \text{model}}{\text{observation}}. \qquad (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 M31 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 M31 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 M51 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 $\beta$ in their pixel-by-pixel modified black body fits of M31. However, a variable $\beta$ 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 $\beta$ 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 M31. 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.

## 5. A face-on view of M31

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 M31 becomes one swirl of clumpy rings.

## 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).
<figure><img src="image_5.png" /><figcaption>Fig. 5. Observed (left) and model (right) colours tracing the FIR dust emission peak. A colour here corresponds to the flux ratio: log(flux₁/flux₂), that is, blue colours indicate a high flux ratio, red colours a low ratio.</figcaption></figure>

## 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 $L_l^{abs}$ is stored in every cell and at every wavelength. Integrating over wavelength gives the total absorbed luminosity $L^{abs}$ 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 (10⁸ 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 $L^{abs}$. 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 $L_{tot}^{abs}$. Including only evolved, young or ionising stellar populations yields $L_{evolved}^{abs}$, $L_{young}^{abs}$ and $L_{ion}^{abs}$, respectively. $L_{ion}^{abs}$ 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 $L_{ion,SF}^{abs}$ (as opposed to energy absorption due to diffuse ISM dust, $L_{ion}^{abs}$).

To test the sampling of the radiation field in our simulations, we plot the sum of $L_{evolved}^{abs}$, $L_{young}^{abs}$ and $L_{ion}^{abs}$ 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

$$F_{unev.}^{abs} = \frac{L_{young}^{abs} + L_{ion.}^{abs} + L_{ion,SF}^{abs}}{L_{evolved}^{abs} + L_{young}^{abs} + L_{ion.}^{abs} + L_{ion,SF}^{abs}}, \quad (2)$$
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Face-on model images of M 31 in the FUV, 3.6, 8.0, 24, 100 and 250 µm. Colour scale is different 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.</figcaption></figure>

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 $F_{ev}^{abs} = 1 - F_{unev}^{abs}$.

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 $F_{unev}^{abs} = 9\%$ 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 $F_{unev}^{abs}$ 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 $F_{unev}^{abs}$ 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 $F_{unev}^{abs}$ 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 $F_{unev}^{abs}$ 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 $F_{unev}^{abs}$, 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
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Left: 2D histogram comparing the integrated absorbed luminosity per dust cell including all stellar components $L_{\text{tot}}^{\text{abs}}$ 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: $F_{\text{unev.}}^{\text{abs}}$. The histogram is weighted by mass fraction of each dust cell and normalised. The red line indicates the median $F_{\text{unev.}}^{\text{abs}} = 0.10$. Right: 2D histogram showing the radial distribution of $F_{\text{unev.}}^{\text{abs}}$. 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.</figcaption></figure>

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 M31

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 $F_{\text{unev.}}$. Naturally, $F_{\text{evolved}}$ and $F_{\text{unev.}}$ 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 $F_{\text{unev.}} + F_{\text{evolved}} = 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_{\text{\lambda,evolved}}$ and $F_{\text{\lambda,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_{\text{\lambda}}$. De Looze et al. (2014) then define:

$$F'_{\text{\lambda,unev.}} = \frac{1}{2} \frac{S_{\text{\lambda,unev.}} + (S_{\text{\lambda,total}} - S_{\text{\lambda,evolved}})}{S_{\text{\lambda,total}}}, \quad (3)$$

and

$$F'_{\text{\lambda,evolved}} = \frac{1}{2} \frac{S_{\text{\lambda,evolved}} + (S_{\text{\lambda,total}} - S_{\text{\lambda,unev.}})}{S_{\text{\lambda,total}}}, \quad (4)$$

where $S_{\text{\lambda,tot}}$ is the dust emission heated by evolved, young and ionising stellar populations. $S_{\text{\lambda,evolved}}$ is the dust emission when only the evolved stellar components are included in the simulations, and $S_{\text{\lambda,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, $F'_{\text{\lambda,unev.}}$ and $F'_{\text{\lambda,evolved}}$ are the approximate heating fractions by unevolved and evolved stellar populations. They now add up to unity by definition.
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Dust heating fraction as a function of wavelength for the un-evolved stellar populations. The solid black line corresponds to the def-inition 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.</figcaption></figure>

This should be compared with the case where the dust heating fractions are not adjusted, that is

$$F_{\lambda,unev.} = \frac{S_{\lambda,unev.}}{S_{\lambda,unev.} + S_{\lambda,evolved}}, \qquad (5)$$

and

$$F_{\lambda, evolved} = \frac{S_{\lambda, evolved}}{S_{\lambda, unev.} + S_{\lambda, evolved}}. \qquad (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).

The general shape of the $F'_{\lambda,unev.}$ curves is qualitatively the same, however. At almost all wavelengths the unevolved stellar populations are not the main heating source in M 31. From 10–50 µm, $F'_{\lambda,unev.}$ is mostly constant at 40–50%, with several peaks. These peaks correspond to stochastically 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, $F'_{\lambda,unev.}$ 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 $F'_{\lambda,unev.}$ 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.

Going back to the global, aggregate heating fractions (not the approximate, wavelength-dependent quantities), we find that $F_{unev.} = 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 $F'_{\lambda,unev.}$ 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 data<fnref n="3" /> – 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) $F_{unev.}$. 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 $\tau$ (Kendall 1938; Kendall & Gibbons 1990). This is a non-parametric correlation coefficient where $\tau = \pm 1$ denotes a perfect (anti-)correlation, and $\tau = 0$ means there is no correlation.

In Fig. 10, we plot $F_{unev.}$ against $NUV - r$ colour, commonly used in the context of galaxy evolution as an observable tracer

<footnote n="3">³ For 512 wavelengths and $2 \times 10^7$ photon packages per wavelength, running one model takes about 4000 CPU hours and 15 GB of RAM.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Maps of the dust heating fraction by unevolved stellar populations $F'_{\lambda,unev.}$ 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.</figcaption></figure>

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 ($\tau = 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 $F_{unev.}$ 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.

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\alpha$ 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 ($L_{unev.}$). This quantity is similar to the $L^{abs}$ parameters used in Sect. 6.1, but now projected to the 2D pixels on the sky.

However, the trends we observe for $L_{unev.}$ are less convincing then for $F_{unev.}$: NUV-r is still the strongest correlation ($\tau = 0.21$), but still weak and clearly bimodal. For the 250/350 colour, the trend is weak ($\tau = 0.19$) and the colour range is narrow. Again, the 160/250 colour ($\tau = 0.07$) has a large variance and is also not suitable as a predictor for $L_{unev.}$ 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 $F_{unev.}$ There is a highly scattered positive trend with $M_{dust}$ ($\tau = 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 ($\tau = 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, $L_{dust}$, was found either ($\tau = -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 ($\tau = -0.34$) with stellar mass ($M_\star$). 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 $F_{unev.}$ 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 ($M_{dust}/M_\star$). This quantity gives an idea of the dust content, relative to the stellar mass. High $M_{dust}/M_\star$ values point at dusty regions and those can indeed be associated with high $F_{unev.}$ ($\tau = 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 ($\tau = 0.51$). This correlation is likely to be intrinsic as the trends with SFR or $M_\star$ separately show much more scatter. This trend is not completely unexpected, since we already found a positive link between $F_{unev.}$ and NUV-r colour, and this colour correlates with sSFR in M 31 (HELGA IV). sSFR traces the hardness of the
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Density plots of the dust heating fraction from unevolved stellar populations $F_{unev.}$ vs. $NUV - r$ and two commonly used FIR colours. Only pixels within a galactocentric radius of 18 kpc are shown. Kendall's correlation coefficient $\tau$ 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.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Density plots of the projected dust absorption luminosity from unevolved stellar populations $L_{unev.}$ vs. $NUV - r$ and two commonly used FIR colours. Only pixels within a galactocentric radius of 18 kpc are shown. Kendall's correlation coefficient $\tau$ 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.</figcaption></figure>

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 M51. We over-plot their data points and best fit curve in the central, bottom-row panel of Fig. 12. In M51, higher fractions of $F_{unev.}$ can be achieved compared to Andromeda. The regions of M31 follow the average relation for the M51 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 M31 pixels is about 3.5 times smaller than for the M51 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 M31.

The bottom line is that extensive galaxy parameters such as $M_{dust}$, $L_{dust}$, $M_{star}$ 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 $M_{dust}/M_\star$ 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 M31 is fitted well and the resulting attenuation curve is consistent with observations, but with a broader UV bump. We are able to constrain 2 of our 3 free parameters: the dust mass and the luminosity of the young stellar component. The intrinsic luminosity of the ionising stellar population is more difficult to constrain.
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Density plots of the dust heating fraction from unevolved stellar populations $F_{\text{unev.}}$ 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 $\tau$ 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.</figcaption></figure>

- 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, NUV - r 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 NUV - r) 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.

**Acknowledgements.** 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

Aniano, G., Draine, B. T., Calzetti, D., et al. 2012, ApJ, 756, 138

Athanassoula, E., & Beaton, R. L. 2006, MNRAS, 370, 1499

Baes, M., & Camps, P. 2015, Astron. Comput., 12, 33

Baes, M., Fritz, J., Gadotti, D. A., et al. 2010, A&A, 518, L39

Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, 22

Bailey, N. D., & Basu, S. 2014, ApJ, 780, 40

Barmby, P., Ashby, M. L. N., Bianchi, L., et al. 2006, ApJ, 650, L45

Battisti, A. J., Calzetti, D., & Chary, R.-R. 2016, ApJ, 818, 13

Beaton, R. L., Majewski, S. R., Guhathakurta, P., et al. 2007, ApJ, 658, L91

Bendo, G. J., Wilson, C. D., Pohlen, M., et al. 2010, A&A, 518, L65

Bendo, G. J., Boselli, A., Dariush, A., et al. 2012, MNRAS, 419, 1833

Bendo, G. J., Baes, M., Bianchi, S., et al. 2015, MNRAS, 448, 135

Bianchi, S. 2007, A&A, 471, 765

Bianchi, S. 2008, A&A, 490, 461

Block, D. L., Bournaud, F., Combes, F., et al. 2006, Nature, 443, 832

Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111

Boquien, M., Buat, V., Boselli, A., et al. 2012, A&A, 539, A145

Boquien, M., Buat, V., & Perret, V. 2014, A&A, 571, A72

Boquien, M., Calzetti, D., Aalto, S., et al. 2015, A&A, 578, A8

Boquien, M., Kennicutt, R., Calzetti, D., et al. 2016, A&A, 591, A6

Boselli, A., Ciesla, L., Buat, V., et al. 2010, A&A, 518, L61

Boselli, A., Ciesla, L., Cortese, L., et al. 2014, A&A, 540, A54

Bressan, A., Fagotto, F., Bertelli, G., et al. 2001, VizieR Online Data Catalog: VI, 196

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Calzetti, D., Wu, S.-Y., Hong, S., et al. 2010, ApJ, 714, 1256

Camps, P., & Baes, M. 2015, Astron. Comput., 9, 20

Camps, P., Baes, M., & Safly, W. 2013, A&A, 560, A35

Camps, P., Misselt, K., Bianchi, S., et al. 2015, A&A, 580, A87

Chabrier, G. 2003, PASP, 115, 763

Ciesla, L., Boquien, M., Boselli, A., et al. 2014, A&A, 565, A128

Corbelli, E., Lorenzoni, S., Walterbos, R., Braun, R., & Thilker, D. 2010, A&A, 511, A89

Corbelli, E., Bianchi, S., Cortese, L., et al. 2012, A&A, 542, A32

Cortese, L., Boselli, A., Franzetti, P., et al. 2008, MNRAS, 386, 1157

Cortese, L., Ciesla, L., Boselli, A., et al. 2012, A&A, 540, A52

Cortese, L., Fritz, J., Bianchi, S., et al. 2014, MNRAS, 440, 942

Courteau, S., Widrow, L. M., McDonald, M., et al. 2011, ApJ, 739, 20

da Cunha, E., Eminian, C., Charlot, S., & Blaizot, J. 2010, MNRAS, 403, 1894

Daddi, E., Dickinson, M., Chary, R., et al. 2005, ApJ, 631, L13

Dalcanton, J. J., Fouesneau, M., Hogg, D. W., et al. 2015, ApJ, 814, 3

Dale, D. A., Aniano, G., Engelbracht, C. W., et al. 2012, ApJ, 745, 95

De Geyter, G., Baes, M., Camps, P., et al. 2014, MNRAS, 441, 869

De Geyter, G., Baes, M., De Looze, I., et al. 2015, MNRAS, 451, 1728

De Looze, I., Baes, M., Bendo, G. J., et al. 2012a, MNRAS, 427, 2797

De Looze, I., Baes, M., Fritz, J., & Verstappen, J. 2012b, MNRAS, 419, 895

De Looze, I., Fritz, J., Baes, M., et al. 2014, A&A, 571, A69

de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Jr., H. G., et al. 1991, in Third Reference Catalogue of Bright Galaxies (Springer)

Devereux, N. A., Price, R., Wells, L. A., & Duric, N. 1994, AJ, 108, 1667

Dierickx, M., Blecha, L., & Loeb, A. 2014, ApJ, 788, L38

Dong, H., Li, Z., Wang, Q. D., et al. 2014, ApJ, 785, 136

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Dunne, L., Eales, S., Edmunds, M., et al. 2000, MNRAS, 315, 115

Eales, S., Smith, M. W. L., Auld, R., et al. 2012, ApJ, 761, 168

Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320

Ford, G. P., Gear, W. K., Smith, M. W. L., et al. 2013, ApJ, 769, 55

Foyle, K., Natale, G., Wilson, C. D., et al. 2013, MNRAS, 432, 2182

Fritz, J., Gentile, G., Smith, M. W. L., et al. 2012, A&A, 546, A34

Galametz, M., Madden, S. C., Galliano, F., et al. 2010, A&A, 518, L55

Galametz, M., Kennicutt, R. C., Albrecht, M., et al. 2012, MNRAS, 425, 763

Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279

Gordon, K. D., Bailin, J., Engelbracht, C. W., et al. 2006, ApJ, 638, L87

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Groves, B., Dopita, M. A., Sutherland, R. S., et al. 2008, ApJS, 176, 438

Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892

Groves, B. A., Schinnerer, E., Leroy, A., et al. 2015, ApJ, 799, 96

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124

Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852

Hughes, T. M., Baes, M., Fritz, J., et al. 2014, A&A, 565, A4

Indebetouw, R., Whitney, B. A., Johnson, K. E., & Wood, K. 2006, ApJ, 636, 362

Jarrett, T. H., Masci, F., Tsai, C. W., et al. 2013, AJ, 145, 6

Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62

Kendall, M., & Gibbons, J. 1990, in Rank correlation methods, Charles Griffin Book Series (London: E. Arnold)

Kendall, M. G. 1938, Biometrika, 30, 81

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kennicutt, Jr., R. C., Hao, C.-N., Calzetti, D., et al. 2009, ApJ, 703, 1672

Kirkpatrick, A., Pope, A., Alexander, D. M., et al. 2012, ApJ, 759, 139

Kirkpatrick, A., Calzetti, D., Kennicutt, R., et al. 2014, ApJ, 789, 130

Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076

Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782

MacLachlan, J. M., Matthews, L. D., Wood, K., & Gallagher, J. S. 2011, ApJ, 741, 6

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1

Mattsson, L., Gomez, H. L., Andersen, A. C., et al. 2014, MNRAS, 444, 797

McConnachie, A. W., Irwin, M. J., Ferguson, A. M. N., et al. 2005, MNRAS, 356, 979

Mentuch Cooper, E., Wilson, C. D., Foyle, K., et al. 2012, ApJ, 755, 165

Misiriotis, A., Popescu, C. C., Tuffs, R., & Kylafis, N. D. 2001, A&A, 372, 775

Miville-Deschênes, M.-A., & Lagache, G. 2005, ApJS, 157, 302

Montalto, M., Seitz, S., Riffeser, A., et al. 2009, A&A, 507, 283

Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67

Murray, N. 2011, ApJ, 729, 133

Natale, G., Tuffs, R. J., Xu, C. K., et al. 2010, ApJ, 725, 955

Natale, G., Popescu, C. C., Tuffs, R. J., et al. 2015, MNRAS, 449, 243

Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJ, 278, L1

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Planck Collaboration I. 2011, A&A, 536, A1

Planck Collaboration Int. XXV. 2015, A&A, 582, A28

Planck Collaboration Int. I. 2016, A&A, 586, A132

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Popescu, C. C., Misiriotis, A., Kylafis, N. D., Tuffs, R. J., & Fischera, J. 2000, A&A, 362, 138

Popescu, C. C., Tuffs, R. J., Völk, H. J., Pierini, D., & Madore, B. F. 2002, ApJ, 567, 221

Popescu, C. C., Tuffs, R. J., Dopita, M. A., et al. 2011, A&A, 527, A109

Rahmani, S., Lianou, S., & Barmby, P. 2016, MNRAS, 456, 4128

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015, A&A, 582, A121

Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25

Rieke, G. H., Alonso-Herrero, A., Weiner, B. J., et al. 2009, ApJ, 692, 556

Rosenfield, P., Johnson, L. C., Girardi, L., et al. 2012, ApJ, 755, 131

Roussel, H. 2013, PASP, 125, 1126

Safly, W., Baes, M., & Camps, P. 2014, A&A, 561, A77

Saglia, R. P., Fabricius, M., Bender, R., et al. 2010, A&A, 509, A61
Schneider, R., Bianchi, S., Valiante, R., Risaliti, G., & Salvadori, S. 2015, A&A, 579, A60

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41

Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103

Skibba, R. A., Engelbracht, C. W., Dale, D., et al. 2011, ApJ, 738, 89

Smith, M. W. L., Eales, S. A., Gomez, H. L., et al. 2012a, ApJ, 756, 40

Smith, M. W. L., Gomez, H. L., Eales, S. A., et al. 2012b, ApJ, 748, 123

Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63

Tabatabaei, F. S., & Berkhuijsen, E. M. 2010, A&A, 517, A77

Tamm, A., Tempel, E., Tenjes, P., Tihhonova, O., & Tuvikene, T. 2012, A&A, 546, A4

Tempel, E., Tamm, A., & Tenjes, P. 2010, A&A, 509, A91

Tempel, E., Tuvikene, T., Tamm, A., & Tenjes, P. 2011, A&A, 526, A155

Thilker, D. A., Hoopes, C. G., Bianchi, L., et al. 2005, ApJ, 619, L67

Városi, F., & Dwek, E. 1999, ApJ, 523, 265

Viaene, S., Fritz, J., Baes, M., et al. 2014, A&A, 567, A71

Viaene, S., De Geyter, G., Baes, M., et al. 2015, A&A, 579, A103

Viaene, S., Baes, M., Bendo, G., et al. 2016, A&A, 586, A13

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Whelan, D. G., Johnson, K. E., Whitney, B. A., Indebetouw, R., & Wood, K. 2011, ApJ, 729, 111

Witt, A. N., & Gordon, K. D. 1996, ApJ, 463, 681

Witt, A. N., & Gordon, K. D. 2000, ApJ, 528, 799

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Wu, H., Zhu, Y.-N., Cao, C., & Qin, B. 2007, ApJ, 668, 87

Xilouris, E. M., Byun, Y. I., Kylafis, N. D., Paleologou, E. V., & Papamastorakis, J. 1999, A&A, 344, 868

York, D. G., Adelman, J., Anderson, Jr., J. E., et al. 2000, AJ, 120, 1579

Ysard, N., Abergel, A., Ristorcelli, I., et al. 2013, A&A, 559, A133

Ysard, N., Köhler, M., Jones, A., et al. 2015, A&A, 577, A110
# Appendix A: Model components

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.

## A.1. Evolved stellar populations

The 3D structure of the evolved stellar populations (>100 Myr<fnref n="4" />) is a necessary input for our model. We simplify this distribution to two components: a disk and a bulge. For M31, 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 M31 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 M31 (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 M31. 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$:

$$S(r) = S_c \exp \left\{ -d_N \left[ \left( \frac{r}{a_c} \right)^{1/N} - 1 \right] \right\}. \quad (A.1)$$

This prescription has two free parameters: the concentration radius of $a_c = 1.155$ kpc which contains half of the total brightness, and a structure parameter $N = 2.7$. The numerical constant $d_N$ is only dependent on $N$ and ensures that half of the brightness is containted in $a_c$. An intrinsic flattening in the vertical direction of $q = 0.72$ was imposed in the sense that

$$S(x, y, z) = \frac{1}{q} S_c \left( \sqrt{x^2 + y^2 + \frac{z^2}{q^2}} \right). \quad (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 M31 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 M32 (see Sect. 2). This correction is unlikely to affect the radiation field in the nucleus of M31 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 M31. 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 M31. 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.

## 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 $10^6-10^7$ 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 M31 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 $S_{FUV}^{young}$, the surface brightness from young stellar populations in the FUV:

$$S_{FUV}^{young} = S_{FUV} - \alpha_{FUV} S_{IRAC 3.6}. \quad (A.3)$$

In HELGA III this was applied to M31 and $\alpha_{FUV} = 8.0 \times 10^{-4}$ was found, based on the same FUV and IRAC 3.6 µm observations of M31 as used in this work. Note that this correction only

<footnote n="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.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Effect of galaxy thickness on residuals in the FUV band (top row) and at 250 µm (bottom row). The first column shows residuals for the original thickness, middle and right columns for 2/3 and 1/2 of that thickness, respectively.</figcaption></figure>

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.

### 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\alpha$ emission. Devereux et al. (1994) mapped a large part of Andromeda in $H\alpha$. 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 $S_{FUV}^{ion}$, the surface brightness of the ionising stellar component in the FUV:

$$S_{FUV}^{ion} = S_{24} - \alpha_{24} \times S_{IRAC3.6} \qquad (A.4)$$

This was also applied to M 31 in HELGA III. They found that $\alpha_{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 ($\log C = 6$), pressure of the surrounding ISM ($P_0 = 10^{12} \text{ K/m}^3$) 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_{\odot} yr^{-1}$. We take the template luminosity in the FUV band and scale this to $0.2 M_{\odot} yr^{-1}$, which is a good approximation of the SFR in M31 (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.

### 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 $A_{FUV}$ 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 $A_{FUV}$, 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 M31 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 M31. 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 $\le 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 on 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

## Appendix B: Observed fluxes

<table><caption>**Table B.1.** Integrated fluxes for M31 adopted in this paper, listed by increasing central wavelength.</caption><thead><tr><th>Instrument and band</th><th>Flux (Jy)</th><th>Error (Jy)</th></tr></thead><tbody><tr><td>GALEX FUV</td><td>1.325</td><td>0.066</td></tr><tr><td>GALEX NUV</td><td>2.445</td><td>0.073</td></tr><tr><td>SDSS u</td><td>17.77</td><td>0.36</td></tr><tr><td>SDSS g</td><td>80.9</td><td>1.6</td></tr><tr><td>SDSS r</td><td>177.2</td><td>3.5</td></tr><tr><td>SDSS i</td><td>270.8</td><td>5.4</td></tr><tr><td>SDSS z</td><td>328.5</td><td>6.6</td></tr><tr><td>WISE 1</td><td>280.8</td><td>6.7</td></tr><tr><td>IRAC 1</td><td>280</td><td>23</td></tr><tr><td>IRAC 2</td><td>158</td><td>11</td></tr><tr><td>WISE 2</td><td>155</td><td>4.3</td></tr><tr><td>IRAC 3</td><td>215</td><td>48</td></tr><tr><td>IRAC 4</td><td>209</td><td>35</td></tr><tr><td>WISE 3</td><td>203</td><td>9.1</td></tr><tr><td>IRAS 12 µm</td><td>181</td><td>24</td></tr><tr><td>WISE 4</td><td>160</td><td>9.1</td></tr><tr><td>MIPS 24 µm</td><td>118</td><td>47</td></tr><tr><td>IRAS 25 µm</td><td>126</td><td>17</td></tr><tr><td>IRAS 60 µm</td><td>820</td><td>110</td></tr><tr><td>MIPS 70 µm</td><td>1042</td><td>104</td></tr><tr><td>PACS 100 µm</td><td>3389</td><td>339</td></tr><tr><td>IRAS 100 µm</td><td>3570</td><td>470</td></tr><tr><td>PACS 160 µm</td><td>7397</td><td>739</td></tr><tr><td>SPIRE 250 µm</td><td>5858</td><td>410</td></tr><tr><td>SPIRE 350 µm</td><td>3061</td><td>214</td></tr><tr><td>Planck 350 µm</td><td>2900</td><td>300</td></tr><tr><td>SPIRE 500 µm</td><td>1317</td><td>92</td></tr><tr><td>Planck 550 µm</td><td>980</td><td>100</td></tr><tr><td>Planck 850 µm</td><td>291</td><td>14</td></tr></tbody></table>

(see Sect. A.3). The total observed dust mass in M31 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 $\kappa_{350 \mu m} = 0.454 m^2/kg$. In the Jones et al. (2013) model, the average $\kappa_{350 \mu m}$ becomes $0.305 m^2/kg$. This enhances the total dust mass by a factor of 1.5. However, we note that this average $\kappa_{350 \mu 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.

Table B.1 lists the integrated fluxes for M31 used in this paper. These are the integrated fluxes from the region we consider in our model. This an ellipse with centre at $\alpha = 00:43:06.28$ and $\delta = +41:21:12.22$, semi-major axis of $1.425^\circ$, semi-minor axis of $0.400^\circ$ and a position angle of $38.1^\circ$. 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).