© ESO 2018

1. Introduction

Cosmic dust severely obscures astronomical objects at ultraviolet (UV) and optical wavelengths and, thus, impedes our study of these objects in this wavelength range. However, with the beginning of the era of 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-submillimetre (FIR-submm) wavelengths.

Until recently, the poor resolution and limited wavelength coverage of IR instruments posed problems for observations of 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–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 greater detail than ever possible before.

Panchromatic radiative transfer (RT) 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, 2011; Bianchi 2008; Baes et al. 2010; 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 RT code to model the propagation of stellar light and its interaction with dust particles in a galaxy. In turn, the dust emission predicted from the RT 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 requires that dust features must easily be identified from optical as well as IR 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). Therefore, using such a multiwavelength technique, a complex three-dimensional (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 RT models and the observed emission in those wavebands (e.g. Popescu et al. 2000, 2011; Misiriotis et al. 2001; Alton et al. 2004; Dasyra et al. 2005; Baes et al. 2010; De Looze et al. 2012a; De Geyter et al. 2015; Mosenkov et al. 2016). Although RT 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 RT 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 with 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. (2015; 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 oligochromatic¹ model (using optical data), and adding a young stellar disc to match the 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 FIR emission still underestimates the observed fluxes by a factor of three. The proposed reasons for this discrepancy in IC 4225 might include the fact 1) that it is too small to properly fit the dust disc parameters, and 2) that 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, 2014, DG15). It is likely that there is no single origin for the dust energy balance problem, and several mechanisms could 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 D > 4 arcmin), and have Herschel observations in five PACS (Poglitsch 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 RT code based on their optical and nearinfrared (NIR) data (Xilouris et al. 1997, 1999; Bianchi 2007), therefore an indirect comparison with the previous models can be done. However, one should note 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 NIR 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 NIR imaging is used, in addition to optical data. Based on our oligochromatic fit, we constructed a panchromatic² RT 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 NIR observations to constrain their dust models and perform panchromatic simulations. We investigate the degree to which 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 of how common the dust energy balance problem might be in galaxies and will provide clues on a possible reason of its existence.

This paper is organised as follows. Section 2 reviews the sample as well as observation 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 Sect. 4, we present results of our study for each individual galaxy and compare our results with the literature. In Sect. 5, the results of our RT modelling procedure are discussed in the frame of the dust energy balance problem. The main conclusions of our study are summarised in Sect. 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 smb/sma – 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 of two to three. 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), whilst in others they are relatively thin and warped (UGC 4277, IC 2531).

All HEROES galaxies are group members (Peterson 1979; Mathewson et al. 1992; Mahtessian 1998), except for UGC 4277 which is classified as an isolated galaxy (Karachentseva 1973).

