Deriving the intrinsic properties of M51 with radiative transfer models

ABSTRACT

A quantitative derivation of the intrinsic properties of galaxies related to their fundamental building blocks, gas, dust, and stars is essential for our understanding of galaxy evolution. A fully self-consistent derivation of these properties can be achieved with radiative transfer (RT) methods that are constrained by panchromatic imaging observations. Here, we present an axi-symmetric RT model of the UV-optical-FIR/submm spectral and spatial energy distribution of the face-on spiral galaxy M51. The model reproduces reasonably well the azimuthally averaged radial profiles derived from the imaging data available for this galaxy, from GALEX, Sloan Digital Sky Survey, 2MASS, Spitzer, and Herschel. We model the galaxy with three distinct morphological components: a bulge, an inner disc, and a main disc. We derive the length parameters of the stellar emissivity and of the dust distribution. We also derive the intrinsic global and spatially resolved parameters of M51. We find a faint ‘outer disc’ bridging M51 with its companion galaxy M51b. Finally, we present and discuss an alternative model, with dust properties that change within the galaxy.

1 INTRODUCTION

Dust is prevalent within the interstellar medium (ISM) of star-forming galaxies (Trumpler 1930), and, although it constitutes only 1 per cent of the mass of this ISM (Greenberg 1963), it plays an important role in the thermodynamic balance of galaxies (Montier & Giard 2004; Giard et al. 2008; Natale et al. 2010; Popescu & Tuffs 2010), since it regulates the cooling and heating mechanisms of the interstellar gas (Dwek & Werner 1981).

The interstellar dust also affects our view of galaxies. This is because stellar photons, primarily in the ultraviolet (UV)/optical range, are continuously absorbed and scattered by the interstellar dust. Photons that are scattered can travel long distances within the ISM before they are scattered again, absorbed or leave the galaxy. Due to scattering, what would be otherwise an isotropic process related to the emission of stellar light in galaxies, becomes highly anisotropic. As such, along any line of sight through a galaxy, the stellar light may be diminished due to extinction, but also amplified by stellar photons being scattered into the line of sight. Thus, for a extended distribution of emitters and absorbers, the complex effect of absorption and scattering processes is what we call attenuation.

The attenuation is dependent not only on the optical properties of the dust grains (as is the case for the extinction processes), but also on the relative distribution of stars and dust, and on the orientation of a galaxy with respect to the observer (Byun, Freeman & Kylafis 1994; Bianchi, Ferrara & Giovanardi 1996; Kuchinski et al. 1998; Ferrara et al. 1999; Baes & Dejonghe 2001; Fischera & Dopita 2004; Pierini et al. 2004; Tuffs et al. 2004; Natale et al. 2015). Dust attenuation affects the observed spatially integrated stellar luminosity of a galaxy, but also its surface brightness distribution and the geometric parameters derived from surface-brightness photometry, like scale length of discs, effective radii of bulges and Sérsic indexes: Byun et al. (1994), Evans (1994), Cunow (2001), Möllenhoff, Popescu & Tuffs (2006), Gadotti, Baes & Falony (2010), Pastrav et al. (2013a, b), and Savchenko et al. (2023).

Due to dust attenuation the observed UV/optical images of star-forming galaxies are far from providing a direct tracer of the constituent stellar populations distributed throughout the galaxy. In order to derive the intrinsic distribution of stars of all ages, and thus to understand how stars form and evolve in galaxies, one also needs to jointly derive the distribution of dust in galaxies (Popescu & Tuffs 2010). This can be done by considering both the absorption and emission properties of interstellar dust. Thus, the stellar light that is absorbed by dust is re-emitted in the infrared (IR) and the energy balance between direct and re-radiated stellar light is an important constraint on the modelling of transfer of radiation in galaxies. In particular dust emission contains unique information on both the spatial and spectral energy distribution (SED) of the energy density of the radiation fields that heat the dust grains, and on the dust properties and dust opacity. Accounting for both dust attenuation and dust emission is thus needed to provide a self-consistent model of a galaxy.

Galaxy models that are based on radiative transfer (RT) methods (Kylafis & Xilouris 2005; Steinacker, Baes & Gordon 2013) and constrained by panchromatic imaging observations are the most powerful tools in quantifying the attenuation of stellar light in galaxies and in deriving the intrinsic parameters of galaxies (Popescu 2021). This type of modelling was initially performed on edge-on galaxies, since in this orientation it is possible to directly constrain the vertical distribution of stars and dust from imaging observations. The first RT modelling of an edge-on galaxy was performed in the optical range, for NGC 891, by Kylafis & Bahcall (1987). This was the first time an RT model was applied in conjunction with realistic distributions of stars and dust (finite exponential discs as opposed to slabs or sandwich geometries – Disney, Davies & Phillipps 1989 or spherical configurations – Witt, Thronson & Capuano 1992; Witt & Gordon 1996, 2000). Fundamental work on characterizing the relative distributions of stars and dust and their statistical trends were later established in a series of papers dedicated to the modelling of edge-on galaxies in the optical and NIR by Xilouris et al. (1997, 1998, 1999).

The first panchromatic (UV/optical/FIR/submm) RT modelling of an edge-on galaxy was achieved for the same edge-on spiral NGC 891 (see above for the first RT modelling) by Popescu et al. (2000). Further work on modelling edge-on spirals was performed by Misiriotis et al. (2001), Popescu et al. (2004), Bianchi (2008), Baes et al. (2010), Popescu et al. (2011), MacLachlan et al. (2011), Robitaille et al. (2012), De Geyter et al. (2015), Mosenkov et al. (2016, 2018), Popescu et al. (2017), and Natale et al. (2021).

More recently, non-edge-on galaxies have also started to be modelled with RT methods. There are two approaches to this type of modelling. The first approach is to use non-axi-symmetric RT models: De Looze et al. (2014), Viaene et al. (2017), Williams et al. (2019), Verstocken et al. (2020), and Nersesian et al. (2020a, b). Non-axi-symmetric models are, in principle, ideal to capture the detailed geometry seen face-on. However, due to computational limitations, they are implemented to only fit the integrated SED, with typically only four free parameters that describe the total luminosity output for stars and the dust mass, with the geometry being fixed. The second approach is to use axi-symmetric RT models (Thirlwall et al. 2020; Rushton et al. 2022) and fit not only the integrated SED, but also the geometry of stars and dust, albeit with the assumption of axi-symmetry. So far only non-interacting face-on spiral galaxies have been modelled with axi-symmetric models (M33; Thirlwall et al. 2020, NGC 628; Rushton et al. 2022). The goal of this paper is to extend the applicability of the model to a non-isolated spiral. We will apply the model to M51, since, despite its interaction with a close companion, most of the extent of the galaxy has regular spiral structure, and is not perturbed by the interaction, thus making it suitable for the axi-symmetry approximation.

First discovered by Charles Messier in 1773, M51 is a system of interacting galaxies consisting of M51a (=NGC 5194), the ‘Whirlpool Galaxy’ , a grand-design spiral of Hubble type Sbc, and the smaller lenticular galaxy M51b (= NGC 5195), of Hubble type SB0. While interaction between the two galaxies is evident from the ‘bridge’ linking them, there is still debate of whether the galaxies have interacted on a single fly-by collision (Toomre & Toomre 1972; Durrell et al. 2003) or by multiple close passes (Salo & Laurikainen 2000; Dobbs et al. 2010). Despite the interaction, as mentioned before, M51a exhibits a good deal of azimuthal symmetry. This, in conjunction with the close proximity (⁠|$\text{D}=8.58\, \text{Mpc}$|⁠; McQuinn et al. 2016) and the wealth of detailed observational data, makes M51a ideal for studying the energy balance between absorption by dust and re-emission in the far-IR for a non-isolated spiral galaxy that is seen not far from face-on. In the following, we will refer to M51a as M51.

