HERschel Observations of Edge-on Spirals (HerOes)
Key Words.:galaxies: ISM - infrared: ISM - galaxies: fundamental: parameters - dust, extinction
We present results of the detailed dust energy balance study for the seven large edge-on galaxies in the HEROES sample using 3D radiative transfer (RT) modelling. Based on available optical and near-infrared observations of the HEROES galaxies, we derive the 3D distribution of stars and dust in these galaxies. For the sake of uniformity, we apply the same technique to retrieve galaxy properties for the entire sample: we use a stellar model consisting of a Sérsic bulge and three double-exponential discs (a superthin disc for a young stellar population and thin and thick discs for old populations). For the dust component, we adopt a double-exponential disc with the new THEMIS dust-grain model. We fit oligochromatic radiative transfer (RT) models to the optical and near-infrared images with the fitting algorithm fitskirt and do panchromatic simulations with the skirt code at wavelengths ranging from ultraviolet to submillimeter. We confirm the previously stated dust energy balance problem in galaxies: for the HEROES galaxies, the dust emission derived from our RT calculations underestimates the real observations by a factor 1.5–4 for all galaxies except NGC 973 and NGC 5907 (apparently, the latter galaxy has a more complex geometry than we used). The comparison between our RT simulations and the observations at mid-infrared–submillimeter wavelengths shows that most of our galaxies exhibit complex dust morphologies (possible spiral arms, star-forming regions, more extended dust structure in the radial and vertical directions). We suggest that, in agreement with the results from Saftly et al. (2015), the large- and small-scale structure is the most probable explanation for the dust energy balance problem.
Cosmic dust severely obscures astronomical objects at ultraviolet and optical wavelengths and, thus, impedes our study of these objects in this waverange. However, with the beginning of the era of the infrared (IR) astronomy, dust began to play an important role in the study of astrophysical processes, taking place in the interstellar medium (ISM), and galaxy evolution.
According to multiple studies (see e.g. Popescu & Tuffs 2002; Viaene et al. 2016, and references therein), approximately one third of the bolometric luminosity in normal spiral galaxies is attenuated by dust: dust reshapes the spectral energy distribution (SED) of galaxies by absorbing and scattering light at short wavelengths and re-emitting at longer wavelengths. Therefore, a dust component, being mainly concentrated in the galactic plane, in optical bands usually appears as a dimmed lane from the edge-on view (e.g. an outstanding example is our Milky Way Galaxy). This component attenuates galactic starlight which passes through it, and, in turn, looks like a very thin emitting disc at far-infrared-submillimeter (FIR-submm) wavelengths.
Until recently, the poor resolution and limited wavelength coverage of infrared instruments was an obstacle to observe dust in even nearby galaxies. With the launch of the Herschel Space Observatory (Pilbratt et al. 2010), we are now able to observe the entire dust peak, in the 70 to 500 m wavelength range (i.e. tracing both warm and cold dust). The high spatial resolution achieved by Herschel allows us to examine the dust in higher detail than ever possible before.
Panchromatic radiative transfer modelling of galaxies provides a powerful tool to analyse different characteristics of dust in galaxies (i.e. optical properties, dust distribution, clumpiness, etc.) in a self-consistent way (see e.g. Popescu et al. 2000; Bianchi 2008; Baes et al. 2010; Popescu et al. 2011; MacLachlan et al. 2011; De Looze et al. 2012a; Holwerda et al. 2012; Mosenkov et al. 2016). The following method has become a standard practice to perform panchromatic dust energy balance studies. From optical observations, the properties and spatial distribution of stars and dust can be constrained using a radiative transfer code to model the propagation of stellar light and its interaction with dust particles in a galaxy. In its turn, the dust emission predicted from the radiative transfer simulations (based on the obtained RT model) can be compared to the observed thermal dust re-emission at IR-submm wavelengths. Such a complementary study imposes that dust features can easily be identified from optical as well as infrared observations. This requirement has limited the number of galaxies for which detailed dust energy balance studies have been attempted in the past. Edge-on spirals are considered ideal cases for studying since projection effects allow to resolve the distribution of different stellar populations (see e.g. van der Kruit & Searle 1981; Morrison et al. 1997; Dalcanton & Bernstein 2002; Tikhonov & Galazutdinova 2005; Seth et al. 2005; Comerón et al. 2018) and dust vertically (see e.g. Kylafis & Bahcall 1987; Xilouris et al. 1999; Misiriotis et al. 2001; Alton et al. 2004; Bianchi 2008; Baes et al. 2010; Popescu et al. 2011; De Looze et al. 2012a; Mosenkov et al. 2016). Thus, using such a multiwavelength technique, a complex 3D model of an edge-on galaxy can be built, which can include both stellar and dust constituents.
Interestingly, multiple dust energy balance studies of individual edge-on galaxies reveal an inconsistency between the predicted FIR/submm fluxes of their radiative transfer models and the observed emission in those wavebands (e.g. Popescu et al. 2000; Misiriotis et al. 2001; Alton et al. 2004; Dasyra et al. 2005; Baes et al. 2010; Popescu et al. 2011; De Looze et al. 2012a; De Geyter et al. 2015a; Mosenkov et al. 2016). Although radiative transfer models might successfully explain the observed optical attenuation, the modelled dust emission underestimates the observed thermal dust re-emission by a factor of 3–4. In order to reconcile the results of the radiative transfer models with the observations, two scenarios have been proposed: either a significant underestimation of the FIR/submm dust emissivity has been argued (see e.g. Alton et al. 2004; Dasyra et al. 2005; MacLachlan et al. 2011) or, alternatively, the distribution of a sizeable fraction of dust in clumps or a second inner dust disc having a negligible attenuation on the bulk of the starlight (see e.g. Popescu et al. 2000; Misiriotis et al. 2001; Bianchi 2008; Bianchi & Xilouris 2011; MacLachlan et al. 2011; De Looze et al. 2012a; Holwerda et al. 2012).
Recently, De Geyter et al. (2015a, hereafter DG15) studied the dust energy balance in two edge-on galaxies with available Herschel observations: IC 4225 and NGC 5166. Using an earlier obtained oligochromatic111modelling based on a small number of images simultaneously. model (using optical data), and adding a young stellar disc to match the ultraviolet (UV) fluxes, they concluded that for NGC 5166 this additional component is necessary in order to reproduce the observed emission at longer wavelengths. However, for the other galaxy, IC 4225, the modelled far-infrared emission still underestimates the observed fluxes by a factor of three. The proposed reasons for this discrepancy in IC 4225 might include 1) that it is too small to properly fit the dust disc parameters, and 2) there is emission from obscured star-forming regions which are embedded in dense dust clouds and, thus, do not contribute noticeably to the observed UV flux but have a clear impact on the FIR emission (see e.g. De Looze et al. 2012a; De Looze et al. 2014; De Geyter et al. 2015a). It is likely that there is no single origin for the dust energy balance problem, and several mechanisms can be responsible for it.
Several years ago, the HEROES project (Verstappen et al. 2013) was initiated to study a set of seven edge-on spiral galaxies, with the aim to present a detailed, systematic and homogeneous study of the amount, spatial distribution and properties of the interstellar dust in these seven galaxies, and investigate the link between the dust component and stellar, gas and dark matter components. These galaxies are well-resolved (with the optical diameter arcmin), and have Herschel observations in five PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) bands at 100, 160, 250, 350 and 500 m. Also, these galaxies were selected to have a clear and regular dust lane, and they have already been fitted with a radiative transfer code based on their optical and near-infrared (NIR) data (Xilouris et al. 1997; Xilouris et al. 1999; Bianchi 2007), therefore an indirect comparison with the previous models can be done. However, we should notice that the used assumptions and fitting strategies in the mentioned studies are different, and, therefore, a uniform approach is necessary to compare the results of RT modelling for all HEROES galaxies.
Mosenkov et al. (2016, hereafter M16) studied the dust energy balance problem in one of the HEROES galaxies IC 2531. For the first time, we used a novel approach for a simultaneous fitting of a set of optical and near-infrared images (7 bands in total) to retrieve a dust model. In it, we showed that the dust properties of IC 2531 are better constrained if near-infrared imaging is used, in addition to optical data. Based on our oligochromatic fit, we constructed a panchromatic222modelling based on a wide range of wavelengths simultaneously. radiative transfer model for the stars and dust in IC 2531, ranging from UV to submm wavelengths. Following DG15, we added another disc component representing a young stellar population to match the FUV flux. Upon comparison of our panchromatic simulations with the Herschel observations at FIR–submm wavelengths, we demonstrated an apparent excess of observed dust emission at those wavelengths for two different dust-grain models: the BARE-GR-S model (Zubko et al. 2004) and the THEMIS model (Jones et al. 2017). The BARE-GR-S model consists of a mixture of polycyclic aromatic hydrocarbons, graphite and silicate grains, while the new THEMIS model is composed of amorphous carbon and amorphous silicates. We showed that the THEMIS model better reproduces the observed emission at the FIR–submm peak. It is interesting to notice, that as in the case of IC 4225 from DG15, the inclusion of the additional young stellar population did not solve the dust energy problem in IC 2531.
In this work, we aim to continue the work started in M16 and analyse the dust characteristics in the six remaining HEROES galaxies, following the same strategy as has been proposed in M16, in a uniform and self-consistent way. For all galaxies in our sample we use available optical and near-infrared observations to constrain their dust models and perform panchromatic simulations. We investigate to which degree our models can predict the galaxy fluxes in different bands, paying particular attention to the dust emission in the FIR-submm domain. This comparison will allow us to answer the question how common the dust energy balance problem might be in galaxies and will give us a clue on a possible reason of its existence.
This paper is organised as follows. Section 2 reviews the sample as well as observing and data reduction strategies. Section 3 presents a brief description of the approach we use in this work to create oligochromatic and panchromatic models. In Section 4, we present results of our study for each individual galaxy and compare our results with the literature. In Section 5, the results of our radiative transfer modelling procedure are discussed in the frame of the dust energy balance problem. The main conclusions of our study are summed up in Section 6. In the Appendix, we provide our models and simulations for each galaxy.
2 The sample
The construction of the sample of the HEROES galaxies is described in detail in Verstappen et al. (2013). Our sample consists of seven relatively nearby edge-on galaxies: NGC 973, UGC 4277, IC 2531, NGC 4013, NGC 4217, NGC 5529, and NGC 5907 (see Fig. 1, which shows their RGB images, and Table 1, which lists some basic characteristics of these objects). Their morphological types suggest that these are late-type galaxies with rather compact bulges, though some of them (NGC 973, IC 2531, NGC 4013, and NGC 5529) clearly exhibit extended boxy/peanut-shape (B/PS) bulges which are often related to bar structures (see e.g. Bureau & Freeman 1999; Chung & Bureau 2004; Méndez-Abreu et al. 2008) and may reside in more than half of the disc galaxies in the Local Universe (Lütticke et al. 2000; Laurikainen & Salo 2016).
The sample galaxies are rather bright and the difference in absolute magnitudes between the faintest and brightest galaxy is not larger than 1.6 mag. All galaxies have similar sizes (in kpc) except for NGC 4013 and NGC 4217 which are about half the size of the others.
As follows from Table 1, all HEROES galaxies are oriented almost exactly edge-on; only NGC 5529 and NGC 5907 have slightly lower inclination (their dust lane is more offset from the major axis). Interestingly, these galaxies differ significantly by apparent axis ratio — IC 2531 and NGC 5907 seem rather thin (0.14 for both) whereas NGC 4013 and NGC 5529 have very thick outermost isophotes (0.38 and 0.30, respectively). Accordingly, the vertical-to-radial extent ratio of stellar discs in these galaxies should differ by a factor 2 to 3. Taking into account that the galaxies are observed at different distances, the angular extent of their dust lanes may also be different. Although the HEROES galaxies were selected by their distinct, regular dust lanes, several objects exhibit somewhat clumpy and patchy dust inclusions. The dust lanes in some galaxies are rather extended and thick (NGC 4217, NGC 5907), in others they are relatively thin and warped (UGC 4277, IC 2531).
For the HEROES galaxies radiative transfer modelling has been done before (though using different modelling strategies): NGC 973 is described in Xilouris et al. (1997, hereafter X97); UGC 4277 – in Bianchi (2007, hereafter B07); NGC 4013 – in X99, B07, and De Geyter et al. (2013); NGC 4217 – in B07; NGC 5529 – in X99, B07; NGC 5907 – in X99 and Misiriotis et al. (2001). The galaxy IC 2531 has been analysed in great detail in M16 (it was also modelled by X99), therefore in this paper we study the remaining six HEROES galaxies. For completeness, the results of our modelling for IC 2531 are also provided.
Notes. The celestial coordinates and morphologies are taken from NASA/IPAC Extragalactic Database (NED). The distances with the corresponding scales are taken from Verstappen et al. (2013), which were found to be the average values of the redshift-independent distance measurements (mostly based on the Tully-Fisher relation). The apparent magnitudes , the semi-major and semi-minor axes ( and ) are measured for the isophote of 25 mag/arcsec. The absolute magnitudes correspond to the adopted distances and apparent magnitudes . For UGC 4277, we used its W1 image to retrieve its magnitudes and semi-major and semi-minor axes. The inclinations of the galaxies are taken from the radiative transfer models from X97 (NGC 973), X99 (IC 2531, NGC 4013, NGC 5529, NGC 5907), and B07 (UGC 4277, NGC 4217). Notice that these inclinations are used in our fitting (see Sect. 3) as a first guess. The fit values of are given in Table 3.
To investigate the dust energy balance problem in a galaxy, we need to collect data over a large wavelength range – from far-ultraviolet to far-infrared/submm regions. For this purpose, we used available databases of optical and NIR sky surveys to extract images of the HEROES galaxies. Here we should notice that we split up our imaging into two sets. The first set of observations was constructed to find fitskirt oligochromatic models, by fitting optical and NIR images (we will call them the reference images) simultaneously. The second set of observations was prepared to supplement galaxy SEDs with fluxes at UV and mid-infrared–far-infrared (MIR–FIR) wavelengths with m and compare the predicted models from our skirt simulations with these observations (see Sect. 3).
Optical observations for NGC 4013, NGC 4217, NGC 5529, and NGC 5907 were retrieved from SDSS DR12 (Eisenstein et al. 2011) in all five bands, with a pixel size of 0.396 arcsec/pixel and an average seeing of 1.1 arcsec. For NGC 973 and UGC 4277, the optical observations were executed on the William Herschel Telescope (WHT) on La Palma, equipped with the Auxiliary-port CAMera (ACAM, Benn & Ellison 1998). The observing run yielded a total of 9 hours, during 3 moonless nights from the 25th till the 27th of December, 2011. The instrument provided a large 8.3 arcmin field of view, with 0.25 arcsec/pixel. The average seeing for the performed observations is 1.32 arcsec. ACAM is equipped with Sloan filters, therefore we used optical observations in the same filters for almost all HEROES galaxies, except for IC 2531 which is located in the southern hemisphere and was observed at the Faulkes Telescope South in the -, - and -band filters (for details, see M16, ).
Four galaxies (NGC 973, UGC 4277, NGC 4013, and NGC 4217) have been observed at the Telescopio Nazionale Galileo (TNG), a 3.58m telescope of the Roque de Los Muchachos Observatory on La Palma. The observations in the band were executed during August and December of 2011 with an average seeing of 1.09 arcsec. The pixel size is 0.25 arcsec/pixel. In addition to the band and for the remaining galaxies we used 2MASS observations in all three bands with a PSF FWHM of 2 to 3 (depending on the atmospheric blurring).
All galaxies, except for UGC 4277, have Spitzer Space Telescope (Werner et al. 2004) observations made by the Infrared Array Camera (IRAC, Fazio et al. 2004) at 3.6 m (IRAC 1). For NGC 973, the mosaic *maic.fits and uncertainties *munc.fits files were taken from post-Basic Calibrated Data, available through the Spitzer Heritage Archive333http://sha.ipac.caltech.edu/. For the other galaxies the observations were downloaded from the Spitzer Survey of Stellar Structure in Galaxies archive (S4G, Sheth et al. 2010; Muñoz-Mateos et al. 2013; Regan 2013; Muñoz-Mateos et al. 2015; Salo et al. 2015)444http://irsa.ipac.caltech.edu/data/SPITZER/S4G/. For UGC 4277, we used the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) in the (3.4 ) band. These optical and NIR observations (from 7 to 9 bands, depending on the galaxy) form our set of reference images. After some processing described in Sect. 2.2 they were used to obtain the galaxy models which can be found in Sect. 3.
The supplementary data include observations from the Galaxy Evolution Explorer (GALEX, Martin et al. 2005; Bianchi et al. 2014) at 0.152 (FUV) and 0.227 (NUV), which were downloaded from the All Sky Imaging Survey reachable trough the GalexView service555http://galex.stsci.edu/GalexView/. WISE observations are available for the entire sky, hence we make use of them in our study as well, in all four (3.4 ), (4.6 ), (12 ) and (22 ) passbands, as an important source of information about the stellar and dust distribution, as well as the star formation rate in galaxies. NGC 973, NGC 5907, and NGC 4013 have additional IRAC 2 (4.5 m), IRAC 3 (5.8 m) and IRAC 4 (8.0 m) observations, while NGC 4217 and NGC 5529 have only IRAC 2 images. For NGC 4013 and NGC 5907, Multiband Imaging Photometer for Spitzer (MIPS, Rieke et al. 2004) observations in all three bands (24, 70, and 160 m) were downloaded from the Spitzer Heritage Archive.
Additionally, we used Herschel PACS and SPIRE imaging to study galaxy SEDs at MIR and FIR wavelengths. To retrieve the Herschel photometer observations, the Herschel Science Archive666http://archives.esac.esa.int/hsa/whsa/ was queried and the observations in all five PACS and SPIRE bands were downloaded.
2.2 Data preparation
All the collected observations described in Sect. 2.1 have been initially reduced and prepared for the subsequent analysis, as it is required by the fitting codes we use in this work (see Sect. 3). The galaxy image processing included several steps: astrometry correction (so that all the frames would have the same astrometry), resampling (rebinning to the frame with the smallest pixel size), determination and subtraction of the sky background, rotating the frames (to align the galactic plane along the horizontal axis), cropping the frames to the minimal size to cover the whole galaxy and at the same time have little empty space. The last step was masking out background and foreground objects, the light of which can contaminate the light of the galaxy. Each step is described in detail in M16.
In short, all these steps are realised with the python Toolkit for skirt (PTS777http://www.skirt.ugent.be/pts/index.html, Camps & Baes 2015, Verstocken et al. in prep.), which contains a lot of functionalities and useful routines. The sky background was fitted by a two-dimensional polynomial of the second order using the unmasked pixels (the initial masks were created by one of the PTS scripts) and then subtracted from the original frames. The borders of the cropped galaxy image were chosen so that the outermost galaxy isophotes, at AB surface brightness (SB) levels of 25.5 mag/arcsec, in all optical and NIR bands would be encompassed. To prepare the images with final masks, we used the segmentation maps created with the sextractor package (Bertin & Arnouts 1996) and then revisited the created masks by hand.
In addition, creating the accurate point spread function (PSF) for every galaxy image in each band is important since atmosphere and telescope blurring can seriously affect the light coming from the galaxy and, thus, change the real SB distribution of the observed object. To determine the PSF in each frame (if applicable), we selected several good stars (not crowded, with a high signal-to-noise ratio). Then they were fitted with a Moffat (Moffat 1969) (for the SDSS, TNG and WHT observations) or Gaussian (for the 2MASS and WISE images) function. For the IRAC images, we used in-flight point response function (PRF) images888http://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/calibrationfiles/ for the centre of the IRAC 3.6 m fields downsampled and re-rotated to correspond to the analysed galaxy frames. For the other bands, we used the kernels provided by Aniano et al. (2011). The final PSF images were rebinned and rotated to match the final galaxy images.
The reduction of the Herschel PACS and SPIRE data was performed in Verstappen et al. (2013), however, we repeated it following the up-to-date recipe described in Davies et al. (2017) and Clark et al. (2018) for the DustPedia sample. For SPIRE, we applied HIPE v13999http://www.cosmos.esa.int/web/herschel/hipe-download and the Bright Galaxy Adaptive Element (BriGAdE; Smith et al. 2012) pipeline. Final SPIRE maps were produced using the HIPE v13 naïve map-maker, with pixel sizes of 6, 8, and 12 arcsec at 250, 350, and 500 m respectively. For PACS data, HIPE v13 and Scanamorphous v24 pipeline (Roussel 2013) were used. The final PACS maps have pixel sizes of 3 and 4 arcsec at 100 and 160 m respectively.
All the images from the reference and supplement sets were prepared in the same way as described above.
Integrated flux densities for all wavelengths were determined by means of aperture photometry from the photutils python package, which was applied to the sky subtracted frames. For each galaxy, we determined an elliptical aperture which is 1.5 times larger than the outermost isophote (at the level) among all the prepared galaxy images. For each galaxy, we used a uniform mask for all reference images and individual masks for supplemental observations. We calculated the total flux inside the aperture by summing the flux value of each pixel; the masked areas were filled with the interpolated values using a Gaussian weights kernel.
To estimate the uncertainties, we calculated the total uncertainty by summing in quadrature the calibration uncertainty (for the surveys we used the values listed in tab. 1 from Clark et al. 2018, for the other images we used a calibration uncertainty of 5%) and the uncertainty on the sky level.
We also added IRAS (Moshir et al. 1990; Sanders et al. 2003; Lisenfeld et al. 2007) flux densities at 25, 60 and 100 m taken from NED, and Planck (Planck Collaboration et al. 2014) APERFLUX values at 353, 545 and 857 GHz (which correspond to 850 , 550 and 350 , respectively) from the Planck Legacy Archive101010http://pla.esac.esa.int/pla. The two brightest galaxies in our sample, NGC 5907 and NGC 4217, were detected by ISO as part of the ISOPHOT 170 m Serendipity Survey (Stickel et al. 2004). For NGC 5907, we found two fluxes at 450 and 850 from the Submillimetre Common-User Bolometer Array (SCUBA, Holland et al. 1999; Stevens et al. 2005). We did not include the Akari/FIS (Murakami et al. 2007) 160 m flux densities found for four HEROES galaxies in our study because the comparison with the PACS 160 fluxes is very poor (to estimate the fluxes they used PSF photometry with a FWHM of about 60 arcsec at 160 , while all the galaxies in our sample are much larger).
We provide all flux densities (measured and taken from literature) and their uncertainties in Table 2.2.
Survey NGC 973 UGC 4277 IC 2531 NGC 4013 NGC 4217 NGC 5529 NGC 5907 () 0.15 — — 0.23 — SDSS 0.36 — — — WHT 0.36 — — — — — FTS 0.45 — — — — — — SDSS 0.47 — — — WHT 0.47 — — — — — FTS 0.55 — — — — — — SDSS 0.62 — — — WHT 0.62 — — — — — FTS 0.66 — — — — — — SDSS 0.75 — — — WHT 0.75 — — — — — SDSS 0.89 — — — WHT 0.89 — — — — — 2MASS 1.25 2MASS 1.65 2MASS 2.20 TNG 2.13 — — — WISE 3.4 IRAC 3.6 — IRAC 4.5 — — WISE 4.60 IRAC 5.8 — — — — IRAC 8.0 — — — — WISE 12.1 IRAS 12 — — WISE 22.2 MIPS 24 — — — — — IRAS 25 — — IRAS 60 — MIPS 70 — — — — — IRAS 100 — PACS 100 MIPS 160 — — — — — PACS 160 ISO 170 — — — — — SPIRE 250 SPIRE 350 Planck 350 — — SCUBA 450 — — — — — — SPIRE 500 Planck 550 — — — Planck 850 — — — SCUBA 850 — — — — — —
3 Fitting strategy and simulations
Our strategy is described in detail in M16. The general idea is as follows. First, we fit a model, which consists of several stellar components, e.g. a Sérsic (Sérsic 1968) bulge and thin and thick double-exponential discs, to a NIR image. IRAC 3.6 observations are ideally suited for solving this task since images at these wavelengths serve as good relative proxies for stellar mass (Elmegreen & Elmegreen 1984; Grosbol 1993; Rix & Zaritsky 1995; Meidt et al. 2012, 2014; Querejeta et al. 2015) and the influence of dust is here significantly diminished. For UGC 4277, there are no Spitzer observations. Therefore, in this case we use the TNG -band image, which also offers a good opportunity to trace the distribution of old stellar populations. The visual inspection of these NIR images (in the 3.6m and bands) did not reveal any apparent traces of a dust lane in these galaxies. Therefore, we can consider, as a first approximation, the attenuation at these wavelengths to be zero.
Once the stellar model is obtained, we can use the fitskirt111111http://www.skirt.ugent.be/skirt/index.html code (De Geyter et al. 2013, 2014). It allows one to find structural properties of the dust component by means of skirt (Baes et al. 2011; Camps & Baes 2015) oligochromatic radiative transfer simulations and with help of a genetic algorithm-based optimisation (Goldberg 1989). At this stage, a set of optical and NIR observations can be fitted simultaneously, together with a dust component which is often represented by a double-exponential disc. Here, we keep fixed the geometrical sizes of the stellar components to the values retrieved in the first step, but allow to change their relative luminosities. This is a relatively fast and effective way to fit the dust model and maximally avoid the problem of the local minimum of , by reducing the number of free parameters when searching for the best model fit.
The third step is performing skirt panchromatic radiative transfer simulations using the obtained fitskirt models. Having them in hand, we can create mock observations (we call them further ‘panchromatic simulations’) and, therefore, predict how the galaxy would look like at different wavelengths, from FUV to submm, and then compare the real and model SEDs, as well as their real and mock images, and conclude how accurately our model is able to describe the observations at specific wavelengths.
Here we should clarify why we use this three-step approach instead of directly fitting all available FUV-submm observations simultaneously to derive the parameters for both stellar and dust components. First, in the present version of fitskirt, all wavelengths are considered to be independent. Therefore, the code is only applicable to wavelengths up to NIR. Further, with increasing wavelength, the dust emission becomes an important part of the simulation, and, hence, this assumption is no longer valid. Second, the computational cost of such multiwavelength fitting is very high because not only the number of free model parameters substantially increases, but also the amount of data, which should be analysed simultaneously during the fitting.
3.1 Decomposition of the IRAC 3.6 images
The surface brightness profile of a stellar bulge is usually fitted with a Sérsic model:
where is the effective surface brightness, i.e. the surface brightness at the half-light radius of the bulge , and is a function of the Sérsic index (Caon, Capaccioli & D’Onofrio 1993; Ciotti & Bertin 1999). An additional free parameter is the apparent axis ratio of the bulge . We also apply a Sérsic profile for describing SB of a halo for one HEROES galaxy, NGC 4013.
The three-dimensional axisymmetric stellar disc can be described by a double-exponential function , which in a cylindrical coordinate system aligned with the disc (where the disc mid-plane has ) is given by
where is the central luminosity density of the disc, is the disc scale length, and is the vertical scale height. In this work we use a modified profile that includes a break radius (it can also be called truncation/anti-truncation radius). Erwin et al. (2005, 2008), Muñoz-Mateos et al. (2013), and M16 have shown that for some galaxies this profile better describes the apparent SB distribution (Erwin et al. 2005; Erwin 2015):
In this formula, parametrises the sharpness of the transition between the inner and outer profiles with the break radius . Large values of the sharpness parameter () correspond to a sharp transition, and small values () set a very gradual break. The dimensionless quantity is a scaling factor, given by
In total, this model contains five free parameters: the scale length of the inner disc , the scale length of the outer disc , the scale height , the break radius , and the sharpness of the break .
Visual inspection of all minor-axis profiles of the HEROES galaxies (see Appendix A) suggests the presence of both a thin and a thick disc. Almost all the galaxies (except for UGC 4277) show breaks on their radial profiles, therefore for them we used the model (3), and for the remaining galaxies the simple model (2) was applied. As such, our uniform model generally consists of a thin and a thick broken disc and a central Sérsic bulge, except that for NGC 4013 the model also includes a Sérsic stellar halo (see Sect. 4.4).
In our study we make use of the imfit121212http://www.mpe.mpg.de/erwin/code/imfit/ code (Erwin 2015) which works with different 3D geometries (including the model of a broken double-exponential disc (3), which is realised in the imfit BrokenExponentialDisk3D function with in Eqs. (40) and (41) in Erwin 2015). To minimise the number of free parameters, we fixed the inclinations for both the thin and thick discs. These inclinations were initially found from the radiative transfer modelling of the HEROES galaxies with a simple ‘disc+bulge’ stellar model and a dust disc (described in the PhD thesis by De Geyter 2015). As shown in Table 3, which lists the results of our RT modelling (see Sect.3.2), the inclinations adopted here are essentially the same as found in our subsequent oligochromatic fitting. Also, we fixed the sharpness of the breaks at , since by varying its value we did not find a significant difference in the decomposition results. In addition, the break radii should have the same value for both the thin and the thick disc.
We wrote a special python wrapper to apply genetic algorithms in our process of searching for the best-fit model. To estimate uncertainties of the free parameters in our model, we applied the genetic algorithm ten times and took the scatter in the fitted parameters as their uncertainties.
|Parameter||unit||NGC 973||UGC 4277||IC 2531||NGC 4013||NGC 4217||NGC 5529||NGC 5907|
|1. Thin disc:|
|2. Thick disc:|
Notes. The inclination angle was fixed for both the thin and thick discs. We used the BrokenExponentialDisk3D function for the thin and/or thick disc if the respective is defined in the table, otherwise the ExponentialDisk3D was used. The bulges and the halo were fitted with the Sérsic function. To compute the total stellar mass from the total luminosity in the respective NIR band, we relied on the mass-to-light ratio at 3.6 of Eskew et al. (2012), except for UGC 4277 for which the mass-to-light ratio in the band from McGaugh et al. (2000) was used.
3.2 fitskirt fitting
fitskirt is a fitting code that combines the output of skirt with a genetic algorithm optimisation library to retrieve best-fitting model parameters for the stellar and dust components in a galaxy. The results retrieved from the 3D fitskirt model should reproduce the fitted images of the galaxy, taking the effects of dust obscuration and multiple anisotropic scattering fully into account.
The automated fitting routine fitskirt has been tested on a mock (skirt simulated) galaxy image (De Geyter et al. 2013). Also, its capabilities have been validated by applying it to a dozen of real edge-on galaxies (see De Looze et al. 2012b; De Geyter et al. 2014; Saftly et al. 2015; Viaene et al. 2015; Mosenkov et al. 2016). In these works, it was found that fitskirt is able to give reasonable constraints on all free parameters describing the stellar disc, stellar Sérsic bulge, dust disc, and even dust ring.
De Geyter et al. (2014) showed that oligochromatic fitting, i.e. fitting applied to a number of bands simultaneously, has clear advantages over monochromatic fitting in terms of accuracy. In particular, the parameters which describe the dust distribution, have a smaller spread as the oligochromatic fitting method is less prone to degeneracies in the free parameters. In addition to that, M16 demonstrated for the example of one HEROES galaxy IC 2531 that the optical+NIR data based model is much better constrained than the model where only the optical images are fitted. This is why in this study we collected available optical and NIR data in a broad range of wavelengths.
To describe the dust constituent in our fitting, we apply a model of a double-exponential dust disc (see e.g. Xilouris et al. 1999; Baes et al. 2010; Verstappen et al. 2013; De Geyter et al. 2013, 2014):
where is the total dust mass, is the radial scale length, and is the vertical scale height of the dust disc. To calculate the central face-on optical depth, one should use the following expression:
where is the extinction coefficient of the dust. The central edge-on optical depth can be calculated as . If the stellar disc has a truncation at some radius, the dust disc is also truncated at this radius.
To find the dust model parameters, we fix the geometry of the HEROES galaxies built upon the IRAC 3.6 data models, which were described in Sect. 3.1. We allow only the dust disc parameters, the alignment of the galaxy centre, the inclination angle , and the galaxy position angle to change during the fitting. At the same time, the luminosities of the stellar components are determined individually at each wavelength. This simulates the wavelength-dependent behaviour of luminosity ratios of different stellar components (e.g. the bulge-to-total luminosity ratio), which are included in the model (De Geyter et al. 2013, 2014).
We should stress here that we do not make any assumptions about the characteristics of the stellar populations, and about how the emission in different wavebands is linked. We merely fit the stellar emission in every band individually (but assume the same dust distribution at each wavelength). By doing so, we significantly reduce the number of free parameters and, at the same time, use a complex model to describe the stellar and dust components of the galaxy.
In our fitskirt fitting we use the new THEMIS model, which is described in detail in Jones et al. (2017) and implemented in skirt by Camps et al. (2015). This model is completely built on the basis of interstellar dust analogue material synthesised, characterised and analysed in the laboratory. In this model, there are two families of dust particles: amorphous carbon and amorphous silicates. For the silicates, it is assumed that 50% of the mass is amorphous enstatite, and that the remaining half is amorphous forsterite. This model is becoming a new paradigm, and is being developed in the frame of the DustPedia project.
The fitskirt computations were done using a multi-core computer and the high performance cluster of the Flemish Supercomputer Center (Vlaams Supercomputer Centrum). For each galaxy, we used about 50 processor units when running fitskirt for one galaxy with the following specification: photon packages, 200 individuals in a population, and 100 generations (this has been tested to be optimal values by De Looze et al. 2012b; De Geyter et al. 2014). To estimate the uncertainties on the parameters, we repeated the fitting five times. Taking into account that the average processing time of one fit was about 70 hours (up to 1 week, depending on the dimensions of the reference images), the total time for the whole sample was approximately 122500 CPU hours, or 2450 hours of work spread among 50 cores.
The results of the fits for all seven HEROES galaxies can be found in Table 3.
|Parameter||unit||NGC 973||UGC 4277||IC 2531||NGC 4013||NGC 4217||NGC 5529||NGC 5907|
3.3 skirt simulations
To perform a uniform and detailed energy balance study of the HEROES galaxies, one important step is needed — we should extend our oligochromatic models from Sect. 3.2 to panchromatic models. We did this to not only reproduce the images at optical and near-infrared wavelengths, but also to restore the entire UV–submm SED and create panchromatic simulations which can be compared to the real ones. As such, our approach consists of two stages.
First, we directly use the oligochromatic fitskirt model obtained earlier in Sect. 3.2 in order to predict the emission of the galaxy in the entire UV to submm domain.
Panchromatic simulations imply that the properties for both the stars and the dust need to be set over the entire wavelength domain. For the stellar components, we assume a Bruzual & Charlot (2003) single stellar population SED with a Chabrier (2003) initial mass function and a solar metallicity (). For all galaxies in our sample we use the same ages of the stellar components as we determined for IC 2531: for both the thin and the thick disc, they are about 5 Gyr, and 8 Gyr for the bulge (we discuss this in Sect. 5.2). The stellar halo for NGC 4013 was assumed to be 10 Gyr old.
As follows from Popescu et al. (2000); De Looze et al. (2012a); De Looze et al. (2012b), DG15 and M16, an additional source of UV luminosity, i.e., a young stellar component, is needed in order to match the observed and model SEDs. As UV radiation is easily absorbed by interstellar dust, the addition of a young population will also affect the dust emission at MIR to submm wavelengths. Therefore, at the second stage we add a young stellar disc, which is completely obscured by dust in the reference images, with the following properties. In our models, we assume a scale height of of the scale height of the dust disc and the same scale length as for the thin disc. However, as shown in De Looze et al. (2014), DG15, M16, and Viaene et al. (2017) varying the thickness of the young stellar disc between one third of the dust scale height to one dust scale height barely affects the resulting SED.
The observed GALEX FUV image traces unobscured star formation regions of the galaxy. The stellar emission spectrum of these stars is described through a starburst99 SED template which represents a stellar population with a constant, continuous star formation rate (SFR) and an evolution up to 100 Myr (Leitherer et al. 1999). In this case, the initial mass function is a Salpeter (1955) IMF with masses between 1 and 100 M and with solar metallicity. The luminosity of this component is constrained by the GALEX FUV (or NUV in the case of UGC 4277) flux density.
The simulated SEDs for both models with (except for NGC 973, which has very noisy GALEX observations) and without a young stellar population are shown in Appendix C. The panchromatic simulations, together with the real observations, are presented in Appendix D.
Below we discuss the results individually for each galaxy.
4.1 Ngc 973
This galaxy classified as an Sb spiral (NED) is the most massive galaxy in our sample, with a radius of the outermost isophote of 25 mag/arcsec up to 40 kpc. It is seen almost exactly edge-on; the dust lane is fairly regular, without strong bending or clumpiness. A clear thick B/PS/X structure in the centre is a distinctive feature of the galaxy (Lütticke et al. 2000) which points to the presence of a strong bar (Patsis & Xilouris 2006). Also, it is classified as a Sy2/LINER object with significant X-ray emission detected by SWIFT and INTEGRAL (see e.g. Sazonov et al. 2007; Weedman et al. 2012). According to Allaert et al. (2015), this galaxy does not show strong disc warping in Hi, nor in optical observations.
Our imfit model (see Fig. 2) for this galaxy includes three stellar components: thin and thick discs with particularly prominent radial breaks and a bulge (which obviously should be interpreted as a ‘bulge+bar’ component). Surprisingly, the thin disc is remarkably faint compared to the thick disc (). The bulge/bar component fitted with the Sérsic function is flat and has an extended surface brightness profile (the Sérsic index is large, ) with an ordinary effective radius. The galaxy shows the largest fraction of the bulge/bar to the total galaxy luminosity among all sample galaxies ().
The fitskirt model, which is based on the obtained IRAC 3.6 model for the stellar components plus the dust exponential disc, describes the observations remarkably well, even despite the X-shape structure in the residual images (see Fig. 3). The ISM appears to be rather optically thick, even in the face-on orientation (); by dust mass this galaxy is one of the most massive galaxies in our sample (excluding NGC 5907). Unfortunately, since this galaxy is rather distant, the GALEX observations are very poor for this galaxy, therefore we did not include and fluxes in the global SED. Thus, for this galaxy we did not add a young stellar disc to the model. As follows from Fig. 4, our panchromatic simulations somewhat underestimate the observed fluxes in the MIR and FIR domain – at 160 m, the model predicts 75 per cent of the observed flux. However, this is not as much as for the other galaxies, as we will see below. The panchromatic simulations and observations look highly similar in the optical and NIR bands (see Fig. 5). The WISE , PACS 100 and PACS 160 images reveal a possible spiral arm, or another structure with a high dust concentration, highlighted by a high emission region, which is hard to discern in the SPIRE images. In comparison with our model, all Herschel observations show an elongated, radially extended structure, which obviously makes a significant contribution to the total emission in these bands.
4.2 Ugc 4277
UGC 4277 is a flat Sc/Scd galaxy, one of the largest galaxies (the semi-major axis kpc) in our sample. It is seen almost exactly edge-on; the dust lane is very prominent in the optical. One of the noticeable features is an integral-shaped optical warp. Allaert et al. (2015) recently showed that the outer regions of the Hi disc demonstrate a small warp as well, with the position angle deviating at most about 4 degrees from its central value. UGC 4277 has a total atomic hydrogen mass of and, similarly to NGC 973, the Hi disc has a clumpy and irregular structure. Interestingly, this galaxy is included in the Catalogue of Isolated Galaxies (CIG, Karachentseva 1973), therefore its evolutionary status may differ from those galaxies residing in a group or a cluster.
For the imfit modelling (which was based on the -band image since there is no Spitzer data for this galaxy), we used an ordinary (with no breaks since they are likely to be smoothed because of insufficient angular resolution) double-exponential profile for both the thin and thick disc and a Sérsic bulge (see Fig. 7). Both discs are comparable by luminosity and radial extent. The bulge has a rather large effective radius with a Sérsic index of about 2.6.
The fitskirt model looks exceptionally good – only the disc warp is noticeable in the residual images (see Fig. 13). The dust attenuation is represented particularly well for the whole set of reference images. It is interesting to note that, photometrically, the bulge shows no signature of an X-shape structure in the centre. However, we can see a prominent B/PS structure in the residual images in the , and bands. This is evidence that UGC 4277 has a bar which is not so extended as in the other HEROES galaxies or it is not viewed purely side-on (Bureau et al. 2006).
Unfortunately, the GALEX FUV observation of this galaxy is too noisy, therefore we used the NUV flux to match the galaxy observation and the model (with an additional young disc) in this band. The model with the young stellar disc (see Fig. 20) describes the observations in the MIR range significantly better, but the discrepancy in the FIR is still present, with an underestimated cold dust emission at wavelengths larger than 200m, by a factor of about 2. Consequently, the panchromatic simulations cannot entirely describe the observations. The panchromatic simulations in the optical and MIR are compatible with the observations (see Fig. 26). However the comparison between the WISE and the Herschel observations with the models indicates that the real observations are again more radially extended, as in the case of NGC 973. Other structural details (e.g. spirals), which would be present in this galaxy, are washed away because of poor resolution.
4.3 Ic 2531
This galaxy has been described in great detail in M16, therefore we refer the interested readers to that article.
4.4 Ngc 4013
NGC 4013 is the most compact disc galaxy in our sample. It is of the Sb type, with several outstanding features which were studied by different authors. The galaxy has a barely visible optical warp (Florido et al. 1991). Martínez-Delgado et al. (2009) found a giant, arcing stellar structure of low surface brightness around the galaxy. Comerón et al. (2011) produced a decomposition considering three flattened stellar components in this galaxy: thin disc + thick disc + halo.
The atomic gas content of NGC 4013 was analysed in great detail by Zschaechner & Rand (2015) with similar techniques as used in Allaert et al. (2015). They found that the gas disc of this galaxy is strongly warped, both along the major axis (in the plane of the sky) and along the line of sight, and shows a lag (decreasing rotation velocity with increasing height above the plane) that is strongest at small radii. Finally, they found no Hi gas at large distances from the plane, despite the presence of extra-planar dust and ionised gas (see e.g, Miller & Veilleux 2003, and the references therein).
Since this galaxy has a bright stellar halo, we used a four-component model ‘thin disc + thick disc + bulge + halo’ (see Fig. 9). The two discs have a broken double-exponential profile, the bulge and the halo are described by a Sérsic function. The thick disc is the most luminous component of the galaxy, whereas the bulge is compact and round. A foreground star almost coincides with the position of the nucleus of the galaxy, which is also covered by a dark dust lane. Ho et al. (1997) found that this galaxy hosts an AGN in its centre, therefore our bulge profile should contain it as well. The halo is indeed noticeably bright (its fraction to the total luminosity is 15%) and extended ( kpc), the Sérsic index is close to 1, therefore it can be called the second thick disc.
Since our fitskirt model does not take into account the X-shape structure, this feature seems a prominent detail in the residuals in Fig. 15. Despite this and some unfitted clumped structures on the left side of the galaxy plane, our fitskirt model describes the observations fairly well (the SDSS band image is obviously too noisy and, thus, the comparison in this band is meaningless).
Our panchromatic simulations (see Fig. 21 and 27) report the same dust energy balance problem – even with an additional young stellar component, the model emission in the MIR/FIR domain is far below the observed one.
Since a lot of data are available for this galaxy, it is interesting to compare the observations and the model in detail. In the FUV and NUV bands, we can see that our galaxy is observed as a patchy disc consisting of a young stellar population with star forming regions, contrary to our model which looks more like an attenuating lane. The band model is closer to the corresponding observation. All other optical–4.6m observations are described exceptionally well by our model. The MIR–FIR observations reveal a probable ring structure (since the inclination angle is almost exactly , we can only assume this by the typical “shoulders” for rings on the SB profile) in the plane of the galaxy which strongly adds to the emission in this wavelength range, for the hot as well as for the cold dust.
4.5 Ngc 4217
NGC 4217 is another Sb-type galaxy in our sample. Visually, it has a noticeably thick (as compared to the other HEROES galaxies) and somewhat warped dust lane, which extends up to the galaxy outskirts.
For the main disc of NGC 4217, Allaert et al. (2015) measured an atomic hydrogen mass of . In addition, they also detected a separate ring of gas outside (and coplanar with) the main gas disc, with an Hi mass of at least . This ring has no optical counterpart and could be the remnant of a recent minor merger. The main gas disc is mostly flat and shows only a small warp in the outer parts where it meets the ring structure.
Our imfit fitting was performed for the model consisting of two broken thin and thick discs and a bulge (see Fig. 10). As one can see from Table 2, the thin disc is less luminous than the thick disc, but they have similar inner disc scale lengths. The bulge is surprisingly compact (the Sérsic index is small) and its fraction to the total galaxy luminosity is minor. The thick disc has the smallest disc scale length-to-scale hight ratio 3.8 – apparently, this galaxy looks thicker than the others.
The fitskirt model is quite adequate with regard to the observations (see Fig. 16). Only some traces of the dust attenuation are left in the residuals in the , and bands which may be related to the patchy nature of the dust disc in this galaxy.
Again, the skirt simulations cannot describe the real observations in the whole MIR-FIR domain (see Fig. 22). The typical factor of flux underestimation is 4 to 5. However, the comparison of the panchromatic simulations and real observations for this galaxy is very good (see Fig. 28), even for the FUV and NUV bands. The WISE and observations reveal some asymmetric structures around the centre, which may be related to a spiral arm. Despite the discrepancy between the observed and model SED, the structure of the emitting cold dust in the Herschel observations is particularly well described by our model. It is likely that the characteristic features, which we saw in the FIR in the previous galaxies, are not spatially resolved here.
4.6 Ngc 5529
NGC 5529 is an Sc edge-on galaxy with slightly lower than inclination (). A prominent optical warp and B/PS structure in the centre are the characteristic features of this galaxy.
It is part of a rich galaxy group with at least 16 other members (Irwin et al. 2007). Kregel et al. (2004) already studied the atomic gas content of this galaxy and discovered that NGC 5529 is connected to two of its satellites through Hi bridges. They also found a bright Hi ridge in the major axis position-velocity diagram, which could indicate the presence of a strong spiral arm. On the approaching (NW) side the atomic gas disc shows the same warp as the stellar disc, while the receding (SE) side of the gas disc appears more strongly warped. Allaert et al. (2015) also found a moderate warp along the line of sight and detected a radial inflow of gas at about 15 km/s in the outer half of the disc. The latter is probably linked to the accretion of gas from the satellites. Allaert et al. (2015) measured a total atomic hydrogen mass of .
Our imfit model consists of a thin broken disc, a thick disc, and a bulge (see Fig. 11). The thin and thick discs have an equal contribution to the total luminosity. The bulge/bar structure is rather extended, with the Sérsic index larger than 2 – a hint that this galaxy may have a classical bulge.
Despite the disc warp (which was masked out during the fitting), the fitskirt decomposition looks quite satisfactory (see Fig. 17), with some traces of the dust content in the residuals in the and bands. However, the skirt simulations (see Fig. 23) reveal a systematic underestimation of flux densities at all wavelengths larger than 10m by a factor of 4 to 5, for the models both with and without the young disc. The comparison of the panchromatic simulations with the real galaxy images (see Fig. 29) indicates that starting from the WISE band with increasing wavelength, the 2D SB distribution of the real observations becomes more complex, with round isophotes (whereas our models have clear discy isophotes at all wavelengths) and apparent existence of spirals (which are distinct emission features on both sides from the centre, from WISE to SPIRE 250). Interestingly, our NUV model does not look similar to what we see in real – the absorption lane in the panchromatic simulation versus the emission stripe in the observation.
4.7 Ngc 5907
NGC 5907 is the nearest and flattest HEROES galaxy, also of Sc morphological type. NGC 5907 is seen within a few degrees of edge-on (). Its tiny bulge is almost hidden behind the opaque dust lane. Martínez-Delgado et al. (2008) discovered extremely faint extended arcing loops around the galaxy, which, as they argue, were produced by low-mass satellite accretion in the hierarchical formation scenario. Another interesting feature of this galaxy is an excess of light at the periphery of the minor axis profile, which was interpreted by Sackett et al. (1994) as an extended faint halo. Sandin (2014), however, showed that this excess can also be explained by diffuse scattered light produced by the extended PSF (we discuss this issue in Sect. 5.1).
The atomic gas disc of NGC 5907 is strongly warped and has a total Hi mass (as measured by Allaert et al. 2015) of . This warp was already noticed by Sancisi (1976). Later, Shang et al. (1998) also discovered a small companion galaxy toward the NW of NGC 5907 in their Hi map. Sancisi (1976) further found that the gas disc of NGC 5907 is also warped along the line of sight and shows strong asymmetries between the approaching and receding side. An extended reservoir of low-density gas was detected around the main disc.
Our stellar model for this galaxy is ‘broken thin disc + thick disc + bulge’ (see Fig. 12). The thin disc dominates over the thick disc, which is indeed particularly vertically extended and faint. The small bulge is exponential and flattened, and therefore can be referred to as a pseudo-bulge (Gadotti 2009; Fisher & Drory 2010).
The fitskirt model and the real observations are plotted together in Fig. 18. In the panel with residuals one can see that for the optical data, our fit does not perfectly restore the observation – a clear residual lane is well visible. However, the overall model looks good and compatible with the observations.
Nevertheless, if we look at the parameters listed in Table 3, the fitskirt model for NGC 5907 differs from the others. From Table 3 we can see that the dust mass and the face-on optical depth are the largest among all HEROES galaxies. It follows from our modelling that the ISM of this galaxy is obviously non-transparent. However, if we look at Fig. 24, we can see that our panchromatic simulations systematically underestimate flux densities in the optical domain, whereas at FIR wavelengths our model with an additional young disc is able to predict the SED behaviour exceptionally well (this is the only HEROES galaxy for which our model follows the observations in the FIR; our panchromatic simulations are compatible with the real observations as well, see Fig. 30). However, the serious discrepancy in the optical bands suggests that the fitskirt model is incorrect due to an overestimation of the dust mass. Perhaps, one of the possible explanations for this erroneous model is that our stellar model is not correct either (e.g. the thick disc is, in fact, not real and produced by scattered light of the PSF).
Interestingly, in the IRAC 4.5 m band image there is a prominent detail above the plane of the disc. However, since it is not seen in the other bands, including the WISE band close to IRAC 4.5 m, we interpret this feature as an artefact.
4.8 Comparison with other works
The direct comparison of our results with other works is impossible because in this study we applied a new technique which is based on the two-step decomposition. This novel strategy has not been applied before, except in M16. However, for uniformity, we provide some indirect comparison with other authors who carried out their own decomposition analysis of the HEROES galaxies in the same band(s).
Our imfit models can be compared to the SG results (Salo et al. 2015) for four HEROES galaxies: NGC 4013, NGC 4217, NGC 5529, and NGC 5907. However, some caution should be made. Firstly, our decomposition is performed for 3D stellar discs, whereas the GALFIT code (Peng et al. 2010), which was used in the SG project, works only with projected intensities. Secondly, for these galaxies the SG models do not contain a bulge component (for NGC 4013, NGC 4217, and NGC 5529, the central component is a nucleus). In our modelling, all our galaxies have a bulge with an effective radius larger than the FWHM of the PSF (where for IRAC 3.6m the Gaussian PSF FWHM is ). So, even being compact, the bulges in our models should have physical meaning. Thirdly, the GALFIT decomposition from Salo et al. (2015) assumes an isothermal disc geometry with the disc height scale , whereas in our decomposition we use a model of a double-exponential disc, although we can compare the height scales of both the exponential and isothermal disc for vertical distances : . The results of our and the SG decomposition are presented together in Table 4. From the table we can see that the retrieved disc scale lengths are similar, whereas the disc scale height may differ by a factor of 0.5–2, which can be explained by the reasons listed above.
Since radiative transfer models have been obtained before for all HEROES galaxies, it is interesting to compare the results of this work with the literature. Table 5 lists the results taken from X97, X99, and B07. In these works, galaxies have a simple ‘bulge+disc’ model plus a dust disc. The shape of the bulge profile (defined by the Sérsic index) was not fitted. The number of bands, which were used in the analysis, is limited – from two to five (optical and near-infrared). The fitting procedures are also different in each work. Therefore, this comparison is indirect and serves only for the completeness and reliability of our work. From the comparison, we can see that the inclination angle is very close for all the works, with the exception of NGC 5907. For this galaxy, our results slightly differ from X99. The comparison between our results and B07 for NGC 4217 and NGC 5529 is very good, for all dust parameters. For NGC 973, our model is very different from that in X97: our dust disc is 50% less extended and thinner, with a higher optical depth. For UGC 4277, the comparison for the dust scales is satisfactory. For NGC 4013, our study, B07, and De Geyter et al. (2013) provide similar results for the dust scale length and scale height. We should notice that the optical depth for almost all galaxies is rather different in all compared works. The same is true for the bulge-to-total luminosity ratio . In addition, it is particularly hard to estimate the uncertainties on these parameters.
Another reason for the differences with previous results must be also due to differences in the dust models for the optical constants of the dust grains which were used previously in the literature and in the current work. For example, Bianchi (2007) used the dust scattering properties from the Milky Way dust grain model of Weingartner & Draine (2001). In this paper, we applied the new THEMIS dust-grain model which is built upon the laboratory-measured properties of interstellar dust analogue materials. Obviously, even if all dust models are tuned to fit the observed extinction curve of the Milky Way, there will be a different contribution from the different grain components, which constitute these models, to the scattering and absorption, resulting in some changes in the output RT models. This makes the comparison of the RT modelling results even more difficult.
On the whole, taking into account that the works in our comparison used different stellar and dust mixture models, the results for the dust component are compatible.
Notes. (1) [Table 2 (de Vacouleurs bulge),]1997A&A…325..135X, (2) X99, (3) B07, (4) De Geyter et al. (2013). Inclination angle is given in degrees. Geometrical scales and are given in kpc (we reconverted original values from these works using the scales from Table 1). For our study, (except IC 2531) were calculated by means of interpolation at m. For (1) and (2), were calculated using the values given in their table 2 and tables 1–6, respectively.
5.1 On our stellar models
Four galaxies in our sample have quite prominent B/PS structures which are not described in our modelling. As was shown in many studies, B/PS features appear because of the buckling instability of a galactic bar (Combes & Sanders 1981; Combes et al. 1990; Pfenniger & Friedli 1991; Bureau & Freeman 1999; Athanassoula & Misiriotis 2002), and, thus, they point to the presence of another stellar component in the galaxy, a stellar bar. In principle, this might affect the results of our decomposition, especially concerning the thick disc and the bulge. However, as we have stressed before, the bulges in our study are, in fact, a superposition of at least two components at the centre – a bulge and a bar (if it is present in the galaxy). Also, as the bulge/bar structure is less extended than the stellar discs (we consider mainly late-type spirals with a small contribution of the central component to the total galaxy luminosity), its influence on the retrieved disc parameters should not be significant.
Recently, a 2D Ferrers profile (Ferrers 1877) with a modulated Fourier mode has been applied to describe B/PS structures in 22 edge-on galaxies (Savchenko et al. 2017). A multicomponent decomposition, using which they were able to retrieve disc, bulge, B/PS and, for some galaxies, even ring structural parameters fairly well, also helped to derive geometrical properties of the X-shape structures. For the three galaxies in our sample IC 2531, NGC 4013, and NGC 5529, the mean length of the rays in the X-shape structure were estimated, and they were found not to exceed 36% of the disc scale length. This proves that, while B/PS structures are not included in our modelling strategy, they should not have a large impact on the decomposition results for the stellar discs. Apart from that, a new 3D luminosity-density function based on the Ferrers ellipsoid has been implemented in the last version 1.5 of imfit. For future applications, this component can be used to accurately model 3D bar (B/PS) structures in galaxies.
Another concern can be the ‘thin disc+thick disc’ model which we used in our fitting. Studies by de Jong (2008) and Sandin (2014, 2015) demonstrated that the extended wings of the PSF can seriously affect our interpretation of different low surface brightness (LSB) features in galaxies, such as LSB haloes and possibly thick discs, which, in fact, may occur as artefacts produced by diffuse scattered light. In particular, Sandin (2014) showed for the HEROES galaxy NGC 5907 that its faint halo (or the thick disc in our study) can be interpreted as scattered light from a thin disc with a high surface brightness. Unfortunately, there was no possibility to create extended PSFs for the reference frames used in our study. Also, the PSF model should be at least as large as the observed galaxy, therefore PSF convolution in this case would strongly increase the time of the modelling, for both the imfit and fitskirt analysis. We should also note that NGC 5907 is the most spatially resolved galaxy in our sample, therefore the effects of the scattered light can be seen for it better than for smaller galaxies. Whether the obtained imfit models are physically valid or not, is hard to say without proper extended PSFs for all the reference images. However, recent work by Comerón et al. (2018) on a sample of 141 edge-on galaxies convincingly proves that the thick discs ‘are not an artefact caused by scattered light as has been suggested elsewhere’. They used a robust technique to fit a model consisting of two stellar discs to the SB profiles of these galaxies perpendicular to the galaxy plane. Also, they used a revised PSF for IRAC 3.6 with extended wings up to 2.5 arcmin. Unfortunately, this work appeared only recently, and it was impossible to use this new kernel earlier. Also, half of our galaxies are larger than their new PSF. Therefore we would not be able to perform a decomposition of all HEROES galaxies with the same level of accuracy. Anyway, since the thick disc is more extended in the vertical direction and faint in the dust disc plane (the average thick disc-to-thin disc scale height ratio for our sample is larger than 5.5), its influence on the dust emission in the infrared should be minimal.
Two galaxies in our sample (NGC 973 and NGC 4013) exhibit some AGN characteristic features. In our analysis, however, we do not include this component because of the following reasons. First, NGC 973 and NGC 4013 host low-luminosity AGNs (they contribute no more than several per cent to the total galaxy luminosities in the IRAC 3.6 m band, Salo et al. 2015). Second, the dusty torus emission has a peak between 10 and 30 m (unless a very large, non-physical torus is considered, see Fritz et al. 2006; Nenkova et al. 2008; Stalevski et al. 2012). After the peak, the emission decreases rapidly with a Rayleigh-Jeans curve. Third, the dust emission from the dust torus, which is heated by photons from the accretion disc, might influence the MIR SED, but in the FIR a very minor contribution is expected from this heating.
In this work, we applied the 3D broken exponential disc geometry to a small sample of edge-on galaxies, probably for the first time. For UGC 4277, the break radius was not robustly determined. Therefore for this galaxy we used a simple exponential disc. We claim that using the model of the broken exponential disc should give more reliable results for measuring the disc scale length and the broken radius since with this model we are able to describe the observed SB distribution in galaxy discs more accurately. For the other six HEROES galaxies, the average radius at which the break occurs is (or kpc), whereas Martín-Navarro et al. (2012) found it to be slightly larger: (or kpc). The discrepancy may be caused by the fact that some of the HEROES galaxies are quite large, of Sb type, whereas most of the galaxies which were considered by Martín-Navarro et al. (2012) are rather small, bulgeless galaxies. The nature of the disc breaks and truncations has been well studied (see e.g. de Grijs et al. 2001; Pohlen et al. 2004, 2007; Pohlen & Trujillo 2006; Florido et al. 2006; Comerón et al. 2012), and as proposed in Martín-Navarro et al. (2012), the first phenomenon may be caused by a threshold in the star formation (and, thus, change in the stellar population) at this radius, while the second one can be a real drop in the stellar mass density of the disc, which is associated with the maximum angular momentum of the stars or the protogalactic cloud (van der Kruit 1987, 1988).
5.2 On our RT models and the dust energy balance problem
In our approach, for all the galaxies we applied the same model of a double-exponential dust disc. However, the comparison of our skirt simulations with the observations in the MIR/FIR reveals that for NGC 4013 and likely NGC 5907 this model should be replaced by a dust ring. Since one of the aims of this study was to apply exactly the same approach to all HEROES galaxies using the same dust component, the fitting with a dust ring is beyond the scope of this paper. The same can be said about the spiral arms found in most HEROES galaxies.
In our oligochromatic modelling we used a fixed geometry of the stellar components, although early studies (Peletier et al. 1994; Corradi et al. 1996; Xilouris et al. 1999) and recent (Casasola et al. 2017) show that the disc scale lengths are wavelength dependent. However, according to Casasola et al. (2017), the -normalised scale-lengths in face-on galaxies are almost identical (0.26 for the band and 0.24 for IRAC 3.6). The dependence of the bulge parameters on wavelength should be stronger (Vulcani et al. 2014). However since only one galaxy in our sample, NGC 973, has a prominent bulge, we expect that our fixed geometrical model should be satisfactorily applicable at optical wavelengths as well.
A promising result of our oligochromatic fitting is that using additional NIR observations, apart from optical ones, can better constrain the output dust model. Thus, we found that for the 12 edge-on galaxies from De Geyter et al. (2014), for which the fitskirt model was fitted to the images, the uncertainties on the parameters are: %, %, %. Whereas, for the HEROES galaxies, we found %, %, %. Apparently, one of the reasons why the output oligochromatic models of the HEROES galaxies are better constrained is that we took into account a larger number of bands within a broader range of wavelengths. In addition, we should notice that the spatial resolution of the HEROES galaxies is much better than that of the 12 CALIFA (Sánchez et al. 2012) galaxies (by major axis, versus ). All in all, our fitskirt models adequately describe the observations we used in our fitting.
In our panchromatic simulations we used the same ages for the same stellar components for all sample galaxies. Obviously, this is not the case since the evolutionary status of each galaxy is different. We tried to vary the ages on a grid of values from 3 Gyr to 8 Gyr for the thin disc and from 5 Gyr to 12 Gyr for the thick disc and the bulge, to find for which ages our model SED better matches the observed SED in the optical and NIR (this is similar to the approach we applied to IC 2531 in M16). However, our attempt to find optimal ages was not successful because the recovered ages were always unrealistic (e.g. the age of the thin disc was larger than the age of the thick disc). The precise determination of the ages and metallicities of different stellar populations should be based on spectral observations of the galaxy, and, therefore, this is outside the scope of this paper. However, in our simulations we used different stellar ages from that grid of possible values and did not find a significant influence on the dust emission in the MIR-submm domain. Therefore we set the same ages of the stellar populations for all HEROES galaxies.
|Parameter||unit||NGC 973||UGC 4277||IC 2531||NGC 4013||NGC 4217||NGC 5529||NGC 5907|
In order to quantify the discrepancy of the dust emission between the observations and our panchromatic models, we applied a standard approach, in which the dust SED is well approximated by a simple single modified black-body (MBB) radiation (Hildebrand 1983). The details, along with the SED fits, are given in the Appendix E. The results of the MBB fitting are presented in Table 6. As one can see, the mass determined from our dust RT fitting (denoted as ) is typically 1.5–2 times smaller than the dust mass derived from the FIR emission. We confirm the result from Verstappen et al. (2013) that for NGC 973 both dust masses are in good agreement. For NGC 5907, the optically determined dust mass is significantly overestimated compared to its FIR mass.
The average dust temperature K and dust luminosity (obtained by integrating from 70 m to 1000 m) are typical values for the interstellar dust medium in spiral galaxies (see e.g. Davies et al. 2012; Galametz et al. 2012; Ciesla et al. 2014).
The ratio shows how well the dust FIR luminosity, retrieved from the MBB fitting, and the dust luminosity, computed from our panchromatic simulations (with a young stellar population included), are consistent with each other. We can see that for NGC 5907, the infrared luminosity from our panchromatic RT modelling and from the MBB fitting are very close. For NGC 973, UGC 4277, and IC 2531, our panchromatic models underestimate the infrared luminosity by a factor 1.3–1.6, whereas for the remaining galaxies we find an underestimation by a factor 3–5. If we carefully look at our panchromatic simulations and compare them with the real observations, we can see that such large discrepancies are found mainly in the galaxies, where a complex dust structure is observed. For example, the underestimation of the infrared luminosity emission in NGC 973 is 30 per cent, whereas in NGC 5529 it is 400%! The latter galaxy has very bright spiral arms which have been detected by Allaert et al. (2015) in its Hi maps. NGC 973, IC 2531, and NGC 4217 also show some prominent features in the MIR/FIR maps, which may be related to fainter spiral structures. NGC 4013 and NGC 5907 are likely to have some ring structures, which should explain their observed profiles in the MIR/FIR.
Another useful test to study the dust energy balance problem in the HEROES galaxies is by quantifying how much of the observed dust emission deficit is due to the model not capturing the observations, and how much is related to sub-pixel structures (clumps and dust emissivities). In Table 6, we also provide the following ratios which characterise these features. For each PACS 160m image (at the peak of the dust emission and where the spatial resolution is high enough), we measure the standard deviation of the image background and determine the outermost isophote of the galaxy at the level . The same level is used to determine the outermost isophote for the skirt simulation, at the same wavelength. For all galaxies in our sample, the latter lies within the former one (the real galaxies are more radially extended than their models, see Figs. 26–30). Thus, to measure the deficit of light in the galaxy outskirts (not captured by the model), we find the ratio between the skirt flux and the observed flux, both measured within an annulus between the skirt outermost isophote and the real one, denoted in Table 6 as . We also measure the ratio of the skirt and observed fluxes within the skirt isophote (inner isophote), to show how consistent our model and the observation are in the inner part of the galaxy. We can see that in the inner part, for the three galaxies NGC 4013, NGC 4217 and NGC 5529, the discrepancy between the model and the observation is 4–5 times; for these galaxies we also see the highest discrepancy between the model and observed FIR luminosities. This discrepancy comes from the small- and large-scale structure (bright dust clumps, spiral arms, bars, rings, etc.). For the outer part of the galaxy, we can see that the model severely underestimates dust emission, by a factor 5–10. However, the fraction of the flux between the inner- and outermost isophotes to the total galaxy flux (where ) varies from 3% (for NGC 5907) up to 20% (for UGC 4277), with an average %. Thus, due to its small fraction, the vertically and radially extended emission (with respect to the emission based on the optical model) cannot be the only reason of the observed dust energy balance problem. Nevertheless, the fact that in some galaxies this emission can contribute quite a different fraction to the total dust luminosity is interesting by itself. For NGC 891, Kamphuis et al. (2007); Bocchio et al. (2016) showed that, apart from the regular dust disc, there is a second spatial dust component which extends up to 1.9 kpc above the galactic mid-plane. Based on GALEX and Swift observations, diffuse UV haloes have recently been discovered around several galaxies (Hodges-Kluck & Bregman 2014; Shinn & Seon 2015; Baes & Viaene 2016). This points to the presence of a substantial amount of diffuse extra-planar dust, which scatters UV radiation emitted from the galactic disc.
An interesting result of our study, which confirms the findings of Verstappen et al. (2013) (see their fig. 12), is a weak trend between the difference in the values of the dust mass, as computed from our fitskirt radiative transfer modelling and as computed from our black-body fitting, and the dust scale height. To plot this correlation, which is presented in Fig. 6 (the regression coefficient ), we added the results for NGC 891 and NGC 4565 (as calculated in Verstappen et al. 2013 using the results from Bianchi & Xilouris 2011 and De Looze et al. 2012b, respectively) and the results from DG15 for IC 4225 and NGC 5166. Of course, the number of points is insufficient to make a strong conclusion, but the physical meaning of it may imply that the thinner the dust disc is, the clumpier its structure is. We suppose, that small clumpy regions do not contribute to the global optical extinction at all, and their presence can only be revealed from the thermal emission of the dust they contain, showing up in the FIR (see e.g. Popescu et al. 2000).
To summarise, in the optical, it is hard to discern any sign of a large or small-scale structure in the ISM, due to line-of-sight projection effects. Smooth, axisymmetric models generally describe the stellar emission and the dust attenuation very well. However, as shown by Saftly et al. (2015), fitting a complex, clumpy 3D structure with simple axisymmetric models can easily underestimate the dust content of a galaxy by a factor of a few. In the MIR/FIR, such small- and large-scale structure can be observable (if the resolution is high enough). Our observations and modelling support the idea that the inhomogeneous, clumpy ISM, along with various large-scale components, is responsible for the dust energy balance problem.
In the frame of the DustPedia project, special attention will be paid to the study of the dust mass, distribution and heating mechanisms in normal galaxies. A set of face-on galaxies will be considered and a new strategy will be applied to them. The unified model will consist of five components: an evolved stellar disc, a young stellar disc, an ionising stellar disc, a dust disc, and a stellar bulge. Each component will be defined in a special way using available multi-wavelength observations. The main difference of this approach from what we used in this study is that the geometry of each component will be constructed directly from the observed images, by means of the deprojection algorithm available in skirt. Using skirt simulations for the defined model, which has the three remaining fit parameters (the luminosity of the young and ionising stellar components and the total dust mass), the fitting of the simulated model SED to the observed global fluxes will be done. This will allow us to determine the fitted parameters in a self-consistent way. Some initial studies on M 51 (De Looze et al. 2014) and the Andromeda galaxy (Viaene et al. 2017) have been done, and the final strategy has been applied to M81 (Verstocken et al. in prep.) and NGC 1068 (Viaene et al. in prep). This will help to study the dust energy balance problem in a different way.
In an attempt to address the dust energy balance problem in Milky Way type galaxies, we have carried out a detailed analysis of the distribution of interstellar dust and stellar constituents in the remaining six HEROES galaxies (the prototype galaxy IC 2531 was described in M16, ) by computing oligochromatic radiative transfer models and panchromatic simulations. By doing so, we constrained both the stellar and dust distributions and studied the famous dust energy balance problem in the sample galaxies.
In our study, we successfully applied the following technique. First, for each galaxy, a photometric model, consisting of a thin disc, a thick disc, and a bulge, was found using the imfit code applied to the NIR/MIR ( band/IRAC 3.6m) images. These stellar models were used in subsequent fitskirt oligochromatic modelling with a fixed stellar geometry and a dust double-exponential disc. In this step, images in different optical and NIR wavebands (from seven to nine bands in total) were fitted simultaneously. The obtained models describe the fitted observations very well and are well-constrained (in comparison to the results from monochromatic or oligochromatic fitting performed in previous works).
The comparison of the results of our modelling with other works on RT fitting of edge-on galaxies shows that, in general, our results are consistent with earlier results.
Using the new THEMIS dust model, we performed skirt simulations with and without an additional young stellar component, which was added to match FUV/NUV observations. A uniform recipe of modelling was applied to all the galaxies. For five HEROES galaxies (UGC 4277, IC 2531, NGC 4013, NGC 4217, and NGC 5529), the constructed SEDs show the same dust energy balance problem even with an additional young stellar component – the dust emission in the FIR/submm is underestimated by a factor of 1.5 to 4. For the galaxy NGC 973, the emission in the FIR/submm is almost completely recovered in our modelling (albeit no young stellar population was added).
For only one galaxy (NGC 5907) in our sample, we found that the oligochromatic model, consisting of two exponential discs to describe the stellar thin and thick discs, another exponential disc for the dust distribution, and a Sérsic bulge, was not able to accurately reproduce the observations of this galaxy in the fitted wavebands which can be seen in the performed panchromatic simulations. In all other cases, based on the dust attenuation in the opical and NIR, we were able to construct oligochromatic models of the galaxies, which match the observations in these bands very well, and constrain their parameters to an acceptable accuracy. In all of the residual frames for our oligochromatic models more than half of the pixels show deviations of at most 20%. Parameters of the stellar components, as well as of the dust disc, are determined within %.
The comparison of the panchromatic simulations with the real ones reveals a more complex nature of the dust constituent in galaxies – it is often more extended in the radial and vertical directions than what can be inferred from our optical modelling. Some traces of spirals, rings, dust at high galactic latitudes, and signatures of a clumpy, inhomogeneous medium are clearly present. Thus, we suggest that these small- and large-scale structures may be responsible for the dust energy balance problem. We should notice, however, that the presence of different heating sources, which can be clearly distinguished at different wavelengths in the NIR–FIR domain (as shown in the Appendix D), suggests that the nature of the dust and stellar constituents and their link in galaxies needs a further investigation.
The average face-on optical depth in the band for the sample galaxies is (for comparison, LSB galaxies with diffuse clumpy dust discs are almost transparent: , MacLachlan et al. 2011). For three galaxies (NGC 973, NGC 4013, and NGC 5907), a significant part of the light is blocked even when the galaxy is seen face-on.
The results of the decomposition, as well as of the dust RT fitting (including uncertainties on individual measurements), of all structural parameters are listed in Tables 2 and 3. The imfit, fitskirt and skirt models (with the constructed SEDs) are presented in Appendices A, B, C, and D (available online).
Acknowledgements.A.V.M. is a beneficiary of a postdoctoral grant from the Belgian Federal Science Policy Office, and also expresses gratitude for the grant of the Russian Foundation for Basic Researches number 14-02-00810 and 14-22-03006-ofi. F.A., M.B., I.D.L. and S.V. gratefully acknowledge the support of the Flemish Fund for Scientific Research (FWO-Vlaanderen). M.B. acknowledges financial support from the Belgian Science Policy Office (BELSPO) through the PRODEX project âHerschel-PACS Guaranteed Time and Open Time Programs: Science Exploitationâ (C90370). T.M.H. acknowledges the CONICYT/ALMA funding Program in Astronomy/PCI Project No:31140020. He also acknowledges the support from the Chinese Academy of Sciences (CAS) and the National Commission for Scientific and Technological Research of Chile (CONICYT) through a CAS-CONICYT Joint Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. J.V. acknowledges support from the European Research Council under the European Unionâs Seventh Framework Programme (FP/2007-2013)/ ERC Grant Agreement nr. 291531. This work was supported by CHARM and DustPedia. CHARM (Contemporary physical challenges in Heliospheric and AstRophysical Models) is a phase VII Interuniversity Attraction Pole (IAP) program organized by BELSPO, the BELgian federal Science Policy Office. DustPedia is a collaborative focused research project supported by the European Union under the Seventh Framework Programme (2007-2013) call (proposal no. 606847). The participating institutions are: Cardiff University, UK; National Observatory of Athens, Greece; Ghent University, Belgium; UniversitÃ© Paris Sud, France; National Institute for Astrophysics, Italy and CEA (Paris), France. The Faulkes Telescopes are maintained and operated by Las Cumbres Observatory Global Telescope Network. We also thank Peter Hill and the staff and students of College Le Monteil ASAM (France), The Thomas Aveling School (Rochester, England), Glebe School (Bromley, England) and St Davidâs Catholic College (Cardiff, Wales). 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. This publication makes 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 publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This study makes use of observations made with the NASA Galaxy Evolution Explorer. GALEX is operated for NASA by the California Institute of Technology under NASA contract NAS5-98034. The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (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, UKSA (UK); and NASA (USA). We acknowledge the use of the ESA Planck Legacy Archive. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for nonHST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This research makes use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration, and the LEDA database (http://leda.univ-lyon1.fr).
- Allaert et al. (2015) Allaert, F., Gentile, G., Baes, M., et al. 2015, A&A, 582, A18
- Alton et al. (2004) Alton, P. B., Xilouris, E. M., Misiriotis, A., Dasyra, K. M., Dumke, M. 2004, A&A, 425, 109
- Aniano et al. (2011) Aniano, G., Draine, B. T., Gordon, K. D., Sandstrom, K. 2011, PASP, 123, 1218
- Athanassoula & Misiriotis (2002) Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35
- Baes et al. (2010) Baes, M., Fritz, J., Gadotti, D. A., et al. 2010, A&A, 518, L39
- Baes & Gentile (2011) Baes, M., & Gentile, G. 2011, A&A, 525, A136
- Baes & Van Hese (2011) Baes, M., & Van Hese, E. 2011, A&A, 534, A69
- Baes et al. (2011) Baes, M., Verstappen, J., De Looze, I., Fritz, J., Saftly, W., Vidal Pérez, E., Stalevski, M., Valcke, S. 2011, ApJS, 196, 22
- Baes & Viaene (2016) Baes, M., & Viaene, S. 2016, A&A, 587, A86
- Benn & Ellison (1998) Benn, C. R., & Ellison, S. L. 1998, New A Rev., 42, 503
- Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
- Bianchi (2007) Bianchi, S. 2007, A&A, 471, 765
- Bianchi (2008) Bianchi, S. 2008, A&A, 490, 461
- Bianchi & Xilouris (2011) Bianchi, S., & Xilouris, E. M. 2011, A&A, 531, L11
- Bianchi et al. (2014) Bianchi, L., Conti, A., & Shiao, B. 2014, Advances in Space Research, 53, 900
- Bocchio et al. (2016) Bocchio, M., Bianchi, S., Hunt, L. K., & Schneider, R. 2016, A&A, 586, A8
- Bottema et al. (1987) Bottema, R., Shostak, G. S., & van der Kruit, P. C. 1987, Nature, 328, 401
- Bruzual & Charlot (2003) Bruzual, G., Charlot, S. 2003, MNRAS, 344, 1000
- Bureau & Freeman (1999) Bureau, M., & Freeman, K. C. 1999, AJ, 118, 126
- Bureau et al. (2006) Bureau, M., Aronica, G., Athanassoula, E., et al. 2006, MNRAS, 370, 753
- Camps & Baes (2015) Camps, P., Baes, M. 2015, Astronomy and Computing, 9, 20
- Camps et al. (2015) Camps, P., Misselt, K., Bianchi, S., et al. 2015, A&A, 580, A87
- Caon, Capaccioli & D’Onofrio (1993) Caon, N., Capaccioli, M., D’Onofrio, M. 1993, MNRAS, 265, 1013
- Casasola et al. (2017) Casasola, V., Cassarà, L. P., Bianchi, S., et al. 2017, A&A, 605, A18
- Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
- Chung & Bureau (2004) Chung, A., & Bureau, M. 2004, AJ, 127, 3192
- Ciesla et al. (2014) Ciesla, L., Boquien, M., Boselli, A., et al. 2014, A&A, 565, A128
- Ciotti & Bertin (1999) Ciotti, L. & Bertin, G. 1999, A&A, 352, 447
- Clark et al. (2018) Clark, C. J. R., Verstocken, S., Bianchi, S., et al. 2018, A&A, 609, A37
- Combes & Sanders (1981) Combes, F., & Sanders, R. H. 1981, A&A, 96, 164
- Combes et al. (1990) Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82
- Comerón et al. (2011) Comerón, S., Elmegreen, B. G., Knapen, J. H., et al. 2011, ApJ, 738, L17
- Comerón et al. (2012) Comerón, S., Elmegreen, B. G., Salo, H., et al. 2012, ApJ, 759, 98
- Comerón et al. (2018) Comerón, S., Salo, H., & Knapen, J. H. 2018, A&A, 610, A5
- Corradi et al. (1996) Corradi, R. L. M., Beckman, J. E., del Rio, M. S., di Bartolomeo, A., & Simonneau, E. 1996, New Extragalactic Perspectives in the New South Africa, 209, 523
- Dalcanton & Bernstein (2002) Dalcanton, J. J., & Bernstein, R. A. 2002, AJ, 124, 1328
- Dale et al. (2012) Dale, D. A., Aniano, G., Engelbracht, C. W., et al. 2012, ApJ, 745, 95
- Dasyra et al. (2005) Dasyra, K. M., Xilouris, E. M., Misiriotis, A., Kylafis, N. D. 2005, A&A, 437, 447
- Davies et al. (2012) Davies, J. I., Bianchi, S., Cortese, L., et al. 2012, MNRAS, 419, 3505
- Davies et al. (2017) Davies, J. I., Baes, M., Bianchi, S., et al. 2017, PASP, 129, 044102
- De Geyter et al. (2013) De Geyter, G., Baes, M., Fritz, J., Camps, P. 2013, A&A, 550, A74
- De Geyter et al. (2014) De Geyter, G., Baes, M., Camps, P., Fritz, J., De Looze, I., Hughes, T. M., Viaene, S., Gentile, G. 2014, MNRAS, 441, 869
- De Geyter et al. (2015a) De Geyter, G., Baes, M., De Looze, I., et al. 2015a, MNRAS, 451, 1728 (DG15)
- De Geyter (2015) De Geyter, G. 2015b, Mapping the interstellar dust in galaxies through radiative transfer modelling. Ghent University. Faculty of Sciences, Ghent, Belgium
- de Grijs et al. (2001) de Grijs, R., Kregel, M., & Wesson, K. H. 2001, MNRAS, 324, 1074
- De Looze et al. (2012a) De Looze, I., Baes, M., Fritz, J., Verstappen, J. 2012a, MNRAS, 419, 895
- De Looze et al. (2012b) De Looze, I., Baes, M., Bendo, G. J., et al. 2012b, MNRAS, 427, 2797
- De Looze et al. (2014) De Looze, I., Fritz, J., Baes, M., et al. 2014, A&A, 571, A69
- de Jong (2008) de Jong, R. S. 2008, MNRAS, 388, 1521
- Dowell et al. (2014) Dowell, C. D., Conley, A., Glenn, J., et al. 2014, ApJ, 780, 75
- Eisenstein et al. (2011) Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72
- Elmegreen & Elmegreen (1984) Elmegreen, D. M., & Elmegreen, B. G. 1984, ApJS, 54, 127
- Erwin et al. (2005) Erwin, P., Beckman, J. E., & Pohlen, M. 2005, ApJL, 626, L81
- Erwin et al. (2008) Erwin, P., Pohlen, M., & Beckman, J. E. 2008, AJ, 135, 20
- Erwin (2015) Erwin, P. 2015, ApJ, 799, 226
- Eskew et al. (2012) Eskew, M., Zaritsky, D., & Meidt, S. 2012, AJ, 143, 139
- Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
- Ferrers (1877) Ferrers, N. M., 1877, Quart. J. Pure Appl. Math., 14, 1
- Fisher & Drory (2010) Fisher, D. B., & Drory, N. 2010, ApJ, 716, 942
- Florido et al. (1991) Florido, E., Prieto, M., Battaner, E., et al. 1991, A&A, 242, 301
- Florido et al. (2006) Florido, E., Battaner, E., Guijarro, A., Garzón, F., & Castillo-Morales, A. 2006, A&A, 455, 467
- Fritz et al. (2006) Fritz, J., Franceschini, A., & Hatziminaoglou, E. 2006, MNRAS, 366, 767
- Gadotti (2009) Gadotti, D. A. 2009, MNRAS, 393, 1531
- Galametz et al. (2012) Galametz, M., Kennicutt, R. C., Albrecht, M., et al. 2012, MNRAS, 425, 763
- Goldberg (1989) Goldberg, D. E. 1989, Reading: Addison-Wesley, 1989
- Griffin et al. (2010) Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3
- Grosbol (1993) Grosbol, P. 1993, PASP, 105, 651
- Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
- Ho et al. (1997) Ho, L. C., Filippenko, A. V., & Sargent, W. L. W. 1997, ApJS, 112, 315
- Hodges-Kluck & Bregman (2014) Hodges-Kluck, E., & Bregman, J. N. 2014, ApJ, 789, 131
- Holland et al. (1999) Holland, W. S., Robson, E. I., Gear, W. K., et al. 1999, MNRAS, 303, 659
- Holwerda et al. (2012) Holwerda, B. W., Bianchi, S., Böker, T., et al. 2012, A&A, 541, L5
- Irwin et al. (2007) Irwin, J. A., Kennedy, H., Parkin, T., & Madden, S. 2007, A&A, 474, 461
- Jones et al. (2017) Jones, A. P., Koehler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, 46
- Kamphuis et al. (2007) Kamphuis, P., Holwerda, B. W., Allen, R. J., Peletier, R. F., & van der Kruit, P. C. 2007, A&A, 471, L1
- Karachentseva (1973) Karachentseva, V. E. 1973, Astrofizicheskie Issledovaniia Izvestiya Spetsial’noj Astrofizicheskoj Observatorii, 8, 3
- Kylafis & Bahcall (1987) Kylafis, N. D., & Bahcall, J. N. 1987, ApJ, 317, 637
- Köhler et al. (2014) Köhler, M., Jones, A., & Ysard, N. 2014, A&A, 565, L9
- Köhler et al. (2015) Köhler, M., Ysard, N., & Jones, A. P. 2015, A&A, 579, A15
- Kregel et al. (2004) Kregel, M., van der Kruit, P. C., & de Blok, W. J. G. 2004, MNRAS, 352, 768
- Laurikainen & Salo (2016) Laurikainen, E., & Salo, H. 2016, Galactic Bulges, 418, 77
- Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3
- Lisenfeld et al. (2007) Lisenfeld, U., Verdes-Montenegro, L., Sulentic, J., et al. 2007, A&A, 462, 507
- Lütticke et al. (2000) Lütticke, R., Dettmar, R.-J., & Pohlen, M. 2000, A&AS, 145, 405
- MacLachlan et al. (2011) MacLachlan, J. M., Matthews, L. D., Wood, K., & Gallagher, J. S. 2011, ApJ, 741, 6
- Mahtessian (1998) Mahtessian, A. P. 1998, Astrofizika, 41, 308
- Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1
- Martínez-Delgado et al. (2008) Martínez-Delgado, D., Peñarrubia, J., Gabany, R. J., et al. 2008, ApJ, 689, 184-193
- Martínez-Delgado et al. (2009) Martínez-Delgado, D., Pohlen, M., Gabany, R. J., et al. 2009, ApJ, 692, 955
- Martín-Navarro et al. (2012) Martín-Navarro, I., Bakos, J., Trujillo, I., et al. 2012, MNRAS, 427, 1102
- Mathewson et al. (1992) Mathewson, D. S., Ford, V. L., & Buchhorn, M. 1992, ApJS, 81, 413
- Mazure & Capelato (2002) Mazure, A., & Capelato, H. V. 2002, A&A, 383, 384
- McGaugh et al. (2000) McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99
- Meidt et al. (2012) Meidt, S. E., Schinnerer, E., Knapen, J. H., et al. 2012, ApJ, 744, 17
- Meidt et al. (2014) Meidt, S. E., Schinnerer, E., van de Ven, G., et al. 2014, ApJ, 788, 144
- Méndez-Abreu et al. (2008) Méndez-Abreu, J., Corsini, E. M., Debattista, V. P., et al. 2008, ApJ, 679, L73
- Miller & Veilleux (2003) Miller, S. T., & Veilleux, S. 2003, ApJS, 148, 383
- Misiriotis et al. (2001) Misiriotis, A., Popescu, C. C., Tuffs, R., & Kylafis, N. D. 2001, A&A, 372, 775
- Moffat (1969) Moffat, A. F. J. 1969, A&A, 3, 455
- Morrison et al. (1997) Morrison, H. L., Miller, E. D., Harding, P., Stinebring, D. R., & Boroson, T. A. 1997, AJ, 113, 2061
- Mosenkov et al. (2016) Mosenkov, A. V., Allaert, F., Baes, M., et al. 2016, A&A, 592, A71
- Moshir et al. (1990) Moshir, M., Kopan, G., Conrow, T., et al. 1990, BAAS, 22, 1325
- Muñoz-Mateos et al. (2013) Muñoz-Mateos, J. C., Sheth, K., Gil de Paz, A., Meidt, S., Athanassoula, E., Bosma, A., Comerón, S., Elmegreen, D. M., & Elmegreen, B. G. 2013, ApJ, 771, 59
- Muñoz-Mateos et al. (2015) Muñoz-Mateos, J. C., Sheth, K., Regan, M., et al. 2015, ApJS, 219, 3
- Murakami et al. (2007) Murakami, H., Baba, H., Barthel, P., et al. 2007, PASJ, 59, S369
- Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160-180
- Patsis & Xilouris (2006) Patsis, P. A., & Xilouris, E. M. 2006, MNRAS, 366, 1121
- Peletier et al. (1994) Peletier, R. F., Valentijn, E. A., Moorwood, A. F. M., & Freudling, W. 1994, A&AS, 108, 621
- Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097
- Peterson (1979) Peterson, S. D. 1979, ApJS, 40, 527
- Pfenniger & Friedli (1991) Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75
- Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1
- Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A1
- Poglitsch et al. (2010) Poglitsch, A. et al. 2010, A&A, 518, L2
- Pohlen et al. (2004) Pohlen, M., Beckman, J. E., Hüttemeister, S., et al. 2004, Penetrating Bars Through Masks of Cosmic Dust, 319, 713
- Pohlen & Trujillo (2006) Pohlen, M., & Trujillo, I. 2006, A&A, 454, 759
- Pohlen et al. (2007) Pohlen, M., Zaroubi, S., Peletier, R. F., & Dettmar, R.-J. 2007, MNRAS, 378, 594
- Popescu et al. (2000) Popescu, C. C., Misiriotis, A., Kylafis, N. D., Tuffs, R. J., & Fischera, J. 2000, A&A, 362, 138
- Popescu & Tuffs (2002) Popescu, C. C., & Tuffs, R. J. 2002, MNRAS, 335, L41
- Popescu et al. (2011) Popescu, C. C., Tuffs, R. J., Dopita, M. A., et al. 2011, A&A, 527, A109
- Querejeta et al. (2015) Querejeta, M., Meidt, S. E., Schinnerer, E., et al. 2015, ApJS, 219, 5
- Regan (2013) Regan, M. W. 2013, American Astronomical Society Meeting Abstracts #221, 221, 230.02
- Riechers et al. (2013) Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329
- Rieke et al. (2004) Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25
- Rix & Zaritsky (1995) Rix, H.-W., & Zaritsky, D. 1995, ApJ, 447, 82
- Roussel (2013) Roussel, H. 2013, PASP, 125, 1126
- Sackett et al. (1994) Sackett, P. D., Morrisoni, H. L., Harding, P., & Boroson, T. A. 1994, Nature, 370, 441
- Saftly et al. (2015) Saftly, W., Baes, M., De Geyter, G., et al. 2015, A&A, 576, A31
- Salo et al. (2015) Salo, H., Laurikainen, E., Laine, J., et al. 2015, ApJS, 219, 4
- Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
- Sánchez et al. (2012) Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al. 2012, A&A, 538, A8
- Sancisi (1976) Sancisi, R. 1976, A&A, 53, 159
- Sanders et al. (2003) Sanders, D. B., Mazzarella, J. M., Kim, D.-C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607
- Sandin (2014) Sandin, C. 2014, A&A, 567, A97
- Sandin (2015) Sandin, C. 2015, A&A, 577, A106
- Savchenko et al. (2017) Savchenko, S. S., Sotnikova, N. Y., Mosenkov, A. V., Reshetnikov, V. P., & Bizyaev, D. V. 2017, MNRAS, 471, 3261
- Sazonov et al. (2007) Sazonov, S., Revnivtsev, M., Krivonos, R., Churazov, E., & Sunyaev, R. 2007, A&A, 462, 57
- Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
- Sérsic (1968) Sérsic, J.L. 1968, Atlas de Galaxias Australes, Observatorio Astronomico, Cordoba
- Seth et al. (2005) Seth, A. C., Dalcanton, J. J., & de Jong, R. S. 2005, AJ, 130, 1574
- Shang et al. (1998) Shang, Z., Zheng, Z., Brinks, E., et al. 1998, ApJ, 504, L23
- Sheth et al. (2010) Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397
- Shinn & Seon (2015) Shinn, J.-H., & Seon, K.-I. 2015, ApJ, 815, 133
- Smith et al. (2012) Smith, M. W. L., Gomez, H. L., Eales, S. A., et al. 2012, ApJ, 748, 123
- Stalevski et al. (2012) Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popović, L. Č. 2012, MNRAS, 420, 2756
- Stevens et al. (2005) Stevens, J. A., Amure, M., & Gear, W. K. 2005, MNRAS, 357, 361
- Stickel et al. (2004) Stickel, M., Lemke, D., Klaas, U., Krause, O., & Egner, S. 2004, A&A, 422, 39
- Tikhonov & Galazutdinova (2005) Tikhonov, N. A., & Galazutdinova, O. A. 2005, Astrophysics, 48, 221
- van der Kruit & Searle (1981) van der Kruit, P. C., & Searle, L. 1981, A&A, 95, 116
- van der Kruit (1987) van der Kruit, P. C. 1987, A&A, 173, 59
- van der Kruit (1988) van der Kruit, P. C. 1988, A&A, 192, 117
- Verstappen et al. (2013) Verstappen, J., Fritz, J., Baes, M., et al. 2013, A&A, 556, A54
- Viaene et al. (2015) Viaene, S., De Geyter, G., Baes, M., et al. 2015, A&A, 579, A103
- Viaene et al. (2016) Viaene, S., Baes, M., Bendo, G., et al. 2016, A&A, 586, A13
- Viaene et al. (2017) Viaene, S., Baes, M., Tamm, A., et al. 2017, A&A, 599, A64
- Vulcani et al. (2014) Vulcani, B., Bamford, S. P., Häußler, B., et al. 2014, MNRAS, 441, 1340
- Weedman et al. (2012) Weedman, D., Sargsyan, L., Lebouteiller, V., Houck, J., & Barry, D. 2012, ApJ, 761, 184
- Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296
- Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1
- Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868-1881
- Xilouris et al. (1997) Xilouris, E. M., Kylafis, N. D., Papamastorakis, J., Paleologou, E. V., & Haerendel, G. 1997, A&A, 325, 135
- Xilouris et al. (1999) Xilouris, E. M., Byun, Y. I., Kylafis, N. D., Paleologou, E. V., Papamastorakis, J. 1999, A&A, 344, 868 (X99)
- Zschaechner & Rand (2015) Zschaechner, L. K., & Rand, R. J. 2015, ApJ, 808, 153
- Zubko et al. (2004) Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211
Appendix A IMFIT models
Appendix B fitskirt models
Appendix C SEDs
Appendix D Panchromatic simulations
Appendix E Global SED fitting
We determined the dust masses for the galaxies in our sample by fitting a simple modified black-body model to the PACS and SPIRE data. The monochromatic luminosity, thus, can be written as
where is the dust mass, is the Planck function, is the dust temperature, is the dust emissivity at , and signifies a power-law dust emissivity in the FIR/sub-mm wavelength range. In our model, we assume an emissivity of m kg at m and adapted to the THEMIS model. The fits were done by using the idl procedure mpcurvefit.