For the HEROES galaxies RT 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 Xilouris et al. (1999; hereafter 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.

2.1. Observations

To investigate the dust energy balance problem in a galaxy, we need to collect data over a large wavelength range – from farultraviolet (FUV) to FIR/submm regions. For this purpose, we used available databases of optical and NIR sky surveys to extract images of the HEROES galaxies. Here, one should note 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 refer to them as the reference images) simultaneously. The second set of observations was prepared to supplement galaxy SEDs with fluxes at UV and MIR-FIR wavelengths with λ > 4 µ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 ugriz 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 Auxiliaryport CAMera (ACAM, Benn & Ellison 1998). The observing run yielded a total of 9 h, during 3 moonless nights from 25 to 27 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 ugriz 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 B-, V- and R-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.58 m telescope of the Roque de Los Muchachos Observatory on La Palma. The observations in the K 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 K band and for the remaining galaxies we used 2MASS observations in all three JHK_s 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 Archive³. 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, 2015; Regan 2013; Salo et al. 2015)⁴. For UGC 4277, we used the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) in the W1 (3.4 µm) 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 µm (FUV) and 0.227 µm (NUV), which were downloaded from the All Sky Imaging Survey reachable trough the GalexView service⁵. WISE observations are available for the entire sky; therefore we make use of them in our study as well, in all four passbands, W1 (3.4 µm), W2 (4.6 µm), W3 (12 µm) and W4 (22 µm), as an important source of information regarding 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 Archive⁶ 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 were initially reduced and prepared for the subsequent analysis, as 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 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), and 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 (PTS⁷, Camps & Baes 2015; Verstocken et al. in prep.), which contains many functionalities and useful routines. The sky background was fitted by a two-dimensional (2D) 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 were 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). These were subsequently fitted with a Moffat (Moffat 1969; for the SDSS, TNG and WHT observations) or Gaussian (for the 2MASS and WISE W1 images) function. For the IRAC images, we used in-flight point response function (PRF) images⁸ 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 v13⁹ 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 2σ 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 Table 1 from Clark et al. 2018, while 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 2003; Lisenfeld et al. 2007) flux densities at 25, 60 and 100 µm taken from NED, and Planck (Planck Collaboration I 2014) APERFLUX values at 353, 545 and 857 GHz (which correspond to 850 µm, 550 µm and 350 µm, respectively) from the Planck Legacy Archive¹⁰. 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 µm and 850 µm 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 µm fluxes is very poor (to estimate the fluxes they used PSF photometry with a FWHM of about 60 arcsec at 160 µm, 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.

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, for example, a Sérsic (Sérsic 1968) bulge and thin and thick double-exponential discs, to a NIR image. IRAC 3.6 µm 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 K-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.6 µm and K 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 FITSKIRT¹¹ code (De Geyter et al. 2013, 2014). This allows one to find structural properties of the dust component by means of SKIRT (Baes et al. 2011; Camps & Baes 2015) oligochromatic RT 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 the geometrical sizes of the stellar components fixed to the values retrieved in the first step, but allow their relative luminosities to change. 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 RT 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 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. Firstly, in the present version of FITSKIRT, all wave-lengths 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 so this assumption is no longer valid. Secondly, the computational cost of such multiwavelength fitting is very high because not only does the number of free model parameters substantially increase, but so does the amount of data, which should be analysed simultaneously during the fitting.

3.1. Decomposition of the IRAC 3.6 µm images

The SB profile of a stellar bulge is usually fitted with a Sérsic model:(1)

where I _e,b is the effective SB, that is, the SB at the half-light radius of the bulge r _e,b, and b _n is a function of the Sérsic index n _b (Caon et al. 1993; Ciotti & Bertin 1999). An additional free parameter is the apparent axis ratio of the bulge q _b. We also apply a Sérsic profile for describing the SB of a halo for one HEROES galaxy, NGC 4013.

The 3D axisymmetric stellar disc can be described by a double-exponential function j(R, z), which in a cylindrical co-ordinate system (R, z) aligned with the disc (where the disc midplane has z = 0) is given by(2)

where j ₀ is the central luminosity density of the disc, h _R is the disc scale length, and h _z 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):(3)

In this formula, s parametrises the sharpness of the transition between the inner and outer profiles with the break radius R _b. Large values of the sharpness parameter (s ≫ 1) correspond to a sharp transition, and small values (s ∼ 1) set a very gradual break. The dimensionless quantity S is a scaling factor, given by(4)

In total, this model contains five free parameters: the scale length of the inner disc h _R,inn, the scale length of the outer disc h _R,out, the scale height h _z, the break radius R _b, and the sharpness of the break s.

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 in their radial profiles, therefore for them we used the model (3); 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 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 n = 100 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 RT modelling of the HEROES galaxies with a simple “disc + bulge” stellar model and a dust disc (described in De Geyter 2015). As shown in Table 4, 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 s = 5, 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.

Table 3 summarises the results of the IMFIT decomposition for all HEROES galaxies. The detailed description of the obtained IMFIT models is given in Sect. 4.

3.2. FITSKIRT fitting

FITSKIRT is a fitting code that combines the output of FITSKIRT 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 real edge-on galaxies (see De Looze et al. 2012b; De Geyter et al. 2014; Saftly et al. 2015; Viaene et al. 2015, M16). 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 the dust ring.

De Geyter et al. (2014) showed that oligochromatic fitting, that is, fitting applied to a number of bands simultaneously, has clear advantages over monochromatic fitting in terms of accuracy. In particular, the parameters that 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. X99; Baes et al. 2010; Verstappen et al. 2013; De Geyter et al. 2013, 2014):(5)