In this work, we will utilize the axisymmetric RT models developed in Popescu et al. (2000), Tuffs et al. (2004), and Popescu et al. (2011), which have been successfully used to account for the properties of star-forming galaxies derived from statistical samples (e.g. van der Giessen et al. 2022), as well as for the detailed spatial and SED of individual nearby galaxies (Popescu et al. 2000; Thirlwall et al. 2020; Rushton et al. 2022) and of the Milky Way (Popescu et al. 2017; Niederwanger et al. 2019; Natale et al. 2021). These models describe the geometries of the stellar emissivity and dust opacity from the UV to submm using parameterized analytic functions. The models consider the absorption and anisotropic scattering of stellar photons with dust grains of various sizes and chemical compositions, with the optical constants taken from Draine & Lee (1984) and Draine & Li (2007). The models also consider the heating of the dust grains and the re-emission in the IR, including explicit calculations for stochastic heating.

Previous work on modelling M51 was done by De Looze et al. (2014) and Nersesian et al. (2020b) using non-axisymmetric RT models. As discussed above, these models have the inherent limitation that they cannot self-consistently derive the geometry of emitters and absorbers. In this study, we aim to model M51, by self-consistently deriving the geometry of stars and dust through fitting the detailed azimuthally averaged UV/optical/FIR/submm surface brightness profiles of M51. As we will show here, our axi-symmetric model produces a good fit to the imaging data within a radial distance of 7 kpc from the centre of M51, within which the galaxy preserves azimuthal symmetry.

The paper is organized as follows: In Section 2, we present the imaging observations used for constraining the model, and the data-reduction steps undertaken in this study. In Section 3, we briefly describe the model and its geometrical components, with the fitting procedure detailed in Section 4. In Section 5, we present the results of the fits for the surface-brightness distributions and for the intrinsic parameters of M51, including the SFR, dust optical depth, dust mass, and the fractional contribution of the different stellar populations in heating the dust and the attenuation curve. We also present the results for the derived morphological components of M51. In Section 6, we compare our model for M51 with previous RT modelling of this galaxy and discuss an alternative to the standard model presented in this paper. We summarize our main results in Section 7.

2 OBSERVATIONAL DATA

In this section, we discuss the data set of M51 used in this study. To start with, basic parameters like inclination and position angle were difficult to determine with standard methods, and, because of this, there has been a lot of uncertainty in the value of these parameters. Determinations of the inclination angle vary from 15° to 42° (Danver 1942; Toomre & Toomre 1972; Tully 1974; Monnet, Paturel & Simien 1981; Considere & Athanassoula 1982; Bersier et al. 1994; Shetty et al. 2007; Hu, Shao & Peng 2013), while determinations for the PA vary from −10° to 36° (Burbidge & Burbidge 1964; Tully 1974; Monnet et al. 1981; Considere & Athanassoula 1982; Shetty et al. 2007; Tamburro et al. 2008; Walter et al. 2008; Hu et al. 2013). We adopted the most recent values from Hu et al. (2013). The distance, inclination, position angle and central coordinates used in this paper are listed in Table 1. We obtained panchromatic images from various missions, as summarized in Table 2.

All data were found on NASA/IPAC Extragalactic Database (NED),¹ except for the GALEX data, which was obtained from the Mikulski Archive for Space Telescopes (MAST).² Before analysis, the units of the surface brightness maps were converted into (MJy sr⁻¹).

2.1 Background removal

We performed a CoG analysis to determine the radius from the centre of M51 at which the emission of the galaxy falls below the background level. This was found to be 25 kpc. We also visually inspect each map to ensure that the emission associated with M51 does not extend beyond this point. The outer sampled radius (30 kpc) was chosen to ensure a large area for measuring the background. This sampling region extended beyond the limit of some maps, thus the full sampling area could not be used in those cases. However, even in the most severe cases, this was a negligible portion of the total sampling area. We then calculate the average background flux by taking the mean surface brightness through six annuli of equal width contained in the region |$25\, \text{kpc}\le R \le 30\, \text{kpc}$|⁠. Although M51 is nearly face-on with an inclination of 20|${_{.}^{\circ}}$|3, the interaction with M51b gives rise to some stripped material in the outer disc, making the outer isophotes more elongated than the inner ones. This change in isophotal shape was taken into account when defining the background. Thus, it was found that the best-fitting annuli for the background had an inclination of 40° (and the same position angle as the galaxy of 12°).

2.2 Foreground star removal

Foreground stars are visible in the UV/optical/NIR wavelengths and can influence the measured values of emission from M51. It is therefore important to mask these stars. For this, we first produce a median map where each pixel of the median map is equal to the median of the 3 × 3 pixels around it. A bias map is also produced, that is principally the same as the median map, but takes the median of the 21 × 21 pixels. This bias map is used as an additional check to mitigate against false positives when masking noisy data. Pixels with values that exceed both the median and bias maps are then masked. After this initial pass, we then manually inspect each mask for regions that may have been missed, or unmask false positives, such as H ii regions. H ii regions appear as point sources (similarly to foreground stars) and thus are often falsely masked by this automated process. These H ii regions were identified by eye and the catalogue of Lee, Hwang & Lee (2011). We also masked the companion galaxy M51b and its diffuse emission.

2.3 Convolutions

The Sloan Digital Sky Survey (SDSS), 2MASS, and Spitzer data have a higher resolution than that of our model, which would correspond to 50 pc at the distance of M51. We therefore degrade the resolution of these data to match the resolution of our model using a two-dimensional Gaussian kernel. The data at all the other wavelengths used in this study have lower resolutions, therefore, for these wavelengths, we degrade the resolution of our model to match the corresponding ones from the data. For this, we used the kernels of Aniano et al. (2011) to convolve the model with the PSF of the observational instrument.

2.4 Galaxy Evolution Explorer

The Galaxy Evolution Explorer (GALEX; Martin et al. 2005) was a NASA Small Explorer mission that performed an all-sky survey in the far-ultraviolet (FUV |$\sim 0.15\, \mu$|m) and near-ultraviolet (NUV |$\sim 0.22\, \mu$|m). Both FUV and NUV data are used in this study and were found on the MAST³ and made available by Gil de Paz et al. (2007). The NUV and FUV data have a pixel scale of 1.5 arcsec (62 pc) per pixel and an FWHM of 4.2 and 5.3 arcsec, respectively. The flux calibration error is |$\frac{\epsilon _{\text{cal}}}{F_\nu }$|  = 10  per cent following Morrissey et al. (2007).

2.5 Sloan Digital Sky Survey

The SDSS (York et al. 2000) has observed M51 with two sets of CCD arrays. Each array has an integration time of 51 s (York et al. 2000) capturing images through u, g, r, i, and z filters. In this study, we use the u, g, r, and i SDSS bands, with images obtained from NED/IPAC and processed by Brown et al. (2014). We do not use the z image as this wavelength is not included in our model fitting. Each band has a pixel scale of 0.396 arcsec (16 pc) and FWHM of 1 arcsec (Gunn et al. 2006). Each band has a 2 per cent rms photometric calibration in all bands (Aihara et al. 2011). The u band required a correction of –0.04 mag (Brown et al. 2014).

2.6 Two-Micron All-Sky Survey

The Two-Micron All-Sky Survey (2MASS; Jarrett et al. 2003) was performed in the NIR bands, J, H, and K_s, and observed M51 as part of the Large Galaxy Atlas Survey. We use the 2MASS data in the J and K_s bands provided on NED/IPAC, and processed by Brown et al. (2014). We do not use the H-band data as this wavelength is not included in our model. Both J and K_s bands have a pixel scale of 1 arcsec (41 pc) and an FWHM of 2.5 arcsec (Skrutskie et al. 2006). The calibration errors for J and K_s are 2 and 3 per cent (Jarrett et al. 2003), respectively.

2.7 Spitzer

2.7.1 Infrared Array Camera