where M _d is the total dust mass, h _R,d is the radial scale length, and h _z,d is the vertical scale height of the dust disc. To calculate the central face-on optical depth, one should use the following expression:(6)

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 De Sect. 3.1. We allow only the dust disc parameters, the alignment of the galaxy centre, the inclination angle i, and the galaxy position angle PA 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 (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 Super-computer Center (Vlaams Supercomputer Centrum). For each galaxy, we used about 50 processor units when running FITSKIRT for one galaxy with the following specifications: 5 × 10⁵ photon packages, 200 individuals in a population, and 100 generations (these have been tested and found 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 h (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 h of work spread among 50 cores.

The results of the fits for all seven HEROES galaxies can be found in Table 4.

3.3. FITSKIRT 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 NIR 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-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 (Z = 0.02). 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 age of the stellar halo for NGC 4013 was assumed to be 10 Gyr.

As follows from Popescu et al. (2000), De Looze et al. (2012a,b), DG15 and M16, an additional source of UV luminosity, that is, 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 one third 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 and 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. Results

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 or 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 (L _t/L _T = 0.13). The bulge/bar component fitted with the Sérsic function is flat and has an extended SB profile (the Sérsic index is large, n _b = 5.48) with an ordinary effective radius. The galaxy shows the largest fraction of the bulge/bar to the total galaxy luminosity among all sample galaxies (L _b/L _tot = 0.56).

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 for it are very poor, and therefore we did not include NUV and FUV fluxes in the global SED; also, we did not, therefore, 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% of the observed flux. However, this is not as much as for the other galaxies, as we see below. The panchromatic simulations and observations look highly similar in the optical and NIR bands (see Fig. 5). The WISE W3, 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 sma = 36.95 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 H I disc demonstrate a small warp as well, with the position angle deviating by about 4 degrees at most from its central value. UGC 4277 has a total atomic hydrogen mass of 2.03 × 10¹⁰ M _⊙ and, similarly to NGC 973, the H I 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 (Fig. C.1).

For the IMFIT modelling (which was based on the K-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. A.1). 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. B.1). 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 g, r and i bands. This is evidence that UGC 4277 has a bar which is not so extended as in the other HEROES galaxies or that it is not viewed purely side-on (Bureau et al. 2006) (Fig. D.1).

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. C.2) 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 200 µm, by a factor of about two. Consequently, the panchromatic simulations cannot entirely describe the observations. The panchromatic simulations in the optical and MIR are compatible with the observations (see Fig. D.3). However the comparison between the WISE W3 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 (2009) found a giant, arcing stellar structure of low SB 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 H I 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; Fig. A.2).

Since this galaxy has a bright stellar halo, we used a four-component model: “thin disc + thick disc + bulge + halo” (see Fig. A.3). 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 (r _e,h ≈ 7.4 kpc), the Sérsic index is close to 1, therefore it can be called the second thick disc (Fig. B.2).

Since our FITSKIRT model does not take into account the X-shape structure, this feature seems a prominent detail in the residuals in Fig. B.3. 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 u band image is obviously too noisy and, therefore, the comparison in this band is meaningless).

Our panchromatic simulations (see Figs. C.3 and D.4) 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 u band model is closer to the corresponding observation. All other optical –4.6 µm observations are described exceptionally well by our model. The MIR-FIR observations reveal a probable ring structure (since the inclination angle is almost exactly 90°, 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 2.50 × 10⁹ M _⊙. 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 3.81 × 10⁸ M _⊙. 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. A.4). As one can see from Table 3, 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-height ratio: 3.8 – apparently, this galaxy looks thicker than the others.

The FITSKIRT model is quite adequate with regard to the observations (see Fig. B.4). Only some traces of the dust attenuation are left in the residuals in the r, i, and z 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. C.4). The typical factor of flux underestimation is four to five. However, the comparison of the panchromatic simulations and real observations for this galaxy is very good (see Fig. D.4), even for the FUV and NUV bands. The WISE W3 and W4 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 90° inclination (87.9°). 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 H I bridges. They also found a bright H I 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 2.69 × 10¹⁰ M _⊙.

Our IMFIT model consists of a thin broken disc, a thick disc, and a bulge (see Fig. A.5). 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 two – 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. B.5), with some traces of the dust content in the residuals in the r and i bands. However, the SKIRT simulations (see Fig. C.5) reveal a systematic underestimation of flux densities at all wavelengths larger than 10 µm by a factor of 4–5, for the models both with and without the young disc. The comparison of the panchromatic simulations with the real galaxy images (see Fig. D.6) indicates that starting from the WISE W3 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 wave-lengths) and apparent existence of spirals (which are distinct emission features on both sides from the centre, from WISE W3 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 (84.9°). 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 1.99 × 10¹⁰ M _⊙. 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 H I 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. A.6). 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. B.6. 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 easily visible. However, the overall model looks good and compatible with the observations.

Nevertheless, if we look at the parameters listed in Table 4, the FITSKIRT model for NGC 5907 differs from the others. From Table 4 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. C.6, 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. D.6). 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 W2 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 S⁴G results (Salo et al. 2015) for four HEROES galaxies: NGC 4013, NGC 4217, NGC 5529, and NGC 5907. However, some caution should be taken. Firstly, our decomposition is performed for 3D stellar discs, whereas the GALFIT code (Peng et al. 2010), which was used in the S⁴G project, works only with projected intensities. Secondly, for these galaxies, the S⁴G 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.6 µm the Gaussian PSF FWHM is 2.1″). 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 z ₀, 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 z ≫ h _z: z ₀ = 2 h _z. The results of ours and the S⁴G decomposition are presented together in Table 5. 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 RT models have been obtained before for all HEROES galaxies, it is interesting to compare the results of this work with the literature. Table 6 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 used in the analysis is limited – from two to five (optical and NIR). 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 i 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 B/T. In addition, it is particularly hard to estimate the uncertainties on these parameters.