The Spitzer Space Telescope contains the Infrared Array Camera (IRAC). IRAC is a four-channel camera that images the 3.6-, 4.5-, 5.8-, and 8.0-μm bands simultaneously. In this study, we use all four bands. The images were found on NED/IPAC and were processed by Brown et al. (2014). Each instrument has a pixel size of 0.75 arcsec (31 pc). The 3.6-, 4.5-, 5.8-, and 8.0-μm data have an FWHM of 1.66, 1.72, 1.88, and 1.98 arcsec, respectively (Fazio et al. 2004). The flux calibration error for the IRAC data are 3 per cent (Fazio et al. 2004).

2.7.2 Multiband Imaging Photometer for Spitzer

The Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) conducted imaging in three wavebands: 24, 70, and 160 μm. In this study, we only use the 24-μm data, since at 70 and 160 μm we use the higher quality data from Herschel. The data were processed by Brown et al. (2014) and made available at NED/IPAC. MIPS24 has a pixel scale of 1.5 arcsec (62 pc) and an FWHM of 6 arcsec (Rieke et al. 2004). The calibration error is 1 per cent (Rieke et al. 2004).

2.8 Herschel

2.8.1 Photo Array Camera & Spectrometer

M51 is observed by the Photo Array Camera & Spectrometer (PACS; Poglitsch et al. 2010) as part of the Very Nearby Galaxy Survey (VNGS) at the 70- and 160-μm bands. Both bands are utilized in this study and images were obtained from NED/IPAC. The pixel scale of the 70- and 160-μm data is 1.4 (58 pc) and 2.85 arcsec (118 pc), respectively. The full width at half-maximum (FWHM) of the 70- and 160-μm data are 5.6 and 11.4 arcsec, respectively (Bendo et al. 2012).

2.8.2 Spectral and Photometric Imaging Receiver

The Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) observed M51 as part of the VNGS at the 250-, 350-, and 500-μm bands. All three bands were used in this study and images were obtained from NED/IPAC and processed by Bendo et al. (2012). The pixel scale of the 250-, 350-, and 500-μm data is 6 (250 pc), 8 (333 pc), and 12 arcsec (500 pc), respectively. The FWHM of SPIRE data for 250, 350, and 500 μm are 18.2, 24.5, and 36 arcsec, respectively (Bendo et al. 2012). Flux calibration error is accurate to |$\pm 15~{{\ \rm per\ cent}}$|(Griffin et al. 2010).

Examples of imaging data used in this analysis are given in the left-hand panels of Fig. 1. They show 15 × 15 kpc images of M51 in the GALEX, NUV, 2MASS K band, and the SPIRE 500 μm.

2.9 Photometry

The observed images were corrected for Galactic extinction using a value of E(B − V) = 0.0308 (Schlafly & Finkbeiner 2011) and R_V = 3.1 (Fitzpatrick 1999). Azimuthally averaged surface brightness profiles and spatially integrated flux densities F_ν were derived using curve of growth (CoG) photometry. Each map was segmented into 250 annuli of a fixed width in which the average brightness was calculated, with the corresponding radius taken as the midpoint of the annulus. The elliptical annuli were defined by the inclination and position angle of M51, as listed in Table 1. Examples of azimuthally averaged images (NUV, K, and 500 μm) produced in this way are given in the middle panels of Fig. 1. Comparison between azimuthally averaged and corresponding original observed images shows that the assumption of axi-symmetry seems to hold up to ∼7 kpc, which we take to represent the limit of applicability of our axi-symmetric model. We will refer to this limit as R_l.

The photometric errors due to background fluctuations, calibration, and configuration noise were derived using the procedure from Thirlwall et al. (2020).

3 THE MODEL

Our model for M51 follows the axisymmetric RT models of Popescu et al. (2011, hereafter PT11), which describe the dust opacity and stellar emissivity for spiral galaxies from the UV to submm wavelengths using parameterized analytic functions. The model considers the absorption and anisotropic scattering of stellar photons with dust grains of various size and chemical composition. A dust composition consisting of silicates, graphites and PAH molecules was used, with a grain size distribution from Weingartner & Draine (2001), and optical constants from Draine & Lee (1984) and Draine & Li (2007). As in our previous work on modelling spiral galaxies, we consider a Milky Way-type dust, since this has proven to provide consistent results for all modelled galaxies. The parameters describing the grain size distribution and grain mixture are taken from Table 1, line 7 of Weingartner & Draine (2001). The RT model of M51 also incorporates explicit calculations for the stochastic heating of the dust grains.

The RT codes used for the modelling of M51 were the PT11 code and the DART-RAY code (Natale et al. 2014, 2015, 2017). The PT11 code is a modified version of the code of Kylafis & Bahcall (1987). Both the PT11 and DART-RAY codes use ray-tracing algorithms.

Following Thirlwall et al. (2020), we adapt the model of PT11 to allow for additional morphological components. In the case of M51, it was found that, in addition to the bulge, an inner and a main disc component were present out to R_l. Each morphological component comprises the generic stellar and dust components from PT11: the stellar disc, the thin stellar disc, the dust disc, and thin dust disc. The stellar disc is made up of all the stars that are old enough to have had time to migrate and form a thicker configuration (a few hundred parsec scale height), while the thin stellar disc is made up by young stars, that are still spatially closely related to the molecular layer from which they formed (50–100 parsec scale height). A clumpy component associated with the star-forming clouds is also incorporated in the model. All these are described in detail in PT11. Following Thirlwall et al. (2020), we also found the need to alter the analytic functions describing the exponential discs, to include an inner truncation radius. Thus, we use equations (1) and (2) from Thirlwall et al. (2020) to describe the stellar volume emissivity and the dust-density distribution for each disc component j:

with

where R and z are the radial and vertical coordinates, h_j and z_j are the scale length and scale height, respectively, A_{0, j} is a constant scaling the amplitude, χ_j describes a linear slope of the radial distributions from the inner radius, and R_in, to the centre of the galaxy. R_tin and R_t are the inner and outer truncation radii, respectively. Equations (1) and (2) are wavelength dependent. The spatial integration of these formulae (equations B1–B3 from Thirlwall et al. 2020) between R_tin and R_t provides the intrinsic stellar luminosity L_ν if j represents a stellar disc component, or the dust mass M_d if j represents a dust disc component.

The emissivity of the bulge is described by the Sérsic distribution

with

η^bulge(λ, 0, 0) is the stellar emissivity at the centre of the bulge, R_e is the effective radius of the bulge, a and b are the semimajor and semiminor axes of the bulge, respectively, and b_s is a constant which depends on the value of the Sérsic index n_s (see equation 8 from Natale et al. 2021). Following Thirlwall et al. (2020), the parameters of the different components are as follows:

The stellar disc: The stellar disc is made up of the old stellar population, defined as a population of stars that had time to migrate from their molecular layer into a thicker configuration (a few hundred parsec scale height). Following from Popescu et al. (2000) and PT11, we make the approximation that the stellar disc only emits in the optical/NIR (no UV counterpart). It is described by the four geometrical parameters |$h_{\rm {s}}^\text{disc}$|⁠, |$z_{\rm {s}}^\text{disc}$|⁠, |$R_{\rm {in,s}}^\text{disc}$|⁠, and |$\chi _{\rm {s}}^\text{disc}$|⁠, and one amplitude parameter L^disc. The free parameters are |$h_{\rm {s}}^\text{disc}$|⁠, |$\chi _{\rm {s}}^\text{disc}$|⁠, and L^disc(λ). |$z_{\rm {s}}^\text{disc}$| is fixed from the model, and |$R_{\rm {in,s}}^\text{disc}$| is fixed from the data. We allow L^disc to be a free parameter for each wavelength corresponding to u, g, r, i, J, K, I1, I2, and I3. We also allow h_s to be wavelength dependent, hence we consider the scale lengths for the g, r, i, J, K, I1, I2, and I3 bands to be free parameters. The scale length in the u band is taken to be the same as in the g band. The model of M51 includes an inner and a main stellar disc.