Another reason for the differences with previous results must also be 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.

5. Discussion

5.1. 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, therefore, 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 from 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 m = 4 Fourier mode has been applied to describe B/PS structures in 22 edge-on galaxies (Savchenko et al. 2017). A multicomponent decomposition, with 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 SB (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 SB. 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 analyses. 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 IR 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 for the following reasons. Firstly, 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). Secondly, 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. Lastly, 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 0.64 ± 0.08 R ₂₅ (or 18.01 ± 5.85 kpc), whereas Martín-Navarro et al. (2012) found it to be slightly larger: 0.77 ± 0.06 R ₂₅ (or 8 ± 1 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. 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 (Peletier et al. 1994; Corradi et al. 1996, X99) and recent (Casasola et al. 2017) studies show that the disc scale lengths are wavelength dependent. However, according to Casasola et al. (2017), the R ₂₅- normalised scale-lengths, h _R/R ₂₅, in face-on galaxies are almost identical (0.26 for the g 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 griz images, the uncertainties on the parameters are: σ(M _d) = 20.7 ± 12.1%, σ(h _R,d) = 18.7 ± 9.2%, σ(h _z,d) = 15.6 ± 8.6%. Whereas, for the HEROES galaxies, we found σ(M _d) = 9.9 ± 3.7%, σ(h _R,d) = 7.0 ± 3.2%, σ(h _z,d) = 11.3 ± 12.8%. 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, one should note 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, 364″ ± 163″ versus 116″ ± 34″). 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 the ages for which our model SED best 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 beyond 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.

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 Appendix E. The results of the MBB fitting are presented in Table 7. 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 〈T _d〉 = (20.6 ± 1.4) K and dust luminosity 〈L _IR〉 = (14.75 ± 5.09) 10⁹ L _⊙ (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 the extent to which 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 IR 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 IR luminosity by a factor of 1.3–1.6, whereas for the remaining galaxies we find an underestimation by a factor of 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 IR luminosity emission in NGC 973 is 30%, 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 7, we also provide the following ratios which characterise these features. For each PACS 160 µm 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 2σ. The same 2σ 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 (the real galaxies are more radially extended than their models, see Figs. D.1–D.6). Therefore, 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 7 as . We also measure the ratio of the SKIRT and observed fluxes within the SKIRT isophote (inner isophote) to show the consistency between our model and observation 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 difference between the model and the observation is a factor of four to five; 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 structures (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 of five to ten. However, the fraction of the flux between the inner- and outermost isophotes to the total galaxy flux F _out/F _full (where F _full = F _in + F _out) varies from 3% (for NGC 5907) up to 20% (for UGC 4277), with an average of 12 ± 6%. Therefore, 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 for 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 in itself. For NGC 891, Kamphuis et al. (2007) and 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 2013). This points to the presence of a substantial amount of diffuse extraplanar 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 RT 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 r = 0.53), 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 structures 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, and the 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.

6. Summary

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 RT 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 (K band/IRAC 3.6 µm) 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–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 ∼10%.

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, and dust at high galactic latitudes are clearly present, as well as signatures of a clumpy, inhomogeneous medium. Therefore, we suggest that these small- and large-scale structures may be responsible for the dust energy balance problem. One should note, however, that the presence of different heating sources, which can be clearly distinguished at different wavelengths in the NIR-FIR domain (as shown in Appendix D), suggests that the nature of the dust and stellar constituents, and the link between them in galaxies, require further investigation.

The average face-on optical depth in the V band for the sample galaxies is 1.07 ± 0.53 (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 3 and 4. The IMFIT, FITSKIRT, and SKIRT models (with the constructed SEDs) are presented in Appendices A–D.

---

Acknowledgments

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 Dust-Pedia. CHARM (Contemporary physical challenges in Heliospheric and AstRo-physical 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).

References

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, therefore, can be written as(E.1)

where M _d is the dust mass, B _v is the Planck function, T _d 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 κ(λ ₀) = 0.64 m² kg⁻¹ at λ ₀ = 250 µm and β = 1.79 adapted to the THEMIS model. The fits were done by using the IDL procedure mpcurvefit.

The results of these modified black-body fits are shown in Fig. E.1 and listed in Table 7.

All Tables

All Figures