The thin stellar disc: The thin stellar disc is made up of the young stellar population of stars that did not have time to migrate from their molecular layer (50–100 pc). This disc emits predominantly in the UV, although we assume a smooth decrease of its emission in the optical/NIR. It is described by the four geometrical parameters |$h_{\rm {s}}^\text{tdisc}$|⁠, |$z_{\rm {s}}^\text{tdisc}$|⁠, |$R_{\rm {in},s}^\text{tdisc}$|⁠, and |$\chi _{\rm {s}}^\text{tdisc}$|⁠, and one amplitude parameter L^tdisc. The free parameters are |$h_{\rm {s}}^\text{tdisc}$|⁠, |$\chi _{\rm {s}}^\text{tdisc}$|⁠, and L^tdisc. |$z_{\rm {s}}^\text{tdisc}$| is fixed from the model and |$R_{\rm {in},s}^\text{tdisc}$| is fixed from the data. The scale length of this disc is assumed to be wavelength independent. The wavelength dependence of L^tdisc is fixed from the model. The model of M51 includes an inner and a main thin stellar disc.

Keeping in line with previous modelling, we express the spectral integrated luminosity of the young stellar disc L^tdisc in terms of a star formation rate SFR, using equations (16)–(18) from PT11.

The dust disc: It is described by the geometrical parameters |$h^{\rm {disc}}_{\rm {d}}$|⁠, |$z_{\rm {d}}^{\rm {disc}}$|⁠, |$R_{\rm {in,d}}^{\rm {disc}}$|⁠, and |$\chi _{\rm {d}}^{\rm {disc}}$| and the amplitude parameter, the face-on optical depth of the B band at the inner radius |$\tau _{\rm {B}}^{\rm {f, disc}}\left(R_{\rm {in,d}}^{\rm {disc}}\right)$|⁠. The free parameters are |$h^{\rm {disc}}_{\rm {d}}$|⁠, |$\chi _{\rm {d}}^{\rm {disc}}$| and |$\tau _{\rm {B}}^{\rm {f,disc}}\left(R_{\rm {in,d}}^{\rm {disc}}\right)$|⁠. |$z_{\rm {d}}^{\rm {disc}}$| is fixed from the model and |$R_{\rm {in,d}}^{\rm {disc}}$| is fixed from the data. Our model of M51 contains an inner, and main dust disc.

The thin dust disc: The thin dust disc is described by the geometrical parameters |$h^{\rm {tdisc}}_{\rm {d}}$|⁠, |$z_{\rm {d}}^{\rm {tdisc}}$|⁠, |$R_{\rm {in,d}}^{\rm {tdisc}}$|⁠, and |$\chi _{\rm {d}}^{\rm {tdisc}}$| and the amplitude parameter, the face-on optical depth of the B band at the inner radius |$\tau _{\rm {B}}^{\rm {f, tdisc}}\left(R_{\rm {in,d}}^{\rm {tdisc}}\right)$|⁠. Following PT11, we fix the geometric parameters of the thin dust disc to that of the thin stellar disc. The model of M51 contains an inner, and main thin dust disc.

The clumpy component: The clumpy component represents the dense dusty clouds associated with star-forming regions that are illuminated by the strong, UV-dominated radiation fields produced by the photons emitted from the young stars inside the clouds. The dust in these regions reaches thermal equilibrium with the intense radiation fields and emits dominantly in the 24- to 70-μm region, often surpassing the emission due to the diffuse component at these wavelengths. Following PT11, we define a ‘clumpiness factor’ F, which is the amplitude parameter of the clumpy component, and represents the fraction of the total luminosity of massive stars that is locally absorbed by the dust in these star-forming clouds. The diffuse component is then illuminated by the escape fraction 1 − F. The F factor has a wavelength dependence that is not determined by the optical properties of the dust grains, but by the evolutionary stage of the birth clouds. This is because in our model the dust cloud is very opaque, thus completely blocking any radiation from inside the cloud. The only way photons escape is through holes in the cloud, produced by the fragmentation in the cloud. The wavelength dependence of the escaping radiation thus arises because lower mass and less blue stars spend a higher proportion of their lifetime radiating when they are further away from their birth-clouds, and because of the progressive fragmentation of the clouds. The formalism is described in Tuffs et al. (2004). The dust emission SED of the birth clouds is taken to follow the prescription from PT11.

4 FITTING THE MODEL

To directly compare the model to the data, we produce azimuthally averaged surface brightness profiles of both. The observed surface-brightness profiles were obtained by carrying out CoG photometry on the maps listed in Table 2, as described in Section 2. For the model profiles, we carried out an identical procedure. Examples of these profiles can be seen in Figs 2 and 3. Examples of model images are given in the right-hand panels of Fig. 1.

From an initial inspection of the surface brightness profiles of the observational data, it was clear that a single disc described by equations (1) and (2) would provide an inadequate fit. Following Thirlwall et al. (2020), we segment the galaxy into distinct morphological components to fully capture the detailed morphology of M51. As mentioned in Section 3, we found that, within R_l, in addition to the bulge, two distinct morphological components were needed to model this galaxy: an inner and a main disc. We allow for each morphological component to be made up of two stellar and two dust components.

We followed the same fitting procedure as in Thirlwall et al. (2020). In brief, this involves an intelligent searching algorithm of the parameter space, that takes into account the fact that different parameters do not equally affect the emission at all wavelengths. This allows one to identify pairs of parameters that predominantly shape emission at key wavelengths and use this to fit one pair (usually made of a geometric parameter and an amplitude parameter) of parameters at a time.

Thus, we begin modelling with some initial guess parameters. These parameters are based on notable features of the surface brightness profiles, as well as informed by previous modelling (Thirlwall et al. 2020). We start by fitting the NUV profile. This is equivalent to constraining the spatial distribution of the young stellar population residing in the thin stellar disc. Although this stellar disc is highly attenuated in the UV, we can still derive an estimate of the corresponding geometrical and amplitude parameters for an initial guess of dust opacity. We thus derive an exponential scale length |$h^{\rm {tdisc}}_{\rm {s}}$| and a luminosity density L^tdisc in the NUV filter, for the inner and main disc, assuming a known wavelength dependence of the luminosity density at the other sampled wavelengths. There is also the free parameter |$\chi _{\rm s}^{\rm tdisc}$|⁠, which is fitted for the inner disc in order to account for the fall of emissivity towards the galaxian centre. For the main disc, |$\chi _{\rm s}^{\rm tdisc}$| is only used to provide a smooth overlap between the inner and the main disc, resembling the observed averaged profiles, instead of simply using delta functions. |$\chi _{\rm s}^{\rm m-tdisc}$| is therefore not a critical parameter in the fit. Inner radii R_in and inner and outer truncation radii R_tin and R_t are fixed from data and are not subject to iterations. Once we obtain a good fit in the NUV (for the trial dust opacity), we then attempt to constrain the dust distribution. To do this, we re-run the RT calculation with the parameters derived from fitting the NUV profile and for an initial guess for the optical emission, and predict the dust emission at all IR wavelengths, including the 500 μm. We then constrain the spatial distribution of the dust by comparing the predicted 500-μm model profiles with the corresponding observed ones, since at this wavelength the dust emission primarily traces dust column density rather than heating sources. At this step, small changes in the dust parameters may be needed for a best fit to the 500-μm data. This results in a fit for the scale length of the dust disc |$h^{\rm {disc}}_{\rm {d}}$| and a dust opacity (plus the parameter |$\chi _{\rm d}^{\rm tdisc}$|⁠). A few iterations are needed between NUV and the 500 μm, until a good fit at both NUV and the 500 μm is achieved. Towards the end of the iteration the FUV luminosity density is also slightly adjusted until the best fit to the FUV profile is also obtained. Once the parameters related to the thin stellar disc and dust disc are constrained, we run again the RT code to derive the distribution of the stellar disc and bulge harbouring the old stellar population, by fitting the optical and NIR profiles at each wavelength for which observational data is available. Since in a face-on galaxy the optical and NIR emission is fairly optically thin to dust, the final fits at these wavelengths do not alter the solution for the dust disc and the thin stellar disc (in the UV) obtained in previous iterations. A clumpy component is finally fitted to account for any remaining 24–70-μm-observed emission not accounted for by the diffuse dust emission. Overall, within a small number of iterations a solution is derived for the whole spectral range. Notably, the 70–350-μm emission is not directly used in the fit, but instead is predicted. To conclude, the fit mainly consists in fitting only a pair of parameters (a scale length and an amplitude) at any given time, for each morphological component at a given wavelength. More details on the fitting procedure can be found in Thirlwall et al. (2020).

The main (free) geometric parameters constrained by the procedure outlined above are listed in Table 3. The free parameters related to stellar luminosity densities are listed in Table A1. The free parameters related to dust opacity are listed in Table A2. Other parameters of the model are listed in Table A3.

The results of the fits are presented in Figs 2 and 3, showing an overall reasonable agreement with the data. To help the comparison, the residuals R between the surface brightness profiles of the observations and of the model are also plotted on the bottom panels of Figs 2and3. The residuals are defined as

In addition, the goodness of fit is also considered with the reduced chi-squared |$\text{chi}^{2}_{r}$| statistics,⁴ which becomes more relevant when finalizing the model solution. We calculated |$\text{chi}^{2}_{r}$| for each wavelength by the following formulae:

where N is the total number of annuli n, m is the number of free parameters fitted at wavelength λ, O_n is the azimuthally averaged flux through annuli n, M_n is the model flux through annuli n, and |$\epsilon ^{2}_{\text{SB,}n}$| is the total error of the azimuthally averaged surface brightness profile corresponding to annuli n. The total number of annuli is 70 while m varies between 2 in the FUV, to 5 (in the NUV and 500 μm), to 8 (in the optical and NIR) for a fitted profile out to 7 kpc radius. At most other wavelengths that are predicted, we can assume zero free parameters. Because N is much larger than m we neglect m in equation (7). Values of the |$\text{chi}^{2}_{\lambda }$|⁠, for selected wavelengths, are listed in Table 4. The total reduced chi-squared statistic over all wavelengths is

The |$\text{chi}^{2}_{r}$| values from Table 4 confirm the trends seen in Figs 2 and 3, proving that the axi-symmetric model provides a reasonable fit to the data.

5 RESULTS

5.1 The surface brightness distributions

Taking the azimuthally averaged surface brightness of the observations produces relatively smooth exponential profiles, as seen in Figs 2 and 3. Despite the slight imperfections in the smoothness, the average profiles are still suitable for description with analytic axi-symmetric functions. The most significant deviations are in the UV and 24-|$\, \mu$|m bands, where strong asymmetries and clumpiness are present. The bumps in these profiles generally coincide with the spiral arms as well as star-forming regions. For example, a significant bump in emission at 1.3 kpc, which can be seen most notably in the FUV profile, is due to bright star-forming regions. There is also a bump in emission between 5 and 6.5 kpc, which coincides with an overlap in the spiral arms. Predictably, the 24-μm profiles exhibit the least smooth curves, as at this wavelength the emission due to the diffuse dust and the dust heated by localized SF regions is roughly equal, hence strong asymmetries arise from these localized sources of emission. This is also seen in the IRAC bands, where there is significant PAH emission.

Overall, the model has reasonable agreement with the data with residuals generally within 20  per cent.

The scale length of the observed profiles in the optical regime generally decreases with increasing wavelength for each morphological component. For wavelengths |$\ge 3.6\, \mu$|m, this trend reverses for the main disc. This is similar with trends found from previous modelling (M33; Thirlwall et al. 2020).

The observed profiles in the IRAC bands are modelled with emission from the old stellar populations as well as dust emission, with a negligible contribution from the young stellar population. At 3.6 and |$4.5\, \mu$|m the emission is dominated by the old stellar population, but at |$5.8\, \mu$|m the dust emission is almost equal in amplitude to that of the old stellar population. At the |$8.0{-}\mu$|m band, the emission is dominated by the dust, primarily due to the PAH emission.

5.2 The global SED of M51

The global SED of M51 was produced by spatially integrating the emission at each wavelength out to R_l. This was done both for observations and for the model and the results are plotted in Fig. 4. The model SED closely resembles the observed SED, with a maximum residual of |$|R|=26.2~{{\ \rm per\ cent}}$| and an average residual of |$\langle |R|\rangle =4.33~{{\ \rm per\ cent}}$|⁠. The model performs worst at |$160\, \mu$|m, where the emission is underestimated by our model.

We also show the contributions from the different morphological components to the global SED in Fig. 5. It is important to note that the dust emission for a given morphological component is not entirely powered by the stellar emission of the same component, as for example, some fraction of the stellar photons emitted from the inner disc will contribute to heating the dust in the main disc. None the less, we still show the separate contributions of the inner and main discs to the dust emission, within the accepted understanding that there is not a one-to-one spatial correspondence between the stellar and the dust emission of these components. As expected due to its extent, the largest contribution to the global SED is from the main disc responsible for 67.4 per cent of the total stellar emission, and 63.7 per cent of the dust emission. The inner disc contributes 30.7 per cent to the total stellar emission and 36.3 per cent to total dust emission. The bulge contributes 1.89 per cent to the total stellar emission. Here, we note that the separation between the inner disc and the bulge is rather degenerate. We found no strong constraints for this separation. We therefore consider the bulge versus inner disc luminosity to be subject to large uncertainties.

The dust emission SED is dominated by the diffuse component. The H ii component has a smaller contribution, of 12.2  per cent. We find a clumpiness factor of F = 0.125. The global dust emission SED peaks at |$130\, \mu$|m, which corresponds to a dust temperature of 22 K. The hottest dust resides within the inner disc, where the SED peaks at |$102\, \mu$|m, corresponding to a dust temperature of 28 K. The main disc has the greatest contribution to the global SED, therefore, as expected, it peaks at a similar wavelength, at around |$133\, \mu$|m, corresponding to a dust temperature of 22 K.

We find that 38.3 ± 4 per cent of the stellar light is absorbed and re-radiated by the dust. This is a typical value for late-type spiral galaxies (Popescu & Tuffs 2002). M33 was found to have a similar value of 35 ± 3 per cent (Thirlwall et al. 2020). We also find that the young stellar population contributes 42.8 per cent to the total intrinsic stellar light, with 55.4 per cent being due to the old stellar population in the disc and the remaining 1.89 per cent due to the bulge.

5.3 Star formation rates in M51

We obtain an |$\text{SFR}=4.11^{+0.41}_{-0.39}\,$| M_⊙yr⁻¹ which is within the broad range of literature values. De Looze et al. (2014) obtained a |$\text{SFR}\sim 3.13 \,$| M_⊙yr⁻¹(⁵) and Nersesian et al. (2020b) found an |$\text{SFR}=4.09\pm 0.35\,$| M_⊙yr⁻¹(⁶). Both these studies use SKIRT, a Monte Carlo RT model. Calzetti et al. (2005) obtain the estimates |$\text{SFR}(24\, \mu \text{m})\sim 3.7\,$| M_⊙yr⁻¹(⁷) and |$\text{SFR(FUV)}\sim 4.7\,$| M_⊙yr⁻¹(⁸).

About two-thirds of the star formation takes place in the main disc with |$\text{SFR}^{\text{m}}=2.73^{+0.28}_{-0.26}\,$|  M_⊙yr⁻¹. This is expected due to both the size and bolometric output of the main disc. For the inner disc, we find |$\text{SFR}^{\text{i}}=1.38^{+0.15}_{-0.14}\,$|  M_⊙yr⁻¹. We also calculate the surface density of star formation rate, Σ_SFR, for each morphological component, by considering the physical extent of each disc, which is shown in Tables B1 and B2. The global surface density in SFR is |$\Sigma _{\text{SFR}}=2.83^{+0.28}_{-0.27}\times 10^{-2}\,$| M_⊙yr⁻¹kpc⁻².

5.4 Dust optical depth and dust mass

We find that the face-on optical depth in the B band has a maximum value of |$\tau _{\text{B}}^{\text{f}}\left(R_{\text{in}}\right)=5.33\pm 0.33$| at the inner radius of the inner disc. The main disc has a face-on optical depth at the inner radius of 2.62 ± 0.08. These values are listed in Table A2. We find that M51 has a global dust mass (out to R_l) of M|$_{\text{d}}=3.74^{+0.34}_{-0.30}\times 10^{7}\, \text{M}_{\odot }$|⁠. The dust masses of the individual morphological components are listed in Table B3.

Using a gas mass of M|$_{\text{G}}=1.1\pm 0.03\times 10^{10}\, \text{M}_{\odot }$| (Mentuch Cooper et al. 2012), we obtain a gas-to-dust ratio of GDR = 294 ± 4.

5.5 Morphological components of M51

In modelling M51, we found three distinct morphological components out to R_l, each requiring different geometrical parameters for the distribution of stars and dust. Although we present the model for the entire galaxy, it is also interesting to study the properties of these individual components. The geometrical parameters of the three morphological components, as constrained by our model are listed in Table 3.

5.5.1 The inner disc

The inner disc extends from 0  to 1.9 kpc and exists as a ring-like structure for the young stars and dust, with an inner radius of 700 and 800 pc, respectively. The distribution of the old stars, however, exponentially increases towards the centre. The distribution of the dust is notably flat with a scale length of 5 kpc, compared to a scale length of 600 pc for the young stellar distribution. The inner disc has an |$\text{SFR}^{\text{i}}=1.38^{+0.15}_{-0.14}\,$| M_⊙yr⁻¹, contributing around 30 per cent to the global SFR, and has the greatest surface density of star formation rate |$\Sigma _{\text{SFR}}^{\text{i}}=1.21^{+0.28}_{-0.11}\times 10^{-1}\,$| M_⊙yr⁻¹kpc⁻². The dust mass within the inner disc is M|$_{\text{d}}^{\text{i}}=6.39^{+0.37}_{-0.21}\times 10^{6}\, \text{M}_{\odot }$|⁠. The temperature of the dust contained within the inner disc is the hottest in the galaxy, at around 28 K.

5.5.2 The bulge

The bulge was found to have an effective radius of 350 pc and was modelled with a Sérsic index of 4. As mentioned in Section 5.2, the bulge parameters are subject to uncertainty, in the sense that there are different configurations in the bulge/inner disc decomposition that could result in a good fit to the data. We based our choice on covering a small dip in the profile, that would otherwise appear in other configurations, and as such we believe that our choice is the most reliable result. None the less, we caution the reader of the subjectivity of the approach. The error estimates on the bulge/inner disc parameters do not take into account the uncertainty due to the other possible combinations.

5.5.3 The main disc

The main disc resides in the region between 1.8 and 6.8 kpc. The interior of the main disc slightly overlaps with the outermost limits of the inner disc. Both the stellar and dust distribution are relatively flat in the main disc, with a scale length of 4.3 kpc for the young stellar disc and 6 kpc for the dust disc. The main disc contributes most significantly to the global SFR, with a SFR|$^{\text{m}}=2.73^{+0.28}_{-0.26}\,$| M_⊙yr⁻¹. The main disc has a surface density SFR of |$\Sigma _{\text{SFR}}^{\text{m}}=2.00^{+0.18}_{-0.17}\times 10^{-2}\,$| M_⊙yr⁻¹kpc⁻². The main disc has a dust mass of M|$_{\text{d}}^{\text{m}}=3.10^{+0.26}_{-0.18}\times 10^{7}\, \text{M}_{\odot }$|⁠, with a dust temperature of 21.8 K.

5.6 The contribution of young and old stellar populations in heating the dust

The fractional contribution of the different stellar populations in heating the dust is an important quantity, since it allows a quantitative understanding of the physical mechanisms related to the interaction between stellar light and interstellar dust. In addition, this is key to the understanding of key correlations, like, for example, the tight and universal correlation between FIR and radio continuum emission of spiral galaxies (de Jong et al. 1985; Helou, Soifer & Rowan-Robinson 1985; Wunderlich, Klein & Wielebinski 1987). The standard interpretation of the so-called FIR-radio correlation is in terms of young and massive star formation activity (Condon 1992). This picture assumes that the dust emission is mainly powered by young stars which are also responsible for the radio emission: the ionizing radiation from the young stars powers the thermal radio emission, and the remnants of the supernova explosions which occur at the end of their lives accelerate the cosmic-ray electrons. None the less both the effect of dust opacity and the contribution of the old stellar populations play a role in shaping the slope of the correlation (Pierini et al. 2003), and as such, knowledge of the fraction of dust heating powered by the young and old stellar populations is a pre-requisite for the analysis of the FIR/radio correlation. Furthermore, knowledge of the different stellar populations in heating the dust is equivalent to understanding the dust attenuation of the different stellar populations, which, in turn, allows population synthesis models to be used for deriving star formation histories in galaxies.

For the global emission of M51, we found that the fractional contribution of the young stellar population in powering the dust emission, |$F^{\rm dust}_{\rm young}$|⁠, is 69 per cent, the rest being attributed to the old stellar populations in the disc and bulge. Previous modelling of M51 with RT methods derived |$F^{\rm dust}_{\rm young}=0.63$| (De Looze et al. 2014) and |$F^{\rm dust}_{\rm young}=0.72$| (Nersesian et al. 2020b). The latter value is more in range with our determination.

Our value of |$F^{\rm dust}_{\rm young}=0.69$| is similar to that derived with the same models for other nearby spirals: |$F^{\rm dust}_{\rm young}=0.71$| for the Milky Way (Natale et al. 2021) and NGC 628 (Rushton et al. 2022), but slightly lower than |$F^{\rm dust}_{\rm young}=0.8$| for M33 (Thirlwall et al. 2020).

Following the formalism from Natale et al. (2015) and Popescu et al. (2000), we use the RT calculations to not only derive |$F^{\rm dust}_{\rm young}$| and |$F^{\rm dust}_{\rm old}$| for the global emission, but also to understand the spatial variation of this quantity. In Fig. 6, we plot |$F^{\rm dust}_{\rm young}$| and |$F^{\rm dust}_{\rm old}$| as a function of radial distance. Within the inner 0.5-kpc dust emission is dominated by the old stellar populations, mainly from the bulge. Beyond 0.5 kpc, it is the young stellar population that contributes a significant fraction of the dust heating. At around 0.7–0.8 kpc (the inner radius of the inner disc), there is a peak in |$F^{\rm dust}_{\rm young}$|⁠, reaching above 70 per cent. Within the extent of the main disc, between 2.7 and 7 kpc, |$F^{\rm dust}_{\rm young}$| is slowly increasing with radial distance, from 60 to 85 per cent. At around 2 kpc, corresponding to the transition region between the inner and the main disc, |$F^{\rm dust}_{\rm young}\approx F^{\rm dust}_{\rm old}$|⁠.

The overall trends are similar to those found for the Milky Way (Natale et al. 2021) and for NGC 628 (Rushton et al. 2022), in the sense that dust heating in the inner regions has a significant contribution from old stellar populations, while the rest of the dust in the disc is mainly heated by young stellar populations, with |$F^{\rm dust}_{\rm young}$| rather constant or slowly increasing towards the outer disc. M51 has larger radial variations in F^dust than the Milky Way or NGC 628.

5.7 The radial variation of Σ_SFR, |$\Sigma _{\text{M}_{\star }}$|⁠, and specific star formation rate

The global SFR of M51 was found to be |$4.11\, \text{M}_{\odot }\, {\rm yr}^{-1}$|⁠. The rate of star formation is not constant throughout the extent of M51, but has a strong radial dependence. To illustrate this, we plot in the top panel of Fig. 7 the surface density of SFR, Σ_SFR, versus radial distance. One can see that there is a peak in Σ_SFR at the inner radius of the inner disc, followed by a steep decrease out to around 2 kpc. Beyond 3 kpc, within the extent of the main disc, the Σ_SFR follows a shallow decrease. Overall, Σ_SFR shows a marked difference between the inner and the main disc. By contrast, the surface density of stellar mass, |$\Sigma _{\text{M}_{\star }}$|⁠, has a smoother radial variation, (middle panel, Fig. 7), showing a shallow decline with increasing radial distance. The stellar mass was calculated using the 3.6- and 4.5-μm flux densities and the calibration from Eskew, Zaritsky & Meidt (2012).

The specific star formation rate (sSFR) is rather constant throughout the disc (bottom panel, Fig. 7), except for the very central region. These trends are similar to the trends found for M33, NGC 628, and the Milky Way.

5.8 The attenuation curve of M51

The attenuation curve of a galaxy depends not only on the optical constants of the dust grains and their grain size distribution (incorporated in the dust extinction), but also on the relative distribution of the stars and dust. Neither the extinction nor the geometry is usually known, making predictions for attenuation rather difficult. It is, however, the strength of radiative-transfer modelling to accurately derive geometries of stars and dust, within the framework of a given dust model.

In Fig. 8, we plot the predicted attenuation curve of M51, as derived from our RT model. For comparison, we also plot the extinction curve of the MW from Fitzpatrick (1999). We fit the attenuation curve of our model with a third-order polynomial plus a Drude profile, which is the functional form presented in Salim, Boquien & Lee (2018). Using this form, we found that the best fit is given by

The attenuation curve of M51 is slightly steeper than the MW extinction curve, becoming increasingly disparate at shorter wavelengths. M51 also exhibits a much greater 2200 -Å  bump compared to the MW extinction. We find that the difference in width of the bumps to be negligible.

6 DISCUSSION

6.1 Comparison with the non-axi-symmetric model of Nersesian et al. (2020b)

As mentioned before, M51 was also modelled with non-axi-symmetric RT codes by De Looze et al. (2014) using the SKIRT code (Baes et al. 2011; Camps & Baes 2020), and more recently by Nersesian et al. (2020b), using a more up-to-date version of the same formalism and code. Here we compare our model to that of Nersesian et al. (2020b). Publicly available model images of Nersesian et al. (2020b) were extracted from DustPedia⁹ and processed in the same way as our model images, to produce azimuthally averaged radial profiles. We also checked that the observed fluxes used in Nersesian et al. (2020b) were consistent with the values used in our work, and, when small differences existed (at per cent level), we adjusted for this. In Fig. 9, we show examples of this comparison for a few selected wavelengths. Overall our model performs better. This is to be expected, as the model of Nersesian et al. (2020b) does not fit the geometry of the system, but only the spatially integrated SED.

6.2 An outer disc

All the results presented so far were confined to the main extent of the galaxy, out to the radial distance R_l = 7 kpc. Note, however, that R_l is not the actual truncation radius of M51. There is very faint and asymmetric emission beyond R_l, also containing the bridge between M51 and its companion M51b. This faint emission does not obey axi-symmetry, therefore it is not a prime target for an axi-symmetric model. None the less, the emission is so faint that it is difficult to trace and model, and one effective way to reveal it is through azimuthally averaging. By averaging, a strong signal becomes apparent, in a form of an exponential decay. We model this outer emission with another morphological component that we call ‘outer disc’ . While we acknowledge that an axi-symmetric model is a poor approximation for this emission, we nevertheless attempt to at least give some characterization of the global properties of this emission.

The model for the outer disc morphology contains the same geometrical components as for the main disc: a stellar disc, a thin stellar disc, a dust disc, and a thin dust disc. The analytic functions we fit are the same, and the whole optimization follows the same procedure as for the main body of the galaxy.

Examples of fits to the outer disc are given in Fig. 10. We find the outer disc to extend from 6.7  to 20 kpc, with a scale length of 3.6 and 1.45 kpc for the dust disc and young stellar disc, respectively. The main geometrical parameters of the outer disc are listed in Table 5.

The outer disc is the morphological component with the smallest SFR, with an SFR|$^{\text{o}}=0.57\pm 0.06\,$| M_⊙yr⁻¹. The surface density of SFR is the least dense of the components, with |$\Sigma _{\text{SFR}}^{\text{o}}=5.11\pm 0.56\times 10^{-4}\,$|M_⊙yr⁻¹kpc⁻².

6.3 A model with a spatial variation in dust properties

Inspection of the fitted profiles from Figs 2 and 3 shows that the axi-symmetric model provides a reasonable fit to the data, within R_l.

However, when deriving the spatially integrated SED (see Fig. 11, top left-hand panel), we noticed that our predicted intrinsic (de-attenuated) stellar SED shows a spike in the NUV, exactly at the position of the well-known ‘2200-Å bump’ of the extinction curve. It can be seen that the inner disc contributes most significantly to this spike, while the main disc has negligible contribution, if any. The intrinsic stellar SED should be smooth, with no feature at this wavelength. Because our model is in good agreement with the data at all wavelengths, as seen in Figs 2 and 3, the spike in the intrinsic SED is indicative that the inner disc is being too heavily attenuated in the NUV. This could be due to a spatial variation in the dust properties, in the sense that the dust properties may be different in the inner disc than assumed in our model. As mentioned before, our model for M51 uses the Milky Way-type dust from Weingartner & Draine (2001), described on line 7 in table 1 from Weingartner & Draine (2001). The extinction curve of the Milky Way dust model along with the contributions from the different grain types is shown in the left-hand panel of Fig. C1, together with its prominent ‘2200 Å  bump’. This feature is primarily due to the PAH molecules, which the model considers to be the carbonaceous grains with sizes |$a\le 0.01\, \mu$|m. The overattenuation of the NUV photons in the inner disc could suggest that the extinction curve corresponding to the dust model that we use has a too strong bump, meaning an overabundance of PAH in the inner disc, than is present in M51. We note that it is extremely hard to constrain the shape of the UV extinction curve of galaxies, especially in the inner region where the attenuation is high and the geometry is complex, with just the GALEX, FUV, and NUV bands. Several authors have shown that these bands alone are insufficient to disentangle the shape of the continuum and the strength of the 2200-Å bump (Hoversten et al. 2011; Hutton, Ferreras & Yershov 2015; Hagen et al. 2017; Decleir et al. 2019).

None the less, within the framework of a model with varying dust properties, our fits at |$8\, \mu$|m (see Fig. 11, bottom right-hand panel) provide additional evidence that we may need a different type of dust in the inner disc. The model profile is clearly overpredicted in the inner region. Since the |$8\, \mu$|m is a strong tracer of PAH abundance, in addition to tracing the radiation fields that are self-consistently calculated in our RT model, the overprediction could again be indicative of an overabundance of PAH in our model of the inner disc.

As an alternative, we sought a dust model with a reduced PAH abundance and chose one of the LMC dust models of Weingartner & Draine (2001), namely the LMC average (line 1 from table 3 of Weingartner & Draine 2001). This LMC-type dust features the same optical constants but different grain size distributions, particularly less PAH grains. Thus, the extinction curve exhibits a less pronounced 2200-Å bump, as seen in the right-hand panel of Fig. C1. As the attenuation in the main and outer discs produced smooth intrinsic SEDs, we chose to modify the dust properties only in the inner disc (using the LMC-type dust). We thus fitted a new model of M51 with variable dust grain properties. We refer to this alternative as the ‘hybrid model’. The new fits produce the intended effect, reducing the NUV peak in the de-attenuated stellar SED of the inner disc, as seen in Fig. 11. The new fits to the |$8\, \mu$|m profiles (see Fig. 11, bottom right-hand panel) also provide the right level of emission at this wavelength. We conclude that using the LMC-type dust for the inner disc could provide an alternative to the standard model.

The finding that PAHs could be less abundant in the inner disc than in the rest of the galaxy raises the question of what physical processes could cause this variation. Theoretical and observational studies suggest that strong UV radiation fields could destroy PAH molecules (e.g. Boulanger et al. 1988; Helou, Ryter & Soifer 1991; Contursi et al. 2000; Siebenmorgen & Krügel 2010). Taking into account that the inner region of M51 has the highest surface density of SFR, and therefore intense UV radiation fields, it is plausible that PAHs are more readily destroyed in this region, thus explaining the reduced abundance inferred by the hybrid model. Several studies have suggested that the fraction of PAHs becomes low in H ii regions (e.g. Pety et al. 2005; Lebouteiller et al. 2007; Thilker et al. 2007; Chastenet et al. 2023) and is a strong function of environment, with higher fractions in less harsh environments due to star formation (Chastenet et al. 2023). It is beyond the scope of our paper to provide a quantitative model for PAH destruction in different ISM environments, but we note the possibility of PAH abundance variation within galaxies.

The fits to the multiwavelength data using the hybrid model do not provide any big differences from the standard model, with the exception of the 8-μm profile, which is largely improved (see examples in the Appendix, in Figs C2 and C3). The energy balance between direct and re-radiated stellar light seems to be the same as in the standard model, with no major improvements in the fit. The values of the chi² are also similar. The individual values at selected wavelengths are listed in Table B4. The values of the derived intrinsic parameters are different, but the quality of the fit is the same. This is perhaps not surprising, since the optical properties of the grains are not changed. There are some other dust models available in the literature, in particular ‘The Heterogeneous dust Evolution Model for Interstellar Solids’ (THEMIS) (Jones et al. 2013; Köhler, Jones & Ysard 2014; Ysard et al. 2015), which consider different carbonaceous grains, with aliphatic rather than aromatic molecules. The THEMIS model has enhanced submm efficiencies for grains, by a factor of about 2.5 with respect to the Draine model. The THEMIS model would produce a different energy balance and therefore different fits.

To conclude, the hybrid model with a reduced PAH abundance in the inner disc, together with a modified grain size distribution in the central region, provides a consistent solution to the panchromatic data of M51, and is based on observational evidence coming from both dust extinction and dust emission. None the less, we cannot prove that these observational constraints are uniquely met by our hybrid model.

7 SUMMARY

Using the generic RT code and formalism of PT11, we derived an axisymmetric model for M51, fitted to multiwavelength imaging data ranging from FUV to submm. We find that, despite the interaction with the companion galaxy M51b, M51 preserves a regular spiral structure within a radial distance of R_l = 7 kpc, making it thus suitable for the assumption of axi-symmetry within R_l. The model fits to the azimuthally averaged surface brightness profiles were found to be in reasonable agreement with the observations.

We find three distinct morphological components out to R_l: a bulge, an inner disc and a main disc. The bulge has an effective radius R_eff = 350 pc and a Sérsic index n = 4. The inner disc is prominent in the UV, with a scale length of the thin stellar disc, |$h_{\rm s}^{\rm i-disc}$|⁠, of 600 pc. The distribution of dust in the inner disc is very flat, with a scale length of the dust disc, |$h_{\rm d}^{\rm disc}$|⁠, of 5  kpc. The main disc resides in the region between ∼1.8 and ∼6.8 kpc. Both the stellar and the dust distributions of this morphological component are relatively flat, with a scale length for the thin stellar disc, |$h_{\rm s}^{\rm t-disc}$| of 4.3 kpc and |$h_{\rm d}^{\rm disc}=6$| kpc for the dust disc.

The main intrinsic properties of M51, as derived from the model are as follows:

• The global star-formation rate is |$\text{SFR}=4.11^{+0.40}_{-0.38}\,$|  M_⊙yr⁻¹. The inner disc has |$\text{SFR}=1.38^{+0.15}_{-0.14}\,$|  M_⊙yr⁻¹and the main disc has |$\text{SFR}=2.73^{+0.28}_{-0.26}\,$|  M_⊙yr⁻¹.

• The global surface density in SFR is |$\Sigma _{\text{SFR}}=2.83^{+0.28}_{-0.27}\times 10^{-2}\,$|  M_⊙yr⁻¹kpc⁻².

• The sSFR is rather constant over much of the extent of M51, except for the very centre.

• The young stars account for 69 per cent of the dust heating. In the inner regions, dust heating has a significant contribution from old stellar populations, while the rest of the dust in the disc is mainly heated by young stellar populations, with their contribution rather constant or slowly increasing towards the outer disc.

We discuss the properties of the faint and non-axi-symmetric emission extending beyond R_l, also containing the bridge between M51 and its companion M51b. We model this outer emission with another morphological component that we call ‘outer disc’ . The outer disc has |$\text{SFR}=0.57\pm 0.06\,$| M_⊙yr⁻¹and |$\Sigma _{\text{SFR}}=5.11\pm 0.56\times 10^{-4}\,$| M_⊙yr⁻¹kpc⁻².

We also show a model with varying dust properties, with the inner disc having a reduced PAH abundance with respect to the main disc. This hybrid model was chosen to solve the overprediction of the standard model at 8 μm in the inner disc, and to alleviate a residual bump in the derived intrinsic SED in the NUV band, at the position of the 2200-Å bump.

ACKNOWLEDGEMENTS

We would like to thank an anonymous referee for very useful and constructive comments, that helped improve the manuscript. Christopher J. Inman acknowledges support from a Science and Technology Facilities Council studentship grant (grant number ST/W507386/1). This work is based in part on observations made with the NASA Galaxy Evolution Explorer. GALEX is operated for NASA by the California Institute of Technology under NASA contract NAS5-98034. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

This work is also based on the Sloan Digital Sky Survey (SDSS) data. Funding for the SDSS IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss4.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics| Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This work has also made use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. We also utilize observations performed with the ESA Herschel Space Observatory (Pilbratt et al. 2010), in particular, to do photometry using the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments.

DATA AVAILABILITY

The data underlying this article are made available at the CDS database via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/MNRAS/

Footnotes

References

APPENDIX A: AMPLITUDE AND GEOMETRICAL PARAMETERS

In our model, we express the amplitude parameters of the stellar discs with spectral luminosity densities L_ν, and the dust discs with central face-on optical depths in the B band, τ_B. Table A1 lists the spectral luminosity densities for each stellar disc and the bulge. Table A2 lists the B band face-on optical depth for each of the morphological components. Table A3 lists other parameters, like the R_in, R_tin, and R_t, which are fixed from data. The table also lists the χ values for the main disc. The latter parameter serves to more cohesively join the inner to the main disc.

APPENDIX B: STAR FORMATION AND DUST MASSES

In this appendix, we present tables with the calculated values for SFR (Table B1), Σ_SFR (Table B2), and dust mass (Table B3). All these values are given for both the global emission and for each morphological component.

APPENDIX C: EXTINCTION CURVES OF THE MILKY WAY AND LMC DUST MODELS

In this paper, we use the Milky Way and LMC dust models given in Weingartner & Draine (2001) in Table 1, line 7 and Table 3, line 1, respectively. Fig. C1 shows the extinction curve for each grain composition as well as the total extinction curve for the Milky Way and LMC dust models. These plots demonstrate that the LMC-type dust exhibits a reduced 2200 -Å bump compared to the Milky Way-type dust. Figs C2 and C3 show the fits at various wavelengths using the hybrid model.