© ESO, 2011

1. Introduction

Copious quantities of interstellar dust are present in the disks of all metal rich star-forming galaxies. This dust pervades all components of the interstellar medium (ISM), ranging from the diffuse ionised and neutral medium filling most of the volume of the gaseous disk, through embedded neutral and molecular clouds of intermediate sizes and densities, down to the dense cloud cores on sub-parsec scales which are the sites of formation of stars. The ubiquity and high abundance (relative to the available pool of interstellar metals) of grains is a fundamental result of the solid state being the favoured repository for refractory elements in all but the coronal component of the ISM, and affects the very nature of galaxies. In particular, dust plays a major role in determining the thermodynamic balance of the ISM, through photoelectric heating and inelastic interactions with gas species, influencing the propensity of galaxies to accrete, cool and condense gas into stars (see e.g. Popescu & Tuffs 2010). In view of the physical role played by dust particles in the process of star formation, it can be regarded as a minor perversity of nature that, by virtue of their strong interaction cross section with stellar photons, the very same particles strongly inhibit and distort our view of the resulting stellar populations, thus preventing a straightforward confrontation of theories for star formation in galaxies with observations.

Statistical studies of the luminosity and colour distributions of large samples of galaxies seen at different orientations (most notably the attenuation-inclination relation) have shown the severity of this effect even in the optical bands, and not just in exceptionally opaque systems, but in the major part of the population of star-forming disk galaxies in the local Universe (Driver et al. 2007, 2008; Shao et al. 2007; Choi 2007; Unterborn & Ryden 2008; Padilla & Strauss 2008; Cho & Park 2009; Maller et al. 2009; Ganda et al. 2009; Masters et al. 2010). This means that we not only need to solve the problem of deriving intrinsic stellar emissivities from newly formed stars emitting predominantly in the ultraviolet, but also need to decode the effect of dust on the amplitude and colours of optical light to derive the SF history of galaxies. This is fundamentally an ill-constrained problem using broad-band optical data alone, due to the age/metallicity/opacity degeneracy (Silva & Elston 1994; Worthey 1994; Li & Han 2008). Even for galaxies bright enough for constraints on dust attenuation to be made through optical spectroscopy, for example from measurements of the emission from hydrogen recombination lines, large uncertainties still prevail due to the complex geometry of dust in galaxies and the different effect of this dust on the light from stellar populations of different ages.

In recent years the advent of spaceborne infrared astronomy has meant that we have now the capability of viewing the absorbed energy of the starlight in galaxies, which would seem to offer a solution to unravelling the effects of dust, at least from an observational perspective. Star-forming disk galaxies are observed to have most of this absorbed energy re-radiated in the far-infrared (FIR), with a significant amount also in the mid-infrared (MIR) and in the Polycyclic Aromatic Carbon (PAH) bands (see review of Sauvage et al. 2005). Simple energy balance measurements indicate that indeed a large fraction of the stellar light is absorbed even in relatively quiescent spiral galaxies in the local universe (Popescu & Tuffs 2002; Xu et al. 2006; Driver et al. 2008). Commonly, the total re-radiated energy is estimated by fitting the (typically sparsely sampled) measurements of PAH/dust emission either with phenomenological SED models (Devriendt et al. 1999; Dale & Helou 2002; Draine & Li 2007) or with SED models that use templates (Xu et al. 1998; Sajina et al. 2006; Marshall et al. 2007; da Cunha et al. 2008). However, this approach is insufficient to derive the intrinsic UV and optical emission of galaxies, even in cases where the total infrared emission can be accurately estimated.

There are several reasons why obtaining information about the relative absorption probabilities of UV/optical photons as a function of wavelength is a very complex problem, which cannot be solved by means of the application of single wavelength dependent correction factors as is the case for extinction corrections for stars¹.

Firstly, as already noted, the ISM is highly structured, with dust present in both large scale and small scale structures. This strongly affects the dependence of attenuation on UV/optical wavelength, even for uniform geometrical distributions of stars, as shown by the study of attenuation of starlight in a turbulent ISM with a log-normal density distribution by Fischera & Dopita (2005). Similarly, various studies have shown that different mass fractions of dust in clumps and diffuse structures can affect the amplitude and wavelength dependence both of the UV/optical attenuation and of the dust emission (Kuchinski et al. 1998; Bianchi et al. 2000a; Witt & Gordon 2000; Misselt et al. 2001; Misiriotis & Bianchi 2002; Pierini et al. 2004; Tuffs et al. 2004; Bianchi 2008). All these studies indicate that models for the attenuation of starlight and its re-emission in the infrared will only have predictive power if independent constraints on the geometry of the dust-bearing structures in the ISM are available.

Secondly, stellar populations of different ages have systematically different geometrical relations to these dust structures (Silva et al. 1998; Charlot & Fall 2000; Popescu et al. 2000) due to the disruption of birth clouds through the action of stellar winds and supernovae, the migration of stars away from these clouds with time and the subsequent dynamical evolution of stellar populations on Gyr timescales from a thin disk to a thick disk configuration (Wielen 1977).

Thirdly, the response of grains to light is not just a question of the wavelength and flux of the photons hitting them but also of their optical properties (dependent on shape, size and composition; see review by Draine 2009). A particularly complex aspect of this is the account that should in principle be made of any systematic differences in composition on location in the ISM, such as for example the transition between pure refractory species in the diffuse ISM irradiated by UV to ices and ice-coated grains in opaque molecular clouds. Furthermore, even supposing the composition and size distribution is known, the stochastic emission of small dust grains responding to impulsive heating has to be calculated throughout the volume, which is computationally challenging.

Arguably, the most fundamental and challenging of these problems is the determination of the geometry of the dust-bearing structures in the ISM, since knowledge of this is a pre-requisite to obtaining a self-consistent solution embracing the amplitude and geometry of the stellar populations and the optical properties of the grains. Where FIR data is available, a powerful way of constraining the geometry of dust on all scales is to use the fact that grains act as test particles with FIR colours characteristic of the intensity and colour of the interstellar radiation fields (ISRF), which in turn are a strong function of morphological structures within galaxies. For example grains locally heated by the strong radiation fields within embedded star formation regions emit an infrared radiation that peaks around 60 μm, while grains heated by the diffuse radiation fields in the disk of galaxies will emit an infrared radiation that typically peaks well longwards of 100 μm, but also exhibit substantial emission in the MIR region due to the stochastic heating of small grains (see Sauvage et al. 2005). In general using dust grains as tracers of the ISRF is a key way to constrain both opaque and translucent components of the ISM on all spatial scales, and the technical implementation of this approach is to use radiative transfer (RT) calculations. Much effort has been put into this problem of developing efficient RT codes, and there are now several available in the literature (see Kylafis & Xilouris 2005, for a comprehensive review on radiative transfer techniques used in galaxy modelling). Specific RT calculations of the integrated MIR/FIR/submm emission from grains heated by starlight in disk galaxies have been made for various specifications of the geometry of stars and dust by Siebenmorgen & Krügel (1992), Silva et al. (1998), Popescu et al. (2000), Efstathiou & Rowan-Robinson (2003), Piovan et al. (2006), Bianchi (2008).

Further empirical constraints on the geometry of dust in galaxies are possible for systems where major dust-bearing components of the ISM are translucent, thus allowing information about the distributions of dust and stars on scales of hundreds of parsecs to be inferred from images of the optical emission. This approach, which also relies on modelling of the data with RT techniques, was adopted by Xilouris et al. (1999) in their investigation of the large scale diffuse components in edge-on spiral galaxies, revealing reproducible trends in the geometry of stars with respect to the dust. The analysis of Xilouris et al. was also notable in that it was also able to measure the extinction law for the visible component of dust for these systems, showing it to be consistent with that of dust in the diffuse ISM of the Milky Way in the observed optical/NIR range. In principle the geometric constraints on the distributions of stellar emissivity and dust derived from optical observations of translucent regions can be extended to more obscured components of the ISM, such as for example the molecular layer in the disk out of which stars are forming, by making use of the physical connection between gas, stars and dust to further constrain the geometry of the problem. Overall, our knowledge of galaxies, at least in the local Universe, suggests that the specification of geometry does not need to be completely ad hoc and can be empirically constrained by the colour of the integrated dust emission, by the optical measurements of structure, and by considerations of physical plausibility. This approach can be directly applied to the panchromatic observations to decode their information into basic physical parameters of galaxies, and we therefore call this approach decoding observed panchromatic SEDs (see Popescu & Tuffs 2010).

The only alternative to invoking empirical constraints on the distributions of dust and stars in galaxies is to calculate these geometry from first principles, using numerical simulations of how galaxies form and evolve. This approach was followed by Chakrabarti et al. (2008), Chakrabarti & Whitney (2009), Jonsson et al. (2009) and we call it encoding predicted physical quantities. A particular advantage of the encoding approach is in its application to high redshift galaxies with more complex geometries due to frequent mergers and interactions, where hydrodynamical simulations are ideal for describing structures on scales greater than ca. 100 pc (Jonsson et al. 2009), at least in a statistical sense. The coupling of cosmological simulations with SED modelling tools is moreover an obvious way to interpret the cosmological evolution of luminosity functions of both the direct and re-radiated components of stellar light, something which has already been done using semi-analytic models of structure formation (Almeida et al. 2010). On the other hand, the encoding process relies on the assumption that the theory used to predict the relation between gas, dust and stars in galaxies is complete. Furthermore, deriving information on galaxies on an object-by-object base is better suited to the decoding process, at least for evolved disk galaxies in the local Universe with more uniform geometries.

In this paper we follow the decoding approach, using an empirical determination of geometry to self-consistently predict the UV/optical attenuation of starlight and the corresponding dust/PAH emission SEDs of star-forming galaxies in the local Universe. To this end we use an updated and enhanced version of the original model of Popescu et al. (2000, hereafter Paper I) which incorporates the geometrical constraints on dust and stars in spiral galaxies found by Xilouris et al. (1999). Specifically, we calculate a comprehensive library of spatially integrated dust/PAH re-emission SEDs² of disk galaxies as a function of a minimal set of physical parameters needed to predict the dust/PAH emission. We further describe how this set of dust/PAH re-emission SEDs can be self-consistently combined with a new library of UV/optical dust attenuations calculated using the same model (an update of the existing library from Tuffs et al. 2004, hereafter Paper III) to invert an observed set of broad-band photometry of a galaxy spanning the UV/optical – FIR/submm range to derive the intrinsic (i.e. as would be observed in the absence of dust) UV/optical emission of the galaxy. In particular we demonstrate that this analysis can, without significant loss in accuracy, be done independent of a priori assumptions of the detailed mathematical form of the dependence of SFR on time. This approach allows the UV/optical SED to be dereddened in a fast and flexible way without recourse to population synthesis models. The latter models can then be compared to the dereddened UV/optical SED to derive the SF history in a further, independent step, avoiding bias due to the age/metallicity-opacity degeneracy.

The primary motivation to develop this approach was to extend the applicability of a fully self-consistent RT solution to the large statistical samples of optically selected local Universe star-forming galaxies spanning a full range of mass, morphology and environment for which, thanks to facilities such as GALEX, WISE, Spitzer, AKARI and Herschel (e.g. Driver et al. 2009; Eales et al. 2010; Martin 2010) integrated photometry is now for the first time becoming routinely available across the full UV-submm range³. In addition to showing how the UV/optical attenuation is constrained through the observed colour and amplitude of the observed PAH/dust emission spectrum, we also describe how to utilise morphological information from high resolution optical observations of galaxies (such as linear sizes of disks and the bulge-disk decompositions) which are expected to come from the next generation of wide field imaging spectroscopic surveys of local Universe galaxies (e.g. Driver et al. 2009). In keeping with the decoding approach that we have adopted, we envisage that the model will be well suited to the establishment of empirical relations between SF activity and SF history of optically selected galaxies with a large range of properties of interest, such as dynamical mass (both of the parent DM halo and in disk), baryonic gas content, dust content, relativistic particle content, and environment. Amongst other applications the model should also be useful for evaluating the relation between broad-band and spectroscopic measures of both SFR and opacity in spiral galaxies so that these quantities may be measured in cases where a full panchromatic coverage of the UV-submm SED is not yet available.

Although we have calibrated the model using the optical observations of the galaxies from Xilouris et al. (1999) which are all spirals in the middle range of the Hubble sequence, we have incorporated a bulge-to-disk luminosity ratio parameter which will extend the applicability of the model to galaxies of earlier and later morphological types. We have also identified a non-parametric test to evaluate the fidelity of this model by virtue of its prediction of the colour-surface-brightness relation for spiral galaxies. This test will allow galaxies with systematic differences from the geometry adopted in this model to be identified.

This paper is organised as follows. In Sect. 2 we describe our model for the dust emission, together with the physical motivation associated with its underlying assumptions, and formulate the mathematical framework underpinning the quantitative comparison of dust emission with the attenuation calculations presented in this paper. We draw particular attention to the changes in the model compared to the original formulation in Paper I. Specifically, we describe the inclusion of PAH molecules, needed to extend the applicability of the model to shorter (mid-IR) wavelengths, and which also led to a revision to the overall grain size distribution to preserve consistency with the extinction and emission properties of interstellar dust. A further major update to the model is the inclusion of the contribution to dust heating of the young stellar population at optical wavelengths, avoiding the artificial cut-off of the emission of this population beyond the UV range in the original version of the model from Paper I. In this section we further describe a revised treatment of the emission from star-formation regions, incorporating a wavelength-dependence of the escaping UV light to achieve consistency with their attenuation properties as formulated in Paper III and including an updated template for the corresponding infrared emission using the recent model of Groves et al. (2008). In Sect. 3 we illustrate our model for the case of the well studied galaxy NGC 891. A more general quantitative analysis of the effects of the main assumptions of the model on the predicted spatially integrated SEDs as well as tests on the fidelity of the model are given in Sect. 4. We then present our library of simulated dust emission SEDs in Sect. 5 and discuss the predicted variation of these SEDs with the main parameters of the model in Sect. 6. We summarize the main results of the paper in Sect. 8.

The updated calculations of attenuation of stellar light for the diffuse component are given in Appendix B. To facilitate the use of this library of attenuation in combination with the dust emission calculations we give in Appendix C formulae for the composite attenuation of stellar light from the different geometrical components of stellar emissivity in terms of the parameters used for the prediction of the dust emission. For an easy access to potential users of the model we give in Appendix D a step-by-step guide on how to use the model for fitting real data, together with the full mathematical formalism. In this appendix we also describe the non-parametric test to objectively judge the efficacy of the model by virtue of its prediction of the colour-luminosity relation for spiral galaxies.

2. Calculation of the infrared SEDs

We calculate the infrared SEDs using an updated version of the model from Paper I which was applied to the interpretation of dust emission from edge-on galaxies by Misiriotis et al. (2001; Paper II) and also used to make predictions for the effects of dust in the optical range in Paper III and by Möllenhoff et al. (2006; Paper IV). In the following we describe the changes in detail, giving their physical motivation, and showing how they interlink with the mathematical formulation of the model.

2.1. The distribution of stars and dust

Star-forming galaxies are fundamentally inhomogeneous, containing highly obscured massive star-formation regions, as well as more extended large scale distributions of stars and dust. It is important to emphasize that the large-scale distribution of diffuse dust plays a major role in mediating the propagation of photons in galaxy disks and dominates the total bolometric output of dust emission. The discovery of this diffuse component was one of the highlights of the Infrared Space Observatory (ISO) (see Tuffs & Popescu 2006; Sauvage et al. 2006, for reviews on the ISO science legacy on normal nearby galaxies). This result was obtained both from analysing the integrated properties of infrared emission from galaxies (Popescu et al. 2002) but also from resolved studies of nearby galaxies (Haas et al. 1998; Hippelein et al. 2003; Tuffs & Gabriel 2003; Popescu et al. 2005) and is now being confirmed by the infrared data from Spitzer (Pérez-González et al. 2006; Hinz et al. 2006; Dale et al. 2007; Bendo et al. 2008; Kennicutt et al. 2009) and AKARI (Suzuki et al. 2007). Our model being empirically motivated we tried to incorporate all available observational constraints provided by the data. Accordingly our model includes both a diffuse component and a localised component representing the star-formation complexes.

A schematic picture of the geometrical components of the model is given in Fig. 1, together with a mathematical prescription of the stellar emissivities and dust opacities used in the model.

The large scale distribution of stars and dust are approximated as continuous spatial functions of stellar emissivity and dust opacity, which we refer to as “diffuse” distributions. We have separate distributions for the old and young stellar populations, and we also consider separate distributions for diffuse dust associated with these populations.

The old stellar population resides in a disk and a bulge, with its emissivity described by a double exponential and a de Vaucouleurs distribution, respectively:

(1)(2)where R and z are the cylindrical coordinates, η^disk(λ,0,0) is the stellar emissivity at the centre of the disk, , are the scalelength and scaleheight of the disk, η^bulge(λ,0,0) is the stellar emissivity at the centre of the bulge, R_e is the effective radius of the bulge, and a and b are the semi-major and semi-minor axes of the bulge.

The total luminosity of the old stellar population in the disk is then given by: (3)whereas for the bulge there is no exact analytical formula for the spatially integrated luminosity.

The dust associated with the old stellar population (“the first dust disk”) is also described by an exponential disk: (4)where is the extinction coefficient at the centre of the disk and and are the scalelength and scaleheight of the dust associated with the old stellar disk.

In our model the parameters describing the geometry of the old stellar population and of the first dust disk are empirically constrained from resolved optical and near-IR images via the results of the modelling procedure of Xilouris et al. (1999). For edge-on systems these calculations completely determine the scale heights and lengths of the exponential disk of old stars and associated diffuse dust, as well as the effective radius and ellipticity of the dustless stellar bulge. This is feasible for edge-on systems since the scale height of the dust is less than that of the stars. By analysing a sample of 5 nearby edge-on galaxies with morphological classification in the range of Sb-Sc, Xilouris et al. found that the scalelength of the dust was larger than that of the stellar disk, that the scaleheight of the stars was larger than that of the dust and that the scalelength of the stellar disk decreases with increasing wavelength. In this way Xilouris et al. was also able to find a general relation between the scaleheights and scalelength of old stars and dust and the dependence of this relation on wavelength. The derived relation enabled us to fix the relative distribution of old stars and dust in our model. In Papers I and II we verified that the derived geometrical parameters from Xilouris et al. correctly predicted the dust emission SEDs of the 5 nearby edge-on galaxies. Observationally it was also confirmed that the dust disk not only has a larger scalelength than the stellar disk (Alton et al. 1998; Davies et al. 1999), but was detected to physically extend well beyond the stellar disk (Popescu & Tuffs 2003; see also Popescu et al. 2002; Hinz et al. 2006, for dwarf galaxies). From resolved studies of galaxies it was also found that there is a large scale distribution of diffuse dust having a face-on opacity that decreases with radius (Boissier et al. 2004; Popescu et al. 2005; Pérez-González et al. 2006; Boissier et al. 2007; Muñoz-Mateos et al. 2009).

The geometrical parameters of our model are listed in Table E.1, where all the length parameters are normalised to the B-band scalelength of the disk, . In our calculations we take  pc, the fixed reference scalelength of our model galaxy, as derived for NGC 891.

The “young” stellar population and associated dust are also specified by exponential disks, which are taken to have small scaleheights (thin disks, thereby the superscript “tdisk”), which we shall refer to as the “young stellar disk” and the “second dust disk”: (5)(6)where η^tdisk(λ,0,0) is the stellar emissivity at the centre of the thin disk, and are the scalelength and scaleheight of the thin disk, is the extinction coefficient at the centre of the thin disk and and are the scalelength and scaleheight of the dust associated with the young stellar disk.

The total luminosity of the young stellar disk is given by: (7)Unlike the old components, the parameters describing the geometry of the young components cannot be constrained from images of stellar light, because the young stellar populations are highly obscured in most cases. It was thus necessary to constrain these parameters from physical considerations. The scale height of the young stars was taken to be 90 pc (the value for the Milky Way) and its scale length is equated to that of the “old stellar disk” in B-band, h^tdisk = h^disk(B). The dust associated with the young stellar population was fixed to have the same scalelength and scaleheight as for the young stellar disk. The reason for this choice is that our thin disk of dust was introduced to mimic the diffuse component of dust which pervades the spiral arms, and which occupies approximately the same volume as that occupied by the young stars. This choice is also physically plausible, since the star-formation rate is closely connected to the gas surface density in the spiral arms, and this gas bears the grains which caused the obscuration. In principle it would be more realistic to place the young stellar population and associated dust into a spiral pattern rather than into a disk. However, in practice, as we show in Sect. 4.2, this makes little difference to the predicted volume integrated dust emission SED and starlight attenuation of a galaxy. Furthermore, the inclusion of spiral arms would necessitate another model parameter to describe this more complex geometry. Also, in most practical applications for the interpretation of the integrated emission of galaxies in large statistical samples, photometric information about the spiral arm pattern is not available. For all these reasons we consider the approximation of the second dust disk to be both reliable and practical. The geometrical parameters of the young stellar disk and second-dust disk are listed in Table E.1.

Apart from the geometrical parameters described above, the distributions of diffuse stellar emissivity and dust are also described in terms of their amplitudes. The amplitudes of the two dust disks , can be expressed in terms of the central face-on opacity in the B band, ,, defined by: (8)(9)The opacity of the first dust disk can be derived jointly with the geometrical parameters from the optimisation technique of Xilouris et al. (1999). The optimisation also determines the extinction law of the diffuse dust empirically, since the calculations are done independently for each optical/NIR image. For the galaxies studied by Xilouris the derived extinction law was consistent with a Milky Way type dust, which is also the type of dust adopted for our dust model (see Sects. 2.4 and A). The opacity of the second dust disk was a free parameter in the calculations from Papers I and II, and is strongly constrained by the level of the submm emission. To minimise the number of free parameters, we fix (10)to the value 0.387 found for our proto-type galaxy NGC 891, which is also close to what was found for a second edge-on galaxy with submm data, NGC 5907, which we modeled in Paper II. We note here that the attenuation-inclination relation predicted for this fixed ratio of opacities in the two dust disks was found to successfully reproduce the observed attenuation-inclination relation of a large and statistically complete sample of galaxies from the Millennium Galaxy Catalogue Survey (Driver et al. 2007). We thus adopt as a free parameter of the model the total central face-on opacity in the B-band : (11)The remaining parameters – the amplitudes of the two stellar disks η^disk, η^tdisk and of the bulge η^bulge – and their link to the free parameters of the model are discussed in Sect. 2.3.

2.2. The clumpy component

An important component of the spatially integrated dust emission from star-forming galaxies arises from dust in the birthclouds of massive stars, as previously modelled by Silva et al. (1998), Popescu et al. (2000), Charlot & Fall (2000), Efstathiou & Rowan-Robinson (2003), Jonsson et al. (2009). Because these clouds are spatially correlated with their progeny on parsec scales, they are illuminated by a strong UV-dominated radiation field of intensity 10 − 100 times that in the diffuse ISM. This gives rise to a localised component of emission from grains in thermal equilibrium with these intense radiation fields, which, despite the tiny filling factor of the SF regions in the galaxy can nevertheless exceed the entire diffuse infrared emission of a galaxy at intermediate wavelengths (ca. 20 to ca. 60 μm). It is therefore particularly important to incorporate a clumpy component of dust associated with the opaque parent molecular clouds of massive stars. Following Popescu & Tuffs (2005) we refer to these clumps as “active clumps”. The active clumps are assumed to have the same spatial distribution as the young stellar disk and the second dust disk. Furthermore, it is assumed that the properties of these clumps do not systematically depend on their radial location within the galaxy. In reality we expect star formation complexes in more pressurised regions (such as the inner disk) to be more compact and therefore have warmer FIR colours than their counterparts in low pressure regions (such as the outer disk), as modelled by Dopita et al. (2005). And indeed this phenomenon has been observed in HII regions in M 33 (Hippelein et al. 2003). However, bearing in mind that we will empirically calibrate our template (see Sect. 2.8) on the whole ensemble of star forming complexes in the Milky Way, the assumption we make in this regard should not significantly affect our predictions for the integrated emission of galaxies.

Since birthclouds of stars are typically fragmented due to the combined effects of supernovae, stellar winds, and the general motion of stars away from the clouds, only a certain fraction of the total luminosity of massive stars in a galaxy will be locally absorbed. Following the original formulation in Paper I we denote this fraction in our model by the “clumpiness factor”, F. Since birthclouds are completely opaque at all UV/optical/near-IR wavelengths (e.g. Sievers et al. 1991), F can be physically identified with the luminosity-weighted mean fraction of directions from the massive stars, averaged over the lifetime of the stars, which intersect the birthcloud. This concept, first introduced by Popescu et al. (2000), allows UV light to freely stream away from star-forming regions in some fraction of directions, allowing the diffuse dust to be illuminated with more UV photons than would have been the case if the young stars had been assumed to be completely cocooned in their parent molecular clouds. This is in qualitative accordance both with the high frequency at which counterparts of star-formation regions seen at 24 μm are detected in the UV, as well as with the predominance of PAH emission from diffuse dust in galaxies correlated with cold dust (Bendo et al. 2008).

A more detailed description of the clumpy component is given in Sects. 2.5.1 and 2.8.

An approximation of our model is that the heating of the grains in the active clumps is dominated by photons from the stellar progeny, and we neglect any external contribution from the diffuse ambient radiation fields in the galaxy. This should be an excellent approximation for spiral galaxies, where the filling factor of star formation regions is small. We note that if this approximation were invalid, we would be forced to perform radiative transfer calculations with parsec resolution to properly describe the heating of grains in the optically thick birth clouds, which would render the calculation of a comprehensive set of SEDs as given in this paper intractable with current computing resources. To date, the best resolution achieved with an adaptive grid code when modeling the observed dust emission SED from an individual galaxy is about 20 pc (Bianchi 2008). In any case, we would only expect collective effects, whereby photons from adjacent SF regions provide a significant fraction of the dust heating, to be significant in galaxies with very high volume densities of SF regions, such as the central region of a starburst galaxy. The low filling factors of opaque clouds in spiral disks is supported by high resolution surveys of the Galactic Plane with large ground-based telescopes in the submm (e.g. the APEX Telescope large Area Survey – ATLAS; Schuller et al. 2009). At these wavelengths the clouds are optically thin, thus directly tracing the total column density of dust, yet only sporadic peaks of dust emission can be seen with local enhancement in visual optical depth much greater than 1. This contrasts with surveys in the CO lines (Matsunaga et al. 2001), which show a more highly filled distribution of emission due to optical depth effects in the radio molecular line. Recently the new Herschel Infrared Galactic Plane Survey (Hi-Gal) (Molinari et al. 2010) has revealed the fainter, more extended emission components which dominate the morphology of the dust emission in the Milky Way, showing again that the optically thick molecular cores have only a small filling factor. In summary, based on all available observational evidence we believe it is a reasonable approximation to ignore the external heating of active clumps.

Another approximation of our model is that the clumpy distribution of dust is exclusively associated with the opaque parent molecular clouds of massive stars – the active clumps. In reality, some of the dust in the diffuse dust disks may also be in clumps without internal photon sources, and thus heated only by the diffuse ambient radiation fields in the disk. Following Paper I (see also Popescu & Tuffs 2005) we refer to such clumps as “passive” or “quiescent” clumps. Provided passive clumps are optically thin, the transfer of radiation through the disks will be virtually identical to that in a homogeneous disk. However, once the passive clumps become optically thick, the self-shielding of grains will yield a solution with a reduction in both the attenuation of the stellar light and the infrared emission compared to the case where the same mass of grains is diffusely distributed. These effects have been quantified by Bianchi et al. (2000a) who showed that the shape of the infrared SED of a disk galaxy can be strongly affected shortwards of 200   μm (see also Misselt et al. 2001, for passive clumps distributed in a spherical shell geometry). We do not include passive clumps in our model for both physical and empirical reasons. Physically, it seems likely that most passive clumps will be optically thin, since as soon as they become optically thick to the impinging external UV light, they will lose their principal source of heating (the photoelectric effect) and will be prone to collapse and form stars (Fischera & Dopita 2008). Empirically, this presumption is supported by the recent findings of Holwerda et al. (2007a,b), that the structure of the diffuse ISM consists of optically thin dusty clouds. Furthermore, the change in the shape of the predicted dust emission SEDs imposed by the incorporation of passive clumps tends to provide a colder solution than needed to fit real data (Bianchi 2008).

2.3. Constraints on the intrinsic SEDs of the stellar populations

As mentioned in Sect. 2.1, the amplitudes of the stellar populations of the two stellar disks (η^disk(λ,0,0) and η^tdisk(λ,0,0)) are parameters that are yet to be determined. In other words we want to find a solution for the attenuation of stellar light and dust emission SED that can fit observed SEDs to provide resulting intrinsic stellar SEDs. At the same time, in our modelling procedure we deliberately do not want to use population synthesis models to fully fix the intrinsic SED, as we want to be as free as possible of any assumptions about the SF history. Nonetheless, some constraints on the input SEDs need to be made to avoid having the amplitude of the stellar SED a free parameter at each of the sampled wavelengths listed in Table E.2. Since the main factors shaping the dust emission SEDs are the total luminosity of the young stellar population and of the old stellar populations, we choose to have these two quantities as free parameters of our model and produce calculations for all combinations of these two variables. We then assume that each of the two components (young and old stellar populations) have a fixed wavelength dependence (a fixed template SED), thus reducing the number of free parameters describing the intrinsic SEDs of the stellar populations to two. We will show in Sect. 4.3 that these assumptions have a negligible effect on the predictions for the dust emission SEDs and the associated parameters and F. In turn, the solution for the wavelength dependence of the attenuation in the UV-optical range, which depends on , F (and inclination), will be similarly secure, thus allowing the true dereddened stellar emission SED to be recovered. In effect we are using the fact that the dependence of the dust emission SED on the colour of the stellar photon field depends primarily on the ratio between the luminosities of the young and old stellar populations rather than on the detailed colour of the emissions from either of these populations. This way of constructing the model will allow SF histories to be extracted from the dereddened stellar emission SED in a separate step from the extraction of and F.

In constructing the template SEDs of the stellar populations we consider the term “optical”, “UV”, “ionising UV” and “non-ionising UV” to denote the wavelength ranges λ ≥ 4430   Å, λ < 4430   Å, λ ≤ 912   Å, 912 < λ < 4430   Å, respectively, thereby marking the boundary between the UV and optical regime in the B band.

2.3.1. The template SED of the old stellar population in the disk

In defining the fixed shape of the SED of the old stellar population in the disk we only consider optical radiation and neglect any contribution in the UV. The precise wavelength dependence of this SED was fixed to the empirical relation obtained by Xilouris et al. (1999) from fitting simulated images produced by radiative transfer calculations to the observed B, V, I, J, K images⁴ of NGC 891. The SED of the old stellar population in the disk, as “calibrated” on NGC 891 is given as the light blue curve in Fig. 8. In applications to other galaxies this curve needs to be scaled according to the total output of the old stellar population in the disk of the modelled galaxy. For this purpose we use a unitless parameter old defined as (12)where  W, which corresponds to 10 times the luminosity of the non-ionising UV photons produced by a 1   M_⊙/yr young stellar population, as defined in Sect. 2.3.3. Thus old is the normalised luminosity of the old stellar populations in the disk and is adopted to be one of the free parameters of the model. In this scheme the best fit solution for NGC 891 corresponds to old = 0.792.

The intrinsic spectral luminosity densities of the old stellar population in the disk corresponding to old = 1, , used in our calculations are listed in Table E.2. The spectral luminosity density of the old stellar population in the disk of a model galaxy is then given by: (13)where has been previously defined in terms of η^disk(λ,0,0) in Eq. (3), (14)and the value of can be derived by integrating over wavelength⁵ the values of from Table E.2.

As already mentioned, we will show in Sect. 4.3 that, although the template SED of the old stellar population used in the model is fixed to a single shape (here defined empirically on the galaxy NGC 891), it is still valid to use it to model the dust emission SEDs of other galaxies which may have very different stellar emission SEDs.

2.3.2. The template SED of the old stellar population in the bulge

The wavelength dependence of the stellar luminosity produced by the bulge stellar population is simply linked to that of the old stellar population in the disk via a wavelength independent bulge-to-disk ratio B/D, one of the free parameters in our model: (15)

2.3.3. The template SED of the young stellar population

The wavelength dependence of the stellar luminosity produced by the young stellar population in the UV cannot be constrained empirically, as this population is heavily obscured by dust. Because of this we used population synthesis models. We used the models from Kotulla et al. (2009), making standard assumptions for the star-formation history, namely an exponentially declining star-formation rate with a time constant τ = 5   Gyr, solar metallicity and a Salpeter IMF with an upper mass cut-off of 100   M_⊙. We emphasise that, although arbitrary, this choice of parameters will not significantly bias the model dereddening of the observed UV/optical SEDs, since, as we will show in Sect. 4.3, using different assumed shapes for the stellar emissivity SED does not significantly affect the predictions for the dust emission SEDs.

We also consider that a fraction f_ion − uv = 0.3⁶ of the ionising UV photons (emitting shortwards of 912   Å) from the massive stars within the HII regions are locally absorbed by dust. This makes only a minor contribution to the dust heating, an order of magnitude less than the contribution of the non-ionising UV photons.

A new feature of the model is the consideration of the optical emission from the young stellar population embedded in the second dust-disk. Previously the emissivity function of this population had been artificially truncated longwards of the B band. In order to fix the shape of the SED template for the young stellar population in the optical we again make use of the best fit solution to NGC 891. Specifically, the optical radiation from this template at each wavelength sampling point was fixed to be the residual between the prediction of the population synthesis model for the best fit of NGC 891 and the empirical SED of the old stellar population (derived from the best fit of NGC 891 as explained above)⁷. The resulting optical part of the SED template of the young stellar population was then combined with the UV part (which is simply the population synthesis SED for the best fit SFR for NGC 891) to fix the shape of the template over the full UV/optical range. For use in other model galaxies, the overall amplitude of this SED of the young stellar population, as “calibrated” on NGC 891, needs to be scaled according to the SFR of that galaxy, adopting SFR as a free parameter of the model.

For logistical purposes only, we chose as a unit for the luminosity of the young stellar population  W, which is the luminosity produced by a SFR of 1   M_⊙/yr (for the standard assumptions about the SF history described before) in the range (912 − 50   000)   Å. can be given as a sum of the luminosity of the non-ionising UV photons  W and the luminosity of the optical photons emitted by the young stellar disk  W. As already introduced in Sect. 2.3.1, we also defined the unit for the luminosity of the old stellar population to correspond to 10 times the luminosity of the non-ionising UV photons produced by a 1   M_⊙/yr young stellar population. In addition we also define  W for the luminosity of the ionising UV photons. We note that is not included in the definition of .

We then parameterised the calculations in terms of SFR(16)where the SFR is linked to the normalised luminosity of the young stellar population in the thin disk.

In an equivalent way we can link the normalised luminosity of the young stellar population in the UV to the SFR: (17)where  W, as defined previously, and (18)The intrinsic spectral luminosity densities of the young stellar population corresponding to SFR = 1   M_⊙   yr^-1, , used in our calculations are listed in Table E.2.

The spectral luminosity density of the young stellar population in the thin disk of a model galaxy is then given by: (19)where has been previously defined in terms of η^tdisk(λ,0,0) in Eq. (7), (20)and the value of can be derived by integrating over wavelength the values of from Table E.2.

2.4. The dust model

The dust model – i.e. the prescription for the optical properties of the grains, the chemical composition and the grain size distribution – was updated to include PAH molecules, which are regarded as the carriers of the unidentified infrared emission features commonly seen in the mid-IR emission spectrum of star-forming galaxies. As shown by Zubko et al. (2004), there is no unique model incorporating PAHs that can simultaneously fit the main observational constraints on such models, namely the extinction and emission properties of the diffuse cirrus in the Milky Way. Several models have been proposed in the literature, including the models from Zubko et al. (2004), the model of Fischera & Dopita (2008) which considers an extra component of iron grains with the optical properties from Fischera (2004), and the model of Weingartner & Draine (2001; see also Li & Draine 2001; Draine & Li 2007). Although not unique in terms of reproducing the properties of the galactic cirrus, these models do nevertheless differ in some predictions relevant to the modelling of the infrared emission of spiral galaxies, particularly with regard to whether conducting particles are included as a constituent of the dust model. For example in the model of Fischera & Dopita a substantial part of the diffuse 60 μm emission is powered by optical photons absorbed by iron grains in equilibrium with the ambient radiation fields (see also Chlewicki & Laureijs 1988), whereas in the model of Weingartner & Draine (2001) this emission mainly arises from carbonaceous grains stochastically heated by UV photons. Potentially, therefore, the choice of grain model can influence the inferred contributions of young and old stellar populations in heating the grains. We are however not aware of any direct evidence for the existence of pure metallic grains in the diffuse interstellar medium of galaxies. Furthermore, simple considerations of potential redox reactions indicate that any such grains in the diffuse neutral ISM will be susceptible to oxidation and will revert to the properties of non-conducting grains (Duley 1980), like silicates. Here we substituted the model of Draine & Lee (1984) and Laor & Draine (1993) used in Papers I − IV to the latest model from Weingartner & Draine (2001) and Draine & Li (2007) incorporating a mixture of silicate, graphite and PAH molecules.

We note that the adopted dust model is only used for the calculation of transfer of radiation through the diffuse interstellar medium, where the properties of the grains are reasonably well constrained. When including the infrared emission from star-forming complexes, which contain dense clouds where grains may form ices and other complex compounds, we use a template model SED function which is observationally constrained, as described in Sect. 2.8. This sidesteps the uncertainties in the optical properties of grains in dense clouds.

For the parameters of the model we consider Case A and R_V = 3.1 from Weingartner & Draine (2001), with the updates from Draine & Li (2007) (revised size distribution and optical constants for PAHs). The parameters of the model considered here are summarised in Table E.3.

Following Li & Draine (2001) we allow the relative abundance of neutral and ionised PAHs to vary with molecule size but to be independent of position in the galaxy (or intensity of the radiation field). The ionisation fraction (Φ_ion) of PAHs was fixed according to the average over the three phases of the diffuse ISM given in Fig. 7 of Li & Draine (2001). The PAH-size distribution was then multiplied with Φ_ion and 1 − Φ_ion, respectively, to obtain the size distribution of ionised and neutral PAH-molecules, respectively.

Since for the dust and PAH emission we accurately take into account the stochastic heating of the grains by the ambient radiation fields (see Sect. 2.6), we also need to derive the heat capacities for our dust grains. The heat capacities for silicate and graphite grains needed to derive the temperature distributions are summarised in Popescu et al. (2000). The heat capacities for the PAH-molecules were calculated using Eq. (15) of Li & Draine (2001) and Eq. (5) from Li & Draine (2002).

Finally, to facilitate the incorporation of future improvements in the dust model/inclusion of new dust models, we have constructed a flexible interface between the code for calculating dust emission and the input dust model.

Because we have modified the dust model we also had to redo the calculations of attenuation from Paper III. From the perspective of the extinction properties of the two dust models (the old version and the new updated version), they are both designed to reproduce the observed extinction curve of the Milky Way. However the discrepancy arises from the change in the relative contribution of absorption and scattering to the total extinction. This will give rise to a small but non-negligible difference in the overall energy balance (around 10% for NGC 891 – see Appendix B). In Appendix B we briefly describe the attenuation calculations obtained with the revised dust model and compare with the results from Paper III.

2.5. Calculation of radiation fields

In the formulation of our model an explicit calculation of radiation fields is only done for the diffuse component. There are two channels by which the diffuse interstellar medium can be illuminated by stellar light, since, as described in Sect. 2.1, our model has a diffuse and a clumpy distribution of stars and dust. The first channel is through the illumination of the diffusely distributed dust by the smoothly distributed population of old stars. The second channel is through the escape of radiation from the young stars, out of the localised star-forming complexes.

Thus, the first step in the calculation of radiation fields in the diffuse medium is to determine the amplitude and SED of the stellar light which escapes from star-forming complexes.

2.5.1. The escape fraction of stellar light from the clumpy component

In Sect. 2.2 we introduced the factor F, which can be physically identified with the luminosity-weighted mean fraction of directions from the massive stars, averaged over the lifetime of the stars, which intersect the birth-cloud. The high opacity assumed by our model for the birth clouds bestows a grey body absorption characteristic for the locally absorbed radiation in any individual star-formation region. There is nevertheless a wavelength dependence of the probability of absorption of stellar photons because stars of different masses survive for different times, such that lower mass and less blue stars spend a higher proportion of their lifetime radiating when they are further away from their birthclouds, and because of the progressive disruption of the clouds by supernovae and stellar winds. To calculate this wavelength dependence we follow the analysis of Appendix A of Paper III, thus preserving consistency between the model of the attenuation of light from the young stellar population given in Paper III and the model of the infrared emission from star-forming regions given here. Denoting the wavelength dependence of the locally absorbed radiation by the function F_λ, the spectrum and amplitude of the starlight from the young stellar population in the thin disk with spectral luminosity density which is intercepted by a parent cloud is thus given by: (21)where (22)and f_λ is tabulated in Table A.1 of Paper III. As described in Paper III, the tabulated values of f_λ correspond to a secular decrease of the solid angle subtended by the parent clouds on a timescale of 3 × 10⁷ yr. For completeness we remind the reader that is defined in terms of the free parameter of the model, SFR, in Eq. (19) (see Sect. 2.3.3)⁸.

The amplitude and spectrum of the light from the young stellar population which illuminates the diffuse dust, is given by (23)Since F_λ must satisfy F_λ ≤ 1.0 at all wavelengths, our model has an intrinsic constraint on the clumpiness factor of (24)Physically, this constraint corresponds to the requirement that the “porosity factor” p₀ of the birthcloud, used in Eq. (A.2), (A.13) of Paper III to denote the fragmentation of the cloud, is always less than unity. This limit denotes a complete lack of fragmentation of the cloud due to the action of supernovae and stellar winds, such that the probability of escape of photons into the diffuse ISM is determined only by the migration of stars away from their birthclouds.

2.5.2. Combining radiation fields from the young and old stellar populations

Since all the diffuse stellar components are attenuated by the same distribution of dust, the diffuse radiation fields seen by each grain can be considered to be a sum of the radiation fields produced by each stellar component. Thus we calculated the diffuse radiation fields separately for the three main stellar components of our model, the young stellar disk, the old stellar disk, and the bulge. Then the individual diffuse radiation fields can be combined to produce the total diffuse radiation fields of our desired model galaxy. This is the same concept as that used in Paper III for the calculation of the attenuation of stellar light.

To achieve this we first ran the radiative transfer calculations for each stellar component, each wavelength and each given optical depth , using the intrinsic spectral luminosity densities and , as tabulated in Table E.2. This resulted in the calculation of the unit energy densities for the disk , for the thin disk and for the bulge for the unit stellar luminosities (i.e. for the disk or the bulge and for the thin disk).

Then the unit radiation fields were scaled according to the luminosity of each diffuse stellar component, according to the parameters of the model. For the old stellar disk (only diffuse) the resulting radiation fields were obtained from: (25)where old is the normalised luminosity of the old stellar population, defined by Eq. (12). For the young stellar disk is derived from: (26)where SFR is related to the normalised luminosity of the young stellar populations in Eq. (16). For the bulge the radiation fields are given by: (27)The total diffuse radiation fields are given by:(28)\arraycolsep1.75ptThis formulation provides an exact treatment of the spatial variation of the colour of the radiation field in galaxies, which is a crucial factor in shaping the infrared SEDs.

In principle this formalism implies that we would need to create a library of radiation fields calculated for each possible value of the F factor, implying a different wavelength dependence (1−Ff_λ) of the fraction of photons escaping from the clumpy component into the diffuse one, making the problem quite intractable. In practice, however, it is possible to simplify the problem if we make the assumption that the wavelength dependence (1−Ff_λ) is fixed to the (1−F_calf_λ) corresponding to a fiducial calibration value F_cal = 0.35, which we consider to be a typical value for spiral galaxies in the local universe. This intermediate value of F_cal (theoretical values of F can vary from 0 to 0.61; Paper III) is consistent with the observation that 60 − 70% of the dust emission from spiral galaxies emanates from the diffuse component (e.g. Sauvage et al. 2006). Furthermore, detailed analysis of face-on spiral galaxies with resolved star-formation regions in the Spitzer 24 μm and the GALEX FUV bands (e.g. Calzetti et al. 2005) have shown an almost one-to-one correspondence between the direct and re-radiated light coming from these complexes, albeit with a large variance in the ratio of the two emissions. Thus there were essentially no Spitzer sources with no FUV counterparts (meaning that F is unlikely to take the value 0.61), but also no FUV sources with no Spitzer counterparts (meaning that F is unlikely to take the value 0). We note that the value of F_cal is only used to determine the exact wavelength dependence of that part of the light from the young stellar population which escapes the star-formation complexes and actually illuminates the diffuse dust. Although the computed colour of this light will be inaccurate if F differs from F_cal, corresponding inaccuracies in the predictions for the dust emission SEDs are unlikely to be more than a few percent, as demonstrated in Sect. 4.3, where we quantify the effect on the dust emission of changing the spectral shape of the stellar emissivity. We note that, for the same reason, when dereddening the observed UV/optical SED, the wavelength dependence of the component of attenuation from the clumpy component can use the exact solution for F found in the optimisation.

The effective reddening law of the continuum photons escaping the star-formation regions corresponding to F_cal is given in Table E.4.

The use of the calibration factor F_cal means that in practice all the equations that describe the illumination of the diffuse dust by the young stellar disk need to be rescaled to accommodate different values of F than those used in the calibration. For this we define a correction factor for the diffuse component corr^d(F): (29)In this case Eq. (26) becomes: (30)In an analog way we define a correction factor for the localised component corr^l(F): (31)in which case Eq. (21) becomes: (32)

2.5.3. The radiative transfer calculations

One of the key elements of any self-consistent model for the SEDs of galaxies is the incorporation of a radiative transfer (RT) code that can be used to determine the radiation fields in galaxies and their appearance in the UV/optical bands. Four methods have been proposed so far in RT models, though only two of them have been broadly used in the context of galaxy modelling, namely ray-tracing methods (Kylafis & Bahcall 1987; Silva et al. 1998; Semionov & Vansevic˘ius 2005a; Semionov & Vansevic˘ius 2006,b) and Monte-Carlo (MC) techniques (Witt et al. 1992; Bianchi et al. 1996; de Jong 1996; Gordon et al. 2001; Baes & Dejonghe 2002; Baes et al. 2003; Baes et al. 2005; Jonsson 2006; Baes 2008; Bianchi 2008; Chakrabarti & Whitney 2008)⁹.

Since our model makes explicit calculations of radiation fields in galaxies we use an ray-tracing code, which is better suited to this type of problem. As in Paper I we use the Kylafis & Bahcall (1987) code, updated and revised according to the specific needs of our modelling technique.

The code of Kylafis & Bahcall (1987) was originally developed to model the surface brightness distribution of edge-on galaxies with cylindrical symmetry. This utilises a ray-tracing method that exactly calculates the direct and first order scattered-light and uses the method of scattered intensities to approximate the higher order scattered light (see review of Kylafis & Xilouris 2005). Several modifications were incorporated for the present application, the most significant of which we briefly describe here.

Firstly, the code used was generalised and optimised for application to any cylindrically symmetric distribution of emissivity and opacity, such as the thin disk representing the young stellar population and associated dust in our standard calculations, and the annular representation of this disk that we use in the tests made in Sect. 4.2 to evaluate the second disk approximation to spiral structure. To achieve this, we incorporated an in-ray sampling scheme which automatically adapts to variations in emissivity and opacity encountered along the rays. Further, in regions of low emissivity, where the method of scattered intensities becomes less accurate, the maximum order of the scattered light had to be reduced. The reduction was done progressively accordingly to the emissivity encountered, such that only the fully accurate first order scattering term was utilised in the regions with the smallest emissivity. This is a good approximation for the geometries of dust and stars considered here, in which low emissivities are never encountered in regions of high opacities. Lastly, for application to the calculation of radiation energy densities the code was further modified to calculate the vector radiation field sampled in each direction of a predefined set of 216 directions.

To provide a more frequent sampling of the inner part of the galaxy the spatial positions were arranged on an irregular grid with sampling intervals in radial direction ranging from 50 pc in the centre to up to 2 kpc in the outer disk and with sampling intervals in vertical direction ranging from 50 pc in the plane to up to 500 pc in the outer halo. This scheme was necessary to properly sample the radiation in the more optically thick central regions, where, as can be seen from Fig. 4, the dust associated with the thin stellar disk in conjunction with the strong positional dependence of the emissivities of the young stellar population (in z) and the inner emissivity cusp of the bulge (in both R and z) can lead to quite small scale structures even in the diffuse radiation fields. We made some further tests with higher sampling, ranging from 25 pc to up to 1 kpc in radial direction and from 25 pc to up 250 pc in vertical direction. This made a less than 1% difference in the results with a large increase in the computation time. We have therefore concluded that the optimum sampling of the radiation fields is the ΔR = 50 pc−2 kpc and Δz = 50 pc−500 pc scheme.

The database of integrated SEDs presented in this paper does not incorporate the effect of the reabsorption of infrared light emitted by grains in the galaxy, as this process has only a marginal effect on the integrated dust emission, and would entail a significant increase in complexity and decrease in speed of the calculations. This is a good approximation for the present problem. For example, for the most opaque disks we consider (8), only the  ~3   μm dust emission shows a significant depression in the integrated emission. However we do calculate reabsorption of infrared light when we produce maps of galaxies seen at higher inclination, and this process has been incorporated in the calculation of the profiles presented in Sect. 4.1.

The radiation fields illuminating the dust in the diffuse disk are calculated taking into account anisotropic scattering, using the albedo and anisotropy of the scattering phase function as defined in Appendix A and the geometry of the dust as defined by Eqs. (4) and (6) in Sect. 2.1. At each UV/optical wavelength¹⁰ of Table E.2 and each optical depth a separate calculation was made for each of the three stellar components, with smoothly distributed emissivities given by Eq. (1) for the old stellar population in the disk and in the bulge, and Eq. (5) for the young stellar population in the thin disk. This adequately samples the wavelength variation in the optical constants of grains; tests done on this indicate uncertainties at the level of a percent. Concerning the corresponding error due to the finite sampling of the stellar SEDs, the reader is referred to Sect. 4.3, in which we show that the exact colours of either young or old component of the stellar SEDs only affect the resulting dust emission SEDs at the level of a percent. Overall, therefore, we consider that uncertainties introduced by our finite wavelength sampling on the predicted dust emission SEDs are minor.

When calculating the heating of the grains, and therefore of the wavelength-integrated energy absorbed, the radiation fields were interpolated on a fine wavelength grid containing 500 data points logarithmically sampled between 912 and 50 000 Å. Tests done with 1000 data points showed negligible differences in the results, therefore we decided to fix the wavelength sampling to 500. Energy balance tests done within a common wavelength sampling indicate that our models conserve energy to within better than 1 percent.

2.6. Calculation of grain temperature distributions

As described in detail in Popescu et al. (2000) our calculation incorporates an explicit treatment of the temperature fluctuations undergone by small grains whose cooling timescales are shorter than the typical time interval between impacts of UV photons. This so-called “stochastic heating” process determines the amplitude and colour of the bulk of the diffuse emission from most spiral galaxies in the shortest infrared wavelength bands (<100   μm). Thus, since for each position in the galaxy the energy density of the radiation fields heating the grains has been derived (see Sect. 2.5), we can calculate the probability distribution of temperature P_i(a,T) for each grain size a and composition i by equating the rate of absorption of energy to the rate of emission of energy:(33)\arraycolsep1.75ptwhere c is the speed of light and B_λ is the Planck function.

The probability distribution of temperature is calculated following Voit (1991), which combines the numerical integration method of Guhathakurta & Draine (1989) and a stepwise analytical approximation. We sample the temperature distribution in 61 logarithmically spaced points.

2.7. Calculation of infrared emissivities and spatially integrated SEDs in the diffuse component

Once the temperature distribution of a grain is known, the brightness I_λ,   i of a grain of radius a and composition i is then given by: (34)whereI_λ,   iis in units of  [ W   m^-2   sr^-1   Å^-1 ] .

Once the infrared brightnesses I_λ,   i are calculated for each grain size and composition, we can then integrate over the grain size distribution n(a) to obtain the infrared emissivity per H atom of grains of a given composition i, : (35)whereis in units of  [ W   sr^-1   Å^-1   H^-1 ] .

Finally, the total infrared emissivity per H atom, is obtained by summing over the grain composition i(36)In Fig. 2 we show the calculation for the infrared emissivity of grains heated by the local interstellar radiation fields (LIRF), with the energy densities of the radiation fields taken from Mathis et al. (1983). The resulting infrared SED is similar to that presented in Fig. 12 of Draine & Li (2007) and indicates that our method for calculating stochastic heating of grains gives similar results to the method of Draine and Li. This can be taken as a benchmark test for the calculation of infrared emission. In Fig. 2 we also plotted the colour-corrected DIRBE/COBE data for the high Galactic latitude. As expected, there is a good match between the model and the observations. In Fig. 2 we also show the predicted contribution of the different dust compositions to the total infrared emissivity.

Once the total infrared emissivities per H atom have been calculated, these can be then scaled to the volume density of dust grains at each position in the galaxy to obtain the volume luminosity density of the first and second dust disks, and , as a function of position: (37)(38)where and are in units of  [ W   m^-3   Å^-1 ]  and C_ext(B) = 6.38 × 10^-22   cm²   H^-1 is the extinction cross-section in the B band (at 4430   Å) for the dust model considered here. The spatially integrated SEDs (the spectral luminosity density) of the diffuse component of the model galaxy (with the fixed size ), is then: (39)We checked the fundamental energy balance in the calculations, by comparing the energy absorbed from the radiation fields with the energy emitted in the infrared, obtaining consistency to within better than 1%.

2.8. The infrared emission from the star-forming regions

In our model we take the total energy absorbed by the birthclouds to be (40)where is defined in Eq. (32), is defined in Eq. (18) and f_ion − uv = 0.3.

That is, as justified in Sect. 2.2, we take the heating of the grains in the birthclouds to be dominated by photons from the stellar progeny, and neglect any external contribution from the ambient radiation fields in the galaxy.

To determine the SED of the reradiated light from dust grains and PAH molecules in the birthclouds we use a fixed spectral template of a photodissociation region (PDR): (41)where is the template function for a photodissociation region, normalised such that (42)We use the model of Groves et al. (2008) (see also Dopita et al. 2005), fitted to broad-band measurements of dust and PAH emission from galactic star-formation regions, to provide an empirically constrained prediction for the detailed spectral form of the template . Although this model was primarily developed for use in predicting the SED of starburst galaxies, its fundamental constituent is a prediction of the PAH/dust and nebular line emission from the HII region and PDR components of individual SF regions, and so, with suitable choice for the values of the model parameters, is also directly applicable to the prediction of the PDR emission from SF regions in spiral galaxies. For the spectral template used here we only consider the PDR component calculated by Groves et al., ignoring the emission component from the HII region (as justified below). We also only include the PAH and dust emission components in the spectral template, not including the free-free and line emission from the gas phase.

The prime physical motivation for the use of the model of Groves et al. is that it self-consistently calculates the effect on the emergent dust and PAH emission SED of the dynamical evolution of the emitting regions. Specifically, the model considers a spherically symmetric, fully enclosed mass-loss bubble, driven by the mechanical energy input through winds and supernovae of a star cluster as a function of the external density (parameterised in terms of the pressure of the ambient ISM) and a “compactness parameter”, C, which scales according to the mean luminosity-weighted flux of photons onto the PDR. The emission from PAH and dust in the PDR is then self-consistently determined through a separate radiation transfer calculation (Groves et al. 2008), as a function of the grain column density in the PDR (which in the formulation of Groves et al. is actually prescribed through a combination of the gas column in the PDR and the metallicity). Being based on a dynamical model, the calculations of Groves et al. explicitly take into account the expected variation in colour and amplitude of the incident photon flux on the PDR due to the increase in radius and evolution of the stellar population, which acts to broaden the SED of the dust emission from the large ensemble of SF regions at different evolutionary stage expected in a spiral galaxy. The model also considers the photodestruction of PAH, strongly influencing the relative amplitudes of the PAH and far-IR/submm emission components from the PDR.

To empirically verify these predictions of the model, and determine the best model parameter to determine the PAH and dust emission template from an ensemble of localised SF regions in spiral galaxies, we compared the model predictions with observed data (see Table 1) for the radio-selected sample of star-formation regions in the Milky Way of Conti & Crowther (2004). This sample is distributed throughout the disk of the Milky Way, so is likely to be representative of the population of SF regions in the galaxy. In Fig. 3 we plot the mid-IR/far-IR/submm SEDs of 57 galactic star-forming regions from Conti & Crowther, normalised to the 100 μm flux density. Overlaid on this is the best fit model for the PDR emission from Groves et al., found by searching the parameter space in the compactness factor C and the dust column density in the PDR (adjusted by varying the product of metallicity and gas column density). The other major parameter, pressure, was kept fixed as the PAH and dust emission is completely insensitive to this parameter which only affects the emission lines (see Fig. 4 of Groves et al.), which are not required in our template. In comparing the model with the data we did however take into account the free-free emission component from the PDRs, since this emission component will also be represented in the broad-band continuum data of the Milky Way star formation regions. The free-free emission slightly raises the level of the emission near 1 mm above the prediction of the dust emission template, as shown by the divergence of the dotted line, which includes the free-free component from the full line in Fig. 3.

The best fitting PDR model is for compactness parameter log (C) = 6.5, metallicity 1.0, and hydrogen column density log (N) = 22.0. Although there is a large spread in observed colours, indicating the wide range of evolutionary states and environment of individual SF regions, overall this model reproduces the mean colours of the population of sources over the whole spectral range covered by the data. This spectral range extends from 12 μm (including the PAH emission), through 25, 60 and 100 μm (where emission is dominated by large grains in thermal equilibrium with the intense local radiation fields), and into the submm/near-mm range (where the emission is dominated by dust which is cold due to the self-shielding of grains in the PDR. This latter cold dust emission component is difficult to predict theoretically, depending on the otherwise poorly constrained dust column density, so is fundamentally an empirically determined quantity. Likewise, the relatively high value of the compactness factor C is needed to fit the rather warm observed far-IR/mid-IR colours. As directly confirmed through high resolution imaging of the sources this indicates a close proximity (typically parsec scales) of the exciting star clusters to the PDRs, most of which are associated with rather dense fragments of molecular material which will be only slowly pushed out from the cluster. It is interesting to note that the fit of the PDR model to the data leaves no room for a warmer dust emission component from an enclosed HII region, suggesting that the HII regions typically extend to larger spatial scales. On this picture, the cloud fragments carrying the PDR emission are embedded in a diffuse ionised medium which blends in with the diffuse emission from the second dust disk.

Although, as we argue above, the physical values for the compactness parameter C, metallicity and hydrogen column density used for the fit seem plausible, we note that there may be room for some systematic uncertainties in these parameters if the optical properties of the grains differ from the model of Weingartner & Draine used by Groves et al. The most important point however in the context of our use of this template SED is that it gives a good fit to a representative sample of star formation regions in our Milky Way. So even if the physical interpretation of the shape of the SED may be still somewhat imprecise, the template itself will be the correct one. In conclusion, we believe that the PDR model of Groves et al. that best fits the galactic star-forming regions can be taken as a template SED for the clumpy component of our model.

2.9. The free parameters for the calculation of the infrared SEDs

The spatially integrated SED for the dust emission of a galaxy is given by: (43)where is defined by Eq. (39) and is defined by Eq. (41). From Eq. (43) one can see that the free parameters of our model are, to recap:

• the central face-on opacity in the B-band (Eq. (11);Sect. 2.1),

• the clumpiness factor F (Sect. 2.2);

• the star-formation rate SFR (Sect. 2.3.3);

• the normalised luminosity of the old stellar population old (Sect. 2.3.1) and;

• the bulge-to-disk ratio B/D (Sect. 2.3.2).

In cases where detailed modelling of the optical/NIR images is also available, e.g. the study of NGC 891, the old and B/D parameters are independently constrained from the optimisation of the optical images, and thus only 3 free parameters are needed for the model: , SFR, F.

We note that the free parameters , F, and B/D are also free parameters for the attenuation of the UV/optical stellar light given in Paper III (the fourth and final free parameter that affects the attenuation is the inclination i of the disk). Therefore the dust emission SEDs predicted here can be used in conjunction with the predictions for attenuation given in this paper for a self-consistent modelling of the UV/optical-IR/submm SEDs.

The parameter can be related to the total dust mass using the equation: (44)where γ is a constant related to the geometry of the distribution of dust in galaxies and to the dust model. For our model with in parsec, γ = 0.9912   M_⊙   pc^-2.

3. Illustration of the model on NGC 891

We illustrate our model, including description of the intermediate steps in the calculation, using NGC 891. Because of the changes to our model, it is important to see if we can now fit the whole SED of NGC 891, including the MIR emission dominated by the PAH spectral features.

As described in Sect. 2.9, since the luminosity of the old stellar populations has been derived from the optimisation of the optical/NIR images, we only need 3 free parameters to fit the dust emission SED of NGC 891: , SFR, F. For this we ran the calculations for trial combinations of these parameters. We note that the normalised luminosity of the old stellar population derived for NGC 891 is old = 0.792. Throughout this section all the examples shown are for the best fit model¹¹ for NGC 891, which is for , SFR = 2.88 M_⊙/yr and F = 0.41.

3.1. The radiation fields

The first step in the calculation (see Sect. 2.5) is the derivation of the radiation fields illuminating the diffuse dust, which we illustrate in Figs. 4 and 5 through examples of radial and vertical profiles of energy densities u_λ of total radiation fields (not to be confused with profiles of stellar emissivities). The sharp rise in the inner parts of the radial profiles in the K and B band is produced by the dominance of the radiation coming from the bulge. We note here the large variation in the colour of the radiation fields with position, in particular in the radial direction, which, as shown in the next figures, will introduce large differences in the shape of the FIR SEDs. Thus, models that assume radiation fields with the fixed colour of the local interstellar radiation fields (LIRF) are likely to introduce systematic uncertainties in the predictions for the dust emission SEDs.

3.2. The temperature distributions

The next step in the calculations (see Sect. 2.6) is the derivation of the probability distributions of dust temperature. In Fig. 6 we show examples of temperature distributions for grains of different sizes and compositions placed at different positions within the galaxy. The sizes plotted for each composition reflect the range of sizes given by the dust model considered in this paper. For example sizes smaller than 10   Å are only considered for PAH molecules. Also, according to the dust model, the biggest PAHs are those of 100   Å, while graphites have only sizes larger than 100   Å. As expected, small grains exhibit broad probability distributions. For example PAH molecules are only heated stochastically, as they have small sizes. Larger grains have narrower probability distributions, eventually tending towards a delta function, when they start to emit at equilibrium temperature.

The figure also illustrates the dependence of the temperature distributions with position in the galaxy, due to the spatial variation of the intensity and colour of the radiation fields. For example a graphite grain of size 0.0316   μm will exhibit a delta function temperature distribution if placed in the centre of the galaxy (first and third row from the top, middle panels, blue curves in Fig. 6), where the radiation fields have higher energy densities and redder colours, than if placed in the outskirts of the galaxy (second and 4th rows from the top, middle panels). Thus, in the centre of the galaxy the 0.0316   μm graphite grain will emit at equilibrium temperature, while the same grain will start to emit stochastically at R = 15   000   pc.

Finally Fig. 6 shows the dependence of temperature on composition. For example a 0.0316   μm grain placed in the centre of NGC 891, will exhibit a lower temperature if it had a silicate composition (left upper panel) than if it had a graphite composition (middle upper panel), due to the differences in the optical properties of the grains.

3.3. The infrared emissivities

In Fig. 7 we show examples of calculated infrared brightnesses (from Eq. (34)) for grains of different sizes and compositions placed at different positions within the diffuse ISM of our prototype galaxy NGC 891. For the calculations we used the probability distributions of temperature shown in Fig. 6. We also show the same positions and sizes as those from Fig. 6. Following the trends in the temperature distributions, one can see the dependence of the infrared brightnesses with position in the galaxy, due to the variation in the intensity and colour of the radiation fields. One should note that here we show the absolute brightness of each grain, meaning that the emission is not weighted to take into account the abundance of grains according to size or composition; it simply indicates the response of each dust grain to the radiation fields.

By comparing the pairs of positions, first and third row from the top on one hand, and second and 4th row on the other hand, we see that overall the SEDs show a stronger radial variation than a vertical variation and that this trend is independent of grain composition and grain size. This is primarily a result of the fact that the radiation fields show a stronger variation in colour in the radial direction than in the vertical direction (see again Figs. 5 and 4), as described in Sect. 3.1, which, in turn, is a direct consequence of the finite disk plus bulge description of our geometrical distributions of stars and dust. We will discuss here both the change in the peak of the SEDs as well as the change in the overall amplitude of the infrared brightness. As expected, for a given grain size and composition the peak of the SEDs shifts towards longer wavelengths and its amplitude decreases in weaker radiation fields, especially at large galactic radii. However the wavelength shift of this peak and its amplitude as a function of grain size have a more complex behaviour.

i) The wavelength variation of the peak of the infrared brightnesses with grain size for a fixed galaxian position (radiation field).

For the stochastically heated grains there is a strong shift in the peak of the SED towards longer wavelengths with increasing grain size. This is seen for the small grains (Si and PAH; the green curves) placed at large galaxian radii (second and 4th row), which all exhibit strong stochastic heating (as also seen from the broad probability distribution of temperature in Fig. 6, corresponding panels). Since most of the energy is radiated at the highest temperature side of the probability distributions, and since the increase in the grain size will decrease the width of the probability distributions, this means that grains with bigger sizes will reach systematically lower maximum temperatures in the probability distributions, thus radiating at longer wavelengths.

By contrast, grains in equilibrium temperature with the radiation fields will show a small shift in the SED peak with increasing grain size, and in the opposite direction, namely towards shorter wavelengths (see the blue curves from the first and third row from the top, left column of Fig. 7). According to Eq. (34), the peak of the infrared brightness will be determined by the wavelength dependence of Q_abs in the far-infrared and by the equilibrium temperature at which they radiate. For a given radiation field u_rad, the equilibrium temperature depends only on the Q_abs (see Eq. (33)). Since bigger grains are more opaque at shorter wavelengths, tending to a black-body case, they absorb more efficiently, and therefore their equilibrium temperature is higher, providing the Q_abs has the same wavelength dependence in the far-infrared, independent of grain size. This is indeed the case for silicate grains. The graphite grains however have some variation in the wavelength dependence in the far-infrared with grain size, which produce the non-monotonic shift in the peak of their infrared brightness (see first and third row from the top, middle column from Fig. 7, where the blue curves intersect).

A special case is represented by the smallest PAH molecules (the red curves on the right column of Fig. 7). Their infrared brightnesses are dominated by the emission bands due to vibrational transitions, and these features occur at the same wavelengths, independent of the molecule size. This is because the vibrations seen in the mid-infrared correspond to the fundamental stretching or bending modes of the C-H and C-C bonds, and do not involve the molecule as a whole.

ii) The wavelength variation of the peak of the infrared brightnesses with grain size for a variable galaxian position (radiation fields).

The shift in the peak of the infrared brightness towards longer wavelengths with increasing grain size for stochastically heated grains (described at i) is a strong function of galaxian position. By comparing the first and 3rd rows from the top with the second and the 4th rows for the small Si and PAH grains (green curves, Fig. 7) we see that the shift becomes less pronounced for grains at small galaxian radii, where the radiation fields are stronger and redder. For the largest PAH molecules (or corresponding silicate grains) there is almost no shift, their peak remaining constant with wavelength. This shows that these grains, despite being small, are in a transition phase towards equilibrium temperature, as also proven by their narrower temperature distribution from the corresponding panels in Fig. 7. So even PAH molecules can start to emit closer to equilibrium temperature if placed in the centre of the galaxy.

Conversely, the described shift of the peak of the SEDs towards shorter wavelength with increasing grain size for grains heated at equilibrium temperature is also a strong function of galaxian position. By comparing the first and 3rd rows from the top with the second and the 4th rows for the big silicate and graphite grains (blue curves) we see that the shift can be reversed in the weaker and bluer radiation fields at large galactic radii. This means that even big grains can start to emit stochastically in the outer regions of galaxies. Indeed, by looking at the corresponding panels from Fig. 6 we see that the big grains have a non-negligible width in the probability distribution of temperature for the outer disks.

To conclude, the shift of the peak of the infrared brightness as a function of grain size strongly depends on the temperature of the grains, and therefore on the stochastic or non-stochastic nature of the heating mechanism. Since the heating mechanism depends both on grain size and on the intensity and colour of the radiation fields, it is clear that the shift cannot be described in terms of grain size only. This also shows that models that have a fixed grain size for the transition between the main heating mechanisms of dust, irrespective of the radiation fields, will lead to systematic spurious shifts in the mid-infrared to FIR colours with increasing galactocentric radius.

iii) The variation of the amplitude of the infrared brightnesses with grain size for a fixed galaxian position (radiation fields).

As apparent from Fig. 7, the amplitude of the infrared brightness increases with increasing grain size, with a small increase for the small grains and a bigger increase for the bigger grains. Since the width of the SEDs is approximately constant (for a constant set of parameters), the amplitude of the infrared brightness will scale with the area under their SEDs, namely with the energy absorbed (per unit a²). For a given radiation field u_rad, Eq. (33) tells us that this is determined only by the optical properties of the grains in the UV and optical regime. Indeed, overall the Q_abs increases with increasing grain size, with a smaller trend for smaller sizes, and a bigger trend for bigger grains. Especially in the case of silicate, the big grains show a tendency for higher efficiency in absorbing optical photons, which will boost the amplitude of their SEDs due to a higher proportion of red photons being absorbed. Thus, bigger grains will absorb more red photons than the smaller grains, which, for a fixed colour of the radiation fields, will allow big grains to have a higher increase in their infrared brightnesses with increasing grain size.

iv) The variation of the amplitude of the infrared brightnesses with grain size for a variable galaxian position (radiation fields).

As mentioned in iii), the bigger grains have a faster increase in the amplitude of their infrared SEDs with increasing size than the small grains, due to the increase in the efficiency of absorbing red photons. If the radiation fields will also change in colour, due to their spatial variation, this will induce an additional difference in the amplitude of the SEDs of big and small grains. Thus, at small galactocentric radii, where the radiation fields are redder (see Fig. 4), more optical photons are available to boost the amplitude of the SEDs of the big grains. The trend of increasing the relative contribution of big grain emission to the small grain emission with decreasing galactocentric radius is apparent from Fig. 7, especially for the case of silicates, which have a broad range in grain sizes. Obviously the increase in the intensity of the radiation fields with decreasing galactocentric radius will also increase the amplitude of the SEDs, but this will produce an overall boost for both small and big grains.

To conclude, the increase in the relative contribution of big grain emission to the small grain emission is a strong function of the colour of the radiation fields, and, unlike the wavelength dependence, does not depend on the heating mechanism of the grains. This also shows that models that assume a fixed colour of the radiation fields (e.g. that of the local interstellar radiation fields) will incur systematic errors in the mid-infrared to FIR colours with increasing galactocentric radius.

Finally, the increase in the amplitude of the 8   μm PAH features with respect to the 3.3   μm feature at small galactocentric radius is also due to the redder radiation fields, which will provide additional optical photons capable of exciting the vibrational transition at around 8   μm, but not energetic enough to excite the transitions around 3.3   μm.

3.4. The integrated SED

The best fit SED, obtained by spatially integrating the solution described above for the diffuse dust emission component and adding it to the solution for the dust emission from the localised component from the star-formation regions, is shown in Fig. 8, together with the observed FIR/submm data that were used to constrain the model solution. Details on the observed flux densities used in the plot are given in Table 2. The best fit solution corresponds to , SFR = 2.88 M_⊙/yr and F = 0.41. Compared with the solution obtained in Paper I (, SFR = 3.8 M_⊙/yr and F = 0.22) there are changes in all three free parameters. The decrease in is solely attributed to the fact that we used improved observational data to fit the submm points. If we had used the original data the solution for opacity would be unchanged. The decrease in SFR and the increase in F is essentially due to the combination of adding the contribution of the young stellar population in the optical, and of including the wavelength dependence of the fraction of photons escaping the star-forming regions.

With the red line we show the model fit for the total dust emission of NGC 891. We also show the main components of the dust emission, the diffuse component (dark-green line) and the localised emission from the clumpy component (light-green line). One can see that overall most of the dust emission is powered by the diffuse component. Our solution gives for the total dust luminosity  W, of which  W is emitted in the diffuse medium (69%). From the figure it is apparent that the diffuse component dominates the emission longwards of 60 μm and shortwards of 20   μm. Thus, most of the emission in the FIR and in the NIR/MIR (PAH region) is diffuse. It is only at intermediate wavelengths (20−60   μm) where the localised dust emission within the star-forming complexes dominates. In the IRAS 25   μm and Spitzer 24   μm bands almost all the emission is predicted to come from the star-forming complexes, suggesting the efficacy of this band as a direct tracer of SFR. Indeed empirically studies have shown that of the Spitzer bands it is the 24 μm band which is more closely related to star-formation (e.g. Calzetti et al. 2010).

In Fig. 8 we also show the corresponding intrinsic stellar SED (as would be seen in the absence of dust) of NGC 891, together with the stellar emissivity components. In all cases the symbols indicate the wavelengths at which the radiative transfer calculations were performed. The black line represents the total stellar emission produced in the galaxy. In the UV range this is given by the population synthesis models (plotted as the blue line), as explained in Sect. 2.3.3. A fraction of this emission is locally absorbed in the star-forming complexes, while the remaining 1 − F escapes in the diffuse young stellar disk (the mauve line). Thus, the difference between the black line and the mauve line in the UV represents the localised (and wavelength dependent) absorption of stellar light (see Eq. (21)). In the optical range the total emission is given as a sum of the emission from the young stellar disk (mauve line plots only the diffuse emission from the young stellar disk, which is slightly different from the total emission of the young stellar disk, due to the local absorption), the old stellar disk (light blue line) and the bulge (pink line).

One important result of such calculations is to determine the fractional contribution of the different stellar components to the dust heating. For the case of NGC 891 we derive the following fractions: 11% for the bulge, 20% for the old stellar disk, 38% for the young stellar disk and 31% for the star forming complexes. This means that the young stellar populations are responsible for 69% of the dust heating while the old stellar populations account for the remaining 31%.

4. Fidelity of the model

In order to fit observed dust/PAH emission SEDs of real galaxies like NGC 891 it was necessary to make some basic assumptions and approximations, which, however, have implications not only for the SED of the integrated re-radiated light, but also for the geometrical characteristics of this light and of the direct UV/optical stellar light. The main assumptions are the existence of a diffuse dust component associated with the young stellar population, the approximation of this dust in the form of an exponential disk and the utilisation of a fixed spectral emissivity law for the young and old stellar populations. Here we check whether our model is consistent with the available observational constraints beyond those of integrated dust emission SEDs and evaluate the limitations imposed by the approximations.

4.1. The existence of a diffuse dust component associated with the young stellar population

A fundamental aspect of our model was the inclusion of a second disk of dust associated with the young stellar population, which was taken to mimic the diffuse dust that is known to exist in spiral arms through direct observations in the FIR, both in the Milky Way and in nearby well resolved galaxies. The second dust disk was originally introduced in Paper I to provide the observed level of submm emission in the spatially integrated SED. As discussed in Paper I and in Popescu & Tuffs (2005), the additional quantity of dust needed to fit the observed submm data cannot be provided by clumpy optically thick dust, whether such dust is in star formation regions or in passive quiescent clumps. This conclusion has been reinforced by the utilisation of the improved dust emission templates for star formation regions introduced in Sect. 2.8, since these are now empirically constrained by submm and near-mm data.

More fundamentally, although one is in principle free to add dust to the model in any form one likes to reproduce the observed level of emission deep in the submm, if this dust is self-shielded, it will in practice struggle to supply the necessary luminosity to fit the FIR flux density peak of spiral galaxies at around 160   μm. To peak at around 160 micron the dust grains must be heated by radiation fields at around 1 Habing, which is indeed the illumination of the diffuse dust. Heating at around 100 Habing (as at the PDR surface in embedded star-formation regions) would provide the luminosity but with too blue dust FIR colours, while the dust emission from the self-shielded dust would be too red as well as not providing the luminosity¹². Of course, one solution would be to invoke star-forming complexes extending to radii of around 10 times their actual sizes into the diffuse component, so that the dust illumination is reduced to the required levels. However this is akin to spreading the dust around the spiral arms, which is exactly what our second dust disk solution tries to mimic, so the difference to a diffuse dust component then becomes semantic.

In general, the fit to data provided by our two-dust disk model is one of relatively large opacity for the central regions of local universe spiral galaxies (; Paper I; Driver et al. 2007; see also Sect. 3.4), with most of the opacity provided by the second dust disk. This is in contrast to the solution of low opacity () and one disk of dust (obtained by Xilouris et al. 1997; and Bianchi 2007, from modelling the optical data only) which fails by a factor of  ~3 to reproduce the submm emission (Popescu et al. 2000; see also Baes et al. 2010).

Apart from the need to fit the submm data, and the physical considerations already outlined in Sect. 2.1 linking the second dust disk with the corresponding star-forming disk of gas, there is further empirical evidence for the existence of a second dust disk, as derived from comparison of the model predictions with other data, as follows:

i) Comparison of model predictions for the stellar emissivity in the optical with population synthesis models

We checked whether the predicted intrinsic SED of the old stellar populations in NGC 891, as derived from optimisation of the optical data only (Xilouris et al. 1999), together with the corresponding SFR needed to produce the total luminosity emitted by the dust, could be fitted by the population synthesis models of Kotulla et al. (2009). In Fig. 8 one can see that the predictions for the optical emission from the old stellar population (light blue line) in the B, V, and I bands fall severely below the predictions from the population synthesis model (dark blue line). In this plot we already included the second dust disk. In its absence the SFR needed to reproduce the energy emitted in the infrared would be even greater, so the discrepancy would increase. It is clear from here that there is a need for extra stellar luminosity presumably hidden by extra dust (the second dust-disk in our formulation). In Popescu et al. (2000) we assumed that the young stellar population was only emitting in the UV, while in the optical we only had the contribution from the old stellar populations, as derived by Xilouris et al. (1999). This assumption is not supported by the predictions of the population synthesis models, and indicates that the optimisation of the optical images can only reveal information about the old stellar populations and associated dust, but completely hides the information about the young stellar populations and associated dust.

ii) Comparison of model prediction for the spatial distribution of PAH emission with observations

Since our model with the second disk of dust and stellar emissivity was only fitted to the spatially integrated dust emission, it is possible to use the spatial information provided by the observations to check the predictions of the model. We have already used the predictions from Popescu et al. (2000) to compare the predictions of the model for the spatial distribution of infrared emission with the ISOPHOT data of NGC 891 at 170 and 200   μm (Popescu et al. 2004). The comparison was done for the radial profiles, as the observed ISOPHOT images were only resolved in radial direction. The model predictions were found to be in excellent agreement with the observations. More recently Spitzer images of NGC 891 became available, which have a linear resolution of 100 pc at 8   μm, comparable to the thickness of the second dust disk. Since we have now included the PAH features in the model, which dominate the 8   μm emission, these data can be used to test the predictions for the vertical distribution of PAH emission in NGC 891. For this comparison we first subtracted the component of direct stellar light at 8   μm (estimated by multiplying a Spitzer 3.6   μm image by a factor of 0.26; Helou et al. 2004). Figure 9 shows again an excellent agreement. Thus the model containing a second dust disk can predict both the overall level of emission and the shape of the predicted vertical profile.

iii) Comparison of model predictions for the attenuation-inclination relation with observations

A strong independent evidence for the existence of a second dust disk came from the attenuation-inclination relation. It was shown in Popescu & Tuffs (2009) that while a two dust-disk model with higher central face-on opacity can reproduce the observed attenuation-inclination relation, a single dust disk model with completely fails to reproduce the observed data. The attenuation-inclination relation is an especially sensitive test, as the rise in attenuation with inclination will very strongly depend on the relative scaleheights of the assumed dust layers and stellar populations. In particular it is an independent test for the existence of the second component of dust represented by the second dust disk. We should also mention here that one of the strength of this test is that it is completely independent of the assumed dust emission properties. Thus, this would also seem to rule out a one-dust disk model with low opacity but with dust grains having modified optical properties in the submm (e.g. enhanced submm emissivity), as proposed by Alton et al. (2004) and Dasyra et al. (2005).

iv) Comparison of model predictions for the vertical profiles in the optical and NIR in NGC 891 with observations

The solution with a second dust-disk and extra luminosity coming from the young stellar population emitting in the optical bands was further used to predict the appearance of the galaxy in the optical bands. In particular the average vertical profiles have been used to test whether the dip produced in the plane of the galaxy by the dust layer has the right deepness. Figure 10 shows a comparison between the I band averaged vertical profile obtained from the observed images of NGC 891 (Xilouris et al. 1998) and the corresponding model predictions for three cases. The first case (left panel) is the original solution obtained by Xilouris et al. (1999), which only includes the old stellar disk and associated dust. The second case (middle panel) is our two dust-disk model, but without the inclusion of the young stellar disk. Finally, the right panel shows again the predictions for our two-dust disk model, but with the inclusion of extra luminosity coming from the young stellar disk. One can immediately see that the inclusion of a second dust disk which is not accompanied by a stellar luminosity component would produce a stronger dust lane (a larger depth in the vertical profile) than shown in the observed images. This is a clear indication that there is a need for extra luminosity, as also indicated by population synthesis models (see i) above). In fact the two dust disk model with two (old and young) stellar components is able to reproduce the observed data better than the one disk model, especially in the B band.

In the K band the model with the two dust disks and two stellar components predicts a somewhat more prominent dust lane than observed. This would indicate either a shorter scalelength for the second dust disk (though this would be difficult to reconcile with the excellent fit we found in Paper I for the radial profile at 850 μm), or more luminosity in the young stellar disk than predicted by the population synthesis model. It was already noted by Dasyra et al. (2005), Bianchi (2007) and Bianchi (2008) that the second dust disk shows its maximum effect in the K band, since it is only there that the first disk of dust becomes transparent and is therefore not shielding the second disk of dust. A possible alternative reason for the difficulty in fitting the vertical profile in K-band is that, due to its very small scale height, the appearance of the second dust disk at these long wavelengths would be easily blurred in real galaxies if perturbations from co-planarity occur, even if any such perturbations had relatively small amplitudes. Although the scaleheight of the molecular layer in Milky Way is around 90 pc, as adopted in our model for the second dust layer, it is well known that CO and other tracers of star-formation exhibit systematic vertical displacements from the mean place known as “corrugations” (Spicker & Feitzinger 1986, and references therein). The discovery of Matthews & Uson (2008) of a non-planar disk in star-formation tracers in an isolated galaxy other than the Milky Way, is a strong evidence that non-planarity is a rather frequent phenomena. Matthews and Uson found that undulated patterns with amplitude of  ~250 pc are visible in particular in the distribution of the young stellar population and the dust, suggesting that the process leading to the vertical displacements may be linked with the regulation of star formation in galaxies. The effect of corrugations will mean that the edge-on thin disk will not be seen as such in the optical, but only the distribution corresponding to the amplitude of the corrugations. This would be more than enough to blur and hide any dust layer in the NIR. Nevertheless, the young stars will be still perfectly correlated with the thin layer of dust and gas with a (local) scaleheight of 90 pc, since the effect of corrugations is purely gravitational. So the basic vertical geometrical coupling between stars and dust which gives rise to the dust emission from the thin disk, as calculated in this model, still applies. To conclude, if corrugations occur, the second dust-disk will tend to be blurred, making dust lanes in edge-on galaxies in the NIR less prominent.

4.2. Approximating the spiral arm component with an exponential disk

Even if it is accepted that there is a diffuse dust layer associated with the young stellar population, we must still evaluate the effect of artificially distributing this stellar emissivity and dust opacity in an exponential disk instead of a spiral arm pattern. It is clear that this will completely prohibit a comparison of surface brightness distributions in face-on systems, but here we are only concerned with the spatially integrated dust emission SEDs and the effects of this approximation on these.

To do this we ran a simulation with the parameters corresponding to the best fit solution of NGC 891 (see Sect. 3.4), i.e. we kept the same luminosity for the young stellar population and the same amount of dust, but we redistributed the corresponding stellar emissivities and dust opacities in a spiral pattern. We modeled the spiral pattern using 3 circular arms, with a radial distribution as given in Fig. 11. The vertical distribution is the same as for the case of the thin exponential disk. This means that for any given line of sight through the spiral arm the opacity is higher than for the case of a second dust disk, and for any given line of sight through the interarm regions the opacity is lower than for the case of the second dust disk – being just the opacity of the first dust disk. So we will have a solution with a high contrast between arm and interarm regions, that can be still characterised by the same central face-on opacity – which is the effective central face-on opacity if all the dust were distributed in an exponential disk instead of in a spiral arm.

After calculating the radiation fields and the infrared emissivities using the same procedure outlined in Sects. 2 and 3 we obtain a total integrated infrared SED that looks very similar to that obtained for the exponential disk case (see Fig. 12). The integrated dust luminosity is only 5.5% lower than for the standard model. Here we should mention that the attenuation-inclination curve (not plotted in this paper) obtained for the solution with a spiral pattern is also almost identical to that obtained for the case of an exponential disk. The insensitivity of the shape of the attenuation-inclination curve to the inclusion of a spiral pattern was already demonstrated by Semionov et al. (2006,b). In passing, we also note that Misiriotis et al. (2000) showed that the appearance of simulated dusty spiral galaxies seen edge-on, calculated using a spiral structure, does not differ from that calculated using a pure exponential disk.

The similarity in the solutions obtained for both the integrated infrared SEDs and for the attenuation-inclination relation reassures us that the approximation of a second exponential thin disk of stellar emissivity and dust opacity is an excellent one when making predictions for the spatially integrated SEDs. This is in contrast to the large change in the shape and amplitude of the infrared SEDs (as we will show in Sect. 5) and attenuations of stellar light (see Fig. B.1) when changing the main free parameters of our model.

4.3. The approximation of a fixed spectral shape for the SEDs of the old and young stellar populations

As discussed in Sect. 2.3, the dust emission SEDs are calculated by dividing the stellar emission SED into two fixed spectral templates with differing spatial distributions – one UV-dominated corresponding to the young stellar population, weighted by the model parameter SFR, and one optically dominated, corresponding to the old population and weighted by the model parameter old. This is a radically different approach to the handling of stellar emissivity SEDs compared to previous models for the panchromatic UV-submm emission of galaxies. It enables us to obtain the same solution for the dust emission, and therefore the reddening of any given galaxy, without the need to input trial population synthesis solutions for the full UV/optical/NIR stellar emission SED to the calculation of the dust emission. This new approach requires however that the solution to the dust emission is invarient to the assumed shape of the stellar emissivity SED within each of the two templates. Here we test this assumption in turn for the old and the young stellar emission templates.

For the old stellar population the spectral shape was empirically derived from fitting the optical images of NGC 891. The spectral shape was consistent with the SED having a flux density per unit frequency peaking in the J band. Reasonable variations from this shape are to consider spectral templates peaking towards longer wavelengths, e.g. in the K band. We therefore consider a calculation in which the spectral shape was altered to allow for a brightening of the K band luminosity by 58% and a corresponding dimming of the J band luminosity by 48%, change that preserves the spectral integrated luminosity of the old stellar population. The result of this calculation from Fig. 13 shows that the predicted dust emission SEDs are completely insensitive to (possible) changes in the spectral shape of the old stellar populations, as indicated from the overlap of the infrared SEDs.

For the young stellar population the spectral shape was fixed from a combination of population synthesis models and empirical fitting of the optical images. This resulted in a spectral template having a flux density per unit frequency peaking in the I band. Reasonable variations from this shape are to consider spectral templates peaking towards shorter wavelengths, e.g. in the V band. We therefore consider a calculation in which the spectral shape was altered to allow for a brightening of the V band luminosity by 34% and a dimming of the I band luminosity by 25%. As before this change was chosen to preserve the spectral integrated luminosity of the young stellar population. The resulting dust emission SED (not shown in this paper) is essentially indistinguishable from the one calculated using our standard model. There is only an increase in the predicted dust luminosity by 0.2%, which is to be expected due to the shift in the spectral peak of the young stellar population to shorter wavelengths, which resulted in a larger fraction of stellar photons being absorbed by dust. So changing the colours of the young stellar population in the optical by as much as 59% produces changes of less of a percent for the predicted dust emission luminosity.

Finally, it is important to check what the effect of changing the colours of the young stellar population in the UV is, since there it is these photons that make a significant contribution to the dust heating. We therefore made a more drastic change, by changing the slope of the UV SED from nearly flat (in F_ν) to a monotonically decreasing slope between the 1350   Å and 2800   Å spectral sampling points. We therefore brightened the 1350   Å data point by 24% and dimmed the 2800   Å flux by 38%, thus producing a very blue spectrum. As before we kept the overall luminosity of the young stellar population constant. The resulting dust emission SED is shown in the two panels of Fig. 14, plotted in a linear scale and for different cuts in luminosity to allow a better visualisation of the small differences in the SEDs. Again, since a bluer stellar SED was considered, this resulted in more stellar photons being absorbed by dust and a larger predicted dust luminosity, as seen from the plots. The overall increase in the dust luminosity was 1.7%, which is still a minor variation taking into account the dramatic change by 62% in the colours of the UV stellar SED. In fact this is almost surprising, but we should keep in mind that the bluer UV photons have higher probability of absorption in the star forming complexes, and therefore the escaping radiation illuminating the diffuse component would still be strongly reddened due to the local absorption.

By contrast we also did a calculation in which we changed the relative contribution of the young to old stellar populations and kept fixed the total stellar luminosity. We have used the solution for NGC 891 for the diffuse component, set the bulge-to-disk ratio to 0 and used this as a reference SED. We then increased the luminosity of the young stellar population by a factor of 2 and decreased the luminosity of the old stellar populations in the disk by a factor of 1.68, which preserves the total stellar luminosity. The resulting SEDs are shown in Fig. 15. This time there is a significant change in the SEDs, followed by an increase in the total dust luminosity by 33% and an increase in the MIR emission by a factor of  ~2. This shows again that the main factors influencing the dust emission SEDs are the overall luminosities of the old and young stellar populations, and not (possible) variations in the spectral shape of the template stellar SEDs.

4.4. The relative opacity of the first and second dust disk

A further approximation of the model is the fixed ratio between the opacity of the first and second dust disks. While most of the fixed parameters have been calibrated to empirical relations, we acknowledge that this ratio has been only calibrated to the value of NGC 891 (close to NGC 5907). The validity of this assumption has been succesfully tested in a statistical sense on the attenuation-inclination relation in Driver et al. (2007). However this assumption still needs to be proven on an object-by-object basis by fitting the panchromatic SEDs of galaxies. As we will show in Appendix D we have identified a parametric test to potentially flag out galaxies that may not follow the colour-luminosity relation predicted by our model due to other geometrical characteristics.

5. The library of SEDs for the diffuse component

We have created a library of diffuse SEDs that spans the parameter space of 4 parameters: , SFR′, old, and B/D¹³. With the exception of SFR′, these are the main parameters of the model. SFR′ is defined in terms of the primary parameters SFR and F: (45)which is equivalent with the SFR for the case that F = 0. This definition is possible because in our model we made the assumption that the wavelength dependence of the fraction of escape photons from the clumpy component into the diffuse one is fixed, allowing us to separately calculate the SEDs for the diffuse component and for the clumpy component. In this formulation SFR′ is an effective star formation rate powering the diffuse dust emission.

In Table 3 we give the parameter values at which the model for the diffuse component has been sampled. In total we have 7 (for (for SFR′) × 9 (for old) × 5 (for B/D) = 2835 combinations, however those with SFR = 0 and old = 0 (no stellar luminosity at all) cannot exist, therefore the library contains only 2800 simulated SEDs. All the model SEDs are available in electronic form. The choice of the parameter values was done such that it covers the parameter space of local universe galaxies but also the asymptotic values of these parameters. Thus for we considered the range [0.1, 8.], which means galaxies with almost no dust at all to galaxies having over twice as much opacity as the average value found from the statistical analysis of the Millennium Galaxy Catalogue (Driver et al. 2007). For SFR′ we considered the range [0., 20.] M_⊙/yr, which means galaxies with no recent star-formation activity (no heating from the young stellar population) to galaxies having very high star-formation rate (typical of starburst galaxies which should be outside the range of spiral galaxies addressed by this paper). For the old stellar population the parameter old was scanned analogous to SFR′. Finally, we considered galaxies spanning the whole Hubble sequence, from bulgeless galaxies (B/D = 0.0) to bulge dominated galaxies (B/D = 2.0). This choice of parameters allows a smooth overlap with the parameter space of starburst galaxies, very quiescent galaxies, early type galaxies, and higher redshift galaxies. It also allows rare (unexpected) cases to be considered. As we will show in Sect. 6.5, the locus of the FIR colours corresponding to the parameter space defined by the values in Table 3 overlaps quite well with the locus of observed FIR colours of real life spiral galaxies.

To compute the library of diffuse SEDs we first calculated a library of radiation fields, computed for the main stellar components of our model: young stellar disk, old stellar disk and bulge, for the 7 values of opacity used in Table 3 and for the 16 UV-optical wavelengths detailed in Sect. 2.3. For the bulge and the old stellar disk only 6 optical wavelengths were used, as in our model we assume that these stellar components have negligible emission in the UV range. In total we created a library of 196 data cubes of radiation fields, sampled at 22 radial positions and 12 vertical positions (264 spatial points within the model galaxy). We then created the library of temperature distributions, for the 4 grain compositions used in our model (silicate, graphite, PAH⁰, PAH⁺), the 7 values of opacity, 9 values of SFR′, 9 values of old and 5 values of B/D. In total there are 11 340 combinations, however those with SFR = 0 and old = 0 cannot exist. Therefore the library contains only 11 200 data cubes of temperature distributions, each sampled at 264 spatial points within the model galaxy and for all grain sizes contained in the dust model. For each of the temperature distribution data cubes we created corresponding cubes of infrared emissivity. In total a library of 11 200 files of infrared emissivities were calculated.

When used to fit observed panchromatic SEDs, the library of simulated dust and PAH emission SEDs should be used in conjunction with the library of simulated attenuations of stellar light recalculated in this paper, taking into account the B/D ratio of the galaxy and its inclination, as described in Appendix D.

6. Predicted variation of the dust emission SEDs with the main parameters of the model

6.1. Variation of the SEDs with

i) Amplitude of the SED

As expected, the amplitude of the diffuse dust and PAH emission SEDs increases with increasing optical depth for the optically thin cases and tends to a saturation value for the optically thick cases. This is seen in Fig. 16 from the bunching of the blue and black curves on one hand () and from the big gap between the green and red curves on the other hand (). An interesting feature of the SEDs is the fact that the ratio between the FIR and MIR (PAH) amplitudes increases with increasing opacity for the models where the stellar luminosity has a higher contribution from the old stellar population with respect to the young stellar population (but where the young stellar population has still a non-negligible contribution). This can be seen for example on the right bottom panel of Fig. 16, where the green curve () shows almost identical levels in the FIR peak and in the PAH features, while the black curve () shows two order of magnitude difference in the FIR to MIR levels. This change in the FIR to MIR colour is due to the fact that with increasing optical depth the disk becomes first optically thick to the UV radiation (provided by the young stellar population), while still being relatively transparent to the optical photons (mainly provided by the old stellar population). This will have the consequence that the PAH emission (and small grain emission), which is mainly heated by the UV photons will tend to a saturation level, as expected for the optically thick case, while the FIR emission will be still boosted by the optical photons, which have not reached the optically thick limit yet. Apart from this effect, which is simply related to the optical properties of the grains, there is an additional geometrical effect which will boost the FIR-to-MIR colour with increasing opacity. This is due to the fact that the old stellar populations have larger scale-heights than the young stellar populations, which will mean that a larger proportion of the optical photons will be less confined to the regions of higher optical depth, and will therefore provide extra heating in the optically thin regions of the galaxy.

ii) Peak of the SED

Another interesting feature of the plots is the fact that the peak of the SEDs from diffuse dust shifts towards longer wavelengths with increasing opacity for the optically thick regime, for the models where the stellar luminosity has a higher contribution from the young stellar populations with respect to the old stellar population. This can be seen for example in the bottom row of Fig. 16, where the trend of shifting the peak of the SEDs becomes increasingly less strong in moving from the left to the right panel, where the solution changes from a young to an old stellar population dominance of the stellar SED. As before, when moving to a solution where more optical photons are available, these will provide extra heating to the dust, due to their more optically thin regime, both because of their optical properties and because of their spatial distribution (as discussed before), which in turn will cancel the shift towards cooler SEDs.

6.2. Variation of the SEDs with SFR′ (young stellar population)

The increase in the contribution of the young stellar population to the stellar SEDs will have the effect of increasing the MIR to FIR level of the diffuse dust emission SEDs, with essentially no change in the submm level (see Fig. 17). At the same time the peak of the dust emission SEDs will shift towards shorter wavelengths. Overall this will result in warmer SEDs, with both FIR SED peaks and MIR-to-FIR colours becoming systematically bluer. It is interesting to note that the effect of increasing SFR is completely different from the effect of increasing , showing that the two parameters are completely orthogonal. One should also notice that the trend in bluer SEDs with increasing SFR becomes less pronounced for models with higher contribution from the old stellar populations to the stellar SEDs (see trends in moving from the left column to the right column of Fig. 17).

6.3. Variation of the SEDs with old (the old stellar population)

The increase in the luminosity of the old stellar population produces a shift of the peak of the infrared SED of the diffuse dust towards shorter wavelengths (see Fig. 18), similar to the case of increasing the luminosity of the young stellar population. However, the shift is accompanied by a decrease in the MIR-to-FIR ratio, opposite to the effect obtained in the case of increasing the luminosity of the young stellar population. Indeed, the MIR-to-FIR colours become cooler, due to the fact that the additional optical photons will mainly boost the FIR emission and not the PAH emission. Obviously this effect is less pronounced for models with higher luminosities of the young stellar populations (see trends from going from the left to the right hand column in Fig. 18) and is enhanced for models with higher dust opacity (see trends for going from the top to the bottom row in Fig. 18).

6.4. Variation of the SEDs with B/D (bulge-to-disk ratio)

The variation in the bulge-to-disk ratio produces variations in the infrared SEDs with less dynamical range (see Fig. 19), simply due to the smaller dynamical range in the values of the parameter B/D. Nevertheless, the peak of the SEDs is shifted to shorter infrared wavelengths with increasing B/D and the MIR-to-FIR colours get cooler, following the trends expected for a variation in the luminosity of an old stellar population. Thus the variation induced by the B/D ratio follows similar trends with the variation induced by the old parameter. It is however expected that in most cases the value of the B/D to be inputted from observations, and thus not to be a free parameter of the model.

Overall it is interesting to observe that the 3 parameters of the model: , SFR and old are orthogonal parameters, producing quite different effects in shaping the dust emission SEDs. For example it is clear that the only way to increase the level of the submm emission is through a variation of , as neither an increase in the luminosity of the old or of the young stellar populations could account for this. In other words in the submm we are tracing dust column densities, while in the MIR and FIR we are tracing both dust column densities and heating sources. It is also obvious that if we want to have warmer MIR-to-FIR colours we have to increase the luminosity of the young stellar populations, while if we want to have cooler MIR-to-FIR colours we need to increase the luminosity of the old stellar populations, with different degrees of modulations due to variations in . Here we need to remember that for fitting total SEDs the template for the star-forming complexes must be added to the diffuse SEDs, boosting the MIR emission (see Sect. 6.5, below).

6.5. Colour variation of the combined diffuse and localised dust emission SEDs

Adding the emission from the localised dust emission template of star-formation complexes to the solution for the diffuse emission to obtain the total emission adds extra variations to the dust emission SEDs, especially boosting the warm dust emission from larger grains. This effect, controlled by the F parameter, is in practice rather orthogonal to the variations in the SED of the diffuse component discussed above, due to the different behaviour of the UV/optical attenuation as a function of F (see Appendix C) and due to the fact that, over the parameter range of the model, the localised dust emission template is almost always much warmer in the FIR than the predicted FIR colours of the diffuse dust emission. In fact the FIR colours of the diffuse dust emission only approach those of the PDR template for the extreme case SFR′ = 20 and old = 20 in an optically thin galaxy in which the UV illuminates the diffuse dust, as can be seen by comparing the top RH panel of Fig. 18 with the template SED in Fig. 5. A further reason for the orthogonality of the template SED to the SED of the diffuse emission is that dust in the diffuse ISM is more efficient in chanelling absorbed UV energy into the PAH features than dust in the star-formation regions, due to the photo-destruction of PAH in star-formation regions (see Sect. 2.8). Overall, solutions for the total infrared emission require most (though not all) of the emission in the 20−60   μm range to be localised emission from star-formation regions, with the MIR and FIR/submm emission flanking this range being predominantly diffuse in origin.

The combined effect of the 5 main parameters of the model , SFR, old, B/D and F will produce a large variation in the colours of the simulated SEDs. In Fig. 20 we plot the 170/100   μm colour versus the 100/60   μm colour for our model SEDs. The SEDs were derived by combining all the simulated SEDs of the diffuse component from our library with the template SED for star-forming regions, whereas the combination was done for different values of the F factor (including the asymptotic values). This defines the locus in the colour-colour diagram occupied by our model SEDs. To check that the parameter space covered by the models overlaps with the parameter space of real life spiral galaxies, we overplotted the corresponding colours of the spiral galaxies with Hubble type earlier than Sd from the ISOPHOT Virgo Cluster Deep Survey (IVCDS; Tuffs et al. 2002a,b; Popescu et al. 2002) and from the catalog of compact sources of the ISOPHOT Serendipity Survey (Stickel et al. 2000). One can see that except from a few outliers, the observed galaxies lie on the locus defined by our models. One can also observe that the position of NGC 891 in this diagram does not have any preferential place, but lies somewhere towards the very quiescent region of the models (as expected). More detailed comparisons with data and applications of the model will be given in a series of future papers.

From Fig. 20 it is also obvious that our models embrace a somewhat larger area in parameter space than that defined by real life spiral galaxies. This was done on purpose, from the choice of the minimum and maximum values of the main parameters of the model. We tried to push the limited values towards the asymptotic values, in order to cover even rare (unexpected) cases that could occur in real life (see Sect. 5 and Table 3 for the range of parameter space).

Most fundamentally, however, we emphasise that in the context of a RT solution it is the combination of the observed amplitude and colour of the dust emission SED with the self-consistently applied constraint on the attenuation in the UV/optical that leads to an inherently non-degenerate solution.

7. Discussion

In this paper we have assembled the components from which we can combine the predictions for the attenuation of starlight (as originally formulated in Paper III and re-expressed here in Appendix B and Appendix C) with the predictions for the SED of the dust/PAH re-radiated starlight (self-consistently calculated as described in Sect. 2) to extract the physical parameters of a star-forming galaxy from an observed UV/optical/NIR – mid-IR/far-IR/submm SED.

The detailed mathematical prescription of how to combine the library of dust emission SEDs with the library of attenuation to fit UV/optical to infrared/submm SEDs is given in Appendix D. Here we discuss some of the issues regarding the applicability and physical relevance of our model.

Firstly, our model is designed for data sets that have good enough optical images to derive B/D and i parameters. For the i parameter this should be the case for most optical surveys of local universe galaxies. For example, the SDSS, with an angular resolution of 1.6′′, was able to derive inclination angles for galaxies with recession velocities ranging to  ≈ 10   000 km s^-1. With modern surveys routinely having  ≈0.5′′ resolution, inclinations should be available up to z ≈ 0.1, which encompasses the local Universe galaxies. In cases where inclination angles cannot be derived, an expectation value may be assumed instead, and the quartile range of uncertainty in derived physical parameters resulting from a lack of knowledge of i can be derived by running separate optimisations for the quartile values of i. It should be noted that the inclination angle should not be considered as a free parameter, as its effect on the SED would not be orthogonal to the parameter. The other parameter B/D should also be an input parameter for the routine and not a free parameter. In cases where information on B/D is missing, users are advised to use the correlation between B/D and apparent optical-NIR colours to estimate the B/D.

The use in the model of the third parameter accessible through optical observations, θ_gal, is qualitatively different to that of B/D and i. Whereas B/D and i can (for practical applications) only be used as constraints, θ_gal may in practice be used either as a constraint (as in the example in Appendix D) or as a free parameter. When θ_gal is used as a free parameter, the angular surface area effectively serves as free amplitude scaling factor for the amplitude of the diffuse dust emission (see Eq. (D.5)).

Secondly, our model is in its present form only designed to analyse star-forming galaxies. As noted above, when performing an optimisation with the optical structural parameters fixed by direct observation, the optimisation no longer has any degree of freedom in terms of scaling parameters. This makes the model sensitive to a test of the fundamental assumption, that the dust is exclusively heated by stellar photons. The fits will be systematically biased if dust is heated by photons from non-stellar sources like AGN or by some completely different channel. For example, by incorporating the optical and UV constraints, the model best fit parameters for SFR and will be likely to be skewed upwards if an AGN heating dust with a harder UV photon spectrum than expected from a pure stellar source is present. In this case the dereddened UV/optical spectrum from step vii of Appendix D will however provide a direct flag for the presence of a bluer spectral emissivity in the UV than could be provided by stars. Potentially, therefore, our SED fitting technique can provide a method for recognising the presence of dust-obscured AGN activity.

Thirdly, we note that the spectral form of the diffuse component of the dust emission SED given by Eq. (D.5) predicts that galaxies will lie on specific locii in a colour – surface brightness diagram for the dust emission, dependent on the surface density of UV/optical emissivity and . In cases where θ_gal is known from optical measurements, any deviation of the observed positions of galaxies in colour – surface brightness space from the positions predicted by the model can be therefore used as a non-parametric test of the fidelity of the geometry used in the model without having very specific a priori knowledge of the physical properties of the observed sources. This would also serve as the most direct verification of the major role played by emission from diffusely distributed dust in our model.

Finally, we caution the potential user of this model that, although RT solutions applied in this way are in principle a very powerful way to constrain physical properties of galaxies, their predictive power is ultimately reliant on an priori knowledge of the geometry of the system. Although we have specified this geometry using empirical constraints derived from spiral galaxies of intermediate morphological type, and have endeavoured to extend the applicability to galaxies of earlier or later types through the B/D parameter, constrained from optical imaging, the practical range of morphological types where the advantages of the RT treatment outweigh deviations of real from assumed geometry still remain to be determined. This question can only be answered by statistical analysis, such as potentially realisable in a non-parametric way via the colour-surface-brightness relation for the diffuse dust emission component (as proposed above) for different morphological classes of galaxies, and will be the subject of future studies.

8. Summary

In this paper we presented a radiative transfer model for spiral galaxies that self-consistently accounts for the attenuation of stellar light in the UV/optical and for the dust re-emission in the MIR/FIR/submm. We used this model to create a comprehensive library of dust and PAH emission SEDs and corresponding attenuations of stellar light (originally presented and described in Paper  III and given in revised form in this paper) that can be used to routinely fit the panchromatic SEDs of large statistical samples of spiral galaxies.

Our model has been calibrated for local Universe galaxies, and therefore its applicability should mainly lie within low redshift galaxies. The models are also targeted to spiral galaxies. We have not attempted to model elliptical galaxies, starburst galaxies, dwarf galaxies or AGN nuclei.

The free parameters of our model for the calculation of the infrared SEDs are:

• the central face-on opacity in the B-band ;

• the clumpiness factor F which defines the fraction of UV photons locally absorbed in star-formation regions and the corresponding fraction available to illuminate the diffuse dust;

• the star-formation rate SFR;

• the normalised luminosity of the old stellar population old;

• the bulge-to-disk ratio B/D.

For the parameter we have provided a relation for converting to total dust masses (see Eq. (44)).

The free parameters , F, and B/D are common to the set of parameters needed to predict the attenuation of the UV/optical stellar light (the fourth free parameter that affects the attenuation is the inclination i of the disk). Therefore the simulated dust emission SEDs presented here can be used in conjunction with the predictions for attenuation given in this paper for a self-consistent modelling of the UV/optical-IR/submm SEDs. Our model has also a parameter the angular size θ_gal, which can be used either as a constraint (fixed via optical images) or as a free parameter. When θ_gal is used as a free parameter it determines the amplitude scaling of the fit.

Since in our model we made the assumption that the wavelength dependence of the fraction of escape photons from the clumpy component into the diffuse one is fixed, we can separately calculate the SEDs for the diffuse component and for the clumpy component. For the clumpy component we adopted the model of Groves et al. (2008), which we tuned to fit the total MIR/FIR/submm SED of the ensemble of star-forming complexes in the Milky Way (Sect. 2.8). For the diffuse component we defined an effective star-formation rate SFR′ (corresponding to the total illumination by the young stellar population of the diffuse dust) and created a library of diffuse SEDs that spans the parameter space of , SFR′, old, and B/D. In total we sampled 7 values in , 9 in SFR′, 9 in old and 5 in B/D, making a total of 2835 combinations. However those with SFR = 0 and old = 0 (no stellar luminosity at all) cannot exist, therefore the library contains only 2800 simulated SEDs. This corresponds to a library of 196 data cubes of radiation fields, sampled at 22 radial positions and 12 vertical positions (264 spatial points within the model galaxy), a library of 11 200 data cubes of temperature distributions and a library of 11 200 files of infrared emissivities. We described (Appendix D) how in practice one can combine the predictions for the attenuation of starlight with the predictions for the SED of the dust/PAH re-radiated starlight to extract the physical parameters of a star-forming galaxy from an observed UV/optical/NIR – mid-IR/far-IR/submm SED.

The analysis of the library of diffuse integrated SEDs (Sect. 6) showed that the main parameters of the model are quite orthogonal, producing different effects in shaping the dust emission SEDs. Specifically we have shown that:

• The amplitude of the dust and PAH emission SEDs increaseswith increasing optical depth for the optically thin cases and tendsto a saturation value for the optically thick cases.

• The ratio between the FIR and MIR (PAH) amplitudes increases with increasing opacity for the models where the stellar luminosity has a higher contribution from the old stellar population with respect to the young stellar population.

• The peak of the SEDs shifts towards longer wavelengths with increasing opacity for the optical thick regime, for the models where the stellar luminosity has a higher contribution from the young stellar populations with respect to the old stellar population.

• The increase in the contribution of the young stellar population to the stellar SEDs will have the effect of producing warmer SEDs, with both FIR SED peaks and MIR-to-FIR colours becoming systematically bluer.

• The increase in the luminosity of the old stellar population produces SEDs with warmer FIR SED peaks but with redder MIR-to-FIR colours.

• The variation in the bulge-to-disk ratio produces variations in the infrared SEDs with less dynamical range, simply due to the smaller dynamical range in the values of the parameter B/D. Overall the peak of the SEDs is shifted to shorter infrared wavelengths with increasing B/D and the MIR-to-FIR colours get cooler, following the trends expected for a variation in the luminosity of an old stellar population.

Overall we found that the effect of increasing SFR is completely different from the effect of increasing and works in the opposite direction from that produced by an increase in old. We also verified that the only way to markedly increase the level of the submm emission is through a variation of .

The spatially integrated SEDs were derived from simulated images of dust emission, themselves calculated from 3D data cubes of infrared emissivity representing the response of grains (integrated over size distribution and composition) to the radiation fields derived from radiative transfer calculations. Our calculations take full account of the large variations in colour of the radiation fields as a function of position in the galaxy, which we show to be important (Sect. 3.1). Such variations lead to large trends in the grain temperature distributions (Sect. 3.2) used in the calculation of stochastic emission, leading to corresponding trends in the infrared brightnesses (Sect. 3.3) with position in the galaxy. In particular we concluded that:

• i)

  The shift of the peak of the infrared brightness as a function ofgrain size strongly depends on the temperature of the grains, andtherefore on the stochastic or non-stochastic nature of the heatingmechanism. Since the heating mechanism depends both on grainsize and on the intensity and colour of the radiation fields, it is clearthat the shift cannot be described in terms of grain size only. Thisalso shows that models that have a fixed grain size for the transitionbetween the main heating mechanisms of dust, irrespective of theradiation fields, will lead to systematic spurious shifts in the mid-infrared to FIR colours with increasing galactocentric radius.

• ii)

  The increase in the relative contribution of big grain emission to the small grain emission is a strong function of the colour of the radiation fields, and, unlike the wavelength dependence, does not depend on the heating mechanism of the grains. This also shows that models that assume a fixed colour of the radiation fields (e.g. that of the local interstellar radiation fields) will incur systematic errors in the mid-infrared to FIR colours with increasing galactocentric radius.

Basis for the radiative transfer calculations were empirical constraints for the geometry of the large scale distributions of stellar emissivity and dust opacity of the translucent and opaque components of galaxies (Sect. 2.1). The translucent components were constrained using the results from the radiation transfer analysis of Xilouris et al., while the geometry of the optically thick components is constrained from physical considerations with a posteriori checks of the model predictions with observational data. The main checks concern the nature of the diffuse distribution of dust associated with the young stellar population which is needed to predict the observed FIR/submm emission of real life galaxies. These checks (Sect. 4) showed that:

• i)

  Comparison of model predictions for the stellar emissivity inthe optical with population synthesis models is not consistentwith one dust disk model with one stellar disk component, andindicates the need for extra luminosity hidden by extra dust – thesecond dust disk and young stellar disk in our model formulation.

• ii)

  Comparison of model predictions for the vertical profiles of PAH emission with observations at 100 pc linear resolution reveals a good agreement, consistent with the existence of a thin dust layer, as represented by the second dust disk.

• iii)

  Comparison of model predictions for the attenuation-inclination relation with observations is also consistent with the existence of a thin dust layer, as represented by the second dust disk.

• iv)

  Comparison of model predictions for the optical surface brightness distributions of the edge-on galaxy NGC 891 with observations indicates that a two dust disk model with two (old and young) stellar components is able to reproduce the observed data better than the one disk model with one stellar disk. We also showed that a two dust disk model which is not accompanied by an extra stellar luminosity component produces a stronger dust lane than observed in real images. In the K band the model with the two dust disks and two stellar components predicts a somewhat more prominent dust lane than observed, indicating either a shorter scalelength for the dust disk (though this would be difficult to reconcile with the excellent fit we found in Paper I for the radial profile at 850 μm) or more luminosity in the young stellar disk than predicted by the population synthesis models. A possible alternative reason for the difficulty in fitting the vertical profile in K-band is that, due to its very small scaleheight, the appearance of the second dust disk in the K band, where the first dust disk becomes transparent, would be easily blurred in real galaxies if perturbations from co-planarity occur, as has been observed in the form of “corrugations” in the Milky Way and external galaxies.

We also investigated the approximation of the spiral arm component with an exponential disk (Sect. 4.2), and showed that this approximation has negligible effect on the predicted integrated dust emission SEDs when the amount of dust is preserved.

A particular feature of our model is the representation of the stellar emissivity SEDs in terms of just two parameters, SFR and old, since this enables UV/optical SEDs to be dereddened without recourse to population synthesis models, thereby avoiding bias due to the age/metallicity-opacity degeneracy. This feature entailed the approximation of a fixed spectral shape for the SEDs of the old and young stellar populations. We quantified this approximation showing it has a negligible effect on the dust emission SEDs (Sect. 4.3).

We have derived updated and improved solutions for the edge-on galaxy NGC 891 (Sect. 3), one of the galaxies from the set used to constrain the geometry of the model. Since the luminosity of the old stellar populations has been fixed to the values derived by Xilouris et al. (1999) from the optimisation of the optical/NIR images, old = 0.792 and B/D = 0.33, only 3 free parameters were needed to fit the dust emission SED of NGC 891, namely , SFR, F. For the integrated SED of NGC 891 the best fit solution was found for , SFR = 2.88 M_⊙/yr and F = 0.41, which corresponds to a total dust luminosity of  W, of which  W is emitted in the diffuse medium (69%). The results indicate that the diffuse component dominates the emission longwards of 60 μm and shortwards of 20   μm, while the localised dust emission within the star-forming complexes dominates at intermediate wavelengths (20−60   μm). It is only at intermediate wavelengths (20−60   μm) where the localised dust emission within the star-forming complexes dominates. The total dust emission of NGC 891 is predominantly powered by the young stellar populations, which contribute 69% to the dust heating. The detailed input from the different stellar components is as follows: 11% from the bulge, 20% from the old stellar disk, 38% from the young stellar disk and 31% from the star forming complexes.

Online material

Appendix A: Calculation of the extinction cross-section and scattering phase function

A prerequisite for the calculation of radiation fields and infrared emission is knowledge of the extinction cross-section of the population of grains, and, since we consider anisotropic scattering in the radiation transfer calculations, the scattering phase-function needs to be known as well.

The absorption cross-section of grains of composition i =  { Si, Gra, PAH⁰, PAH ⁺  } , C_abs,   i, is obtained by integrating the absorption efficiencies Q_abs,   i over the grain size distribution n(a):

(A.1)where C_abs,   i is given in units of  [cm²   H^-1] , a_min is the minimum grain size and a_max is the maximum grain size.

Similarly, the scattering cross-section of grains of composition i, C_sca,   i, is obtained by integrating the scattering efficiencies Q_sca,   i over the grain size distribution n(a): (A.2)Then, by summing over the grain composition i we obtain the total absorption and scattering cross-sections, C_abs and C_sca: The extinction cross-section C_ext is the sum of the absorption and scattering cross-sections: (A.5)We note that the extinction cross section C_ext is defined as per unit H. In some applications it is useful to define a cross section per unit dust mass, which we denote here : (A.6)where ρ_g is the density of the grain material and is in units of cm²/g.

In turn, is related to the extinction coefficient κ_ext, as used in the mathematical prescription of the dust distribution from Eqs. (4) and (6) using: (A.7)where ρ_dust(R,z) is the dust mass density at position (R,z) in the galaxy in units of g   cm^-3 and κ_ext(λ,R,z) is in units of cm^-1.

Figures A.1 and A.2 show the resulting extinction curve of the dust model adopted here, together with the absorption and scattering components (Fig. A.1) and the components given by the different grain composition (Fig. A.2). As expected, the figures confirm that the model extinction curve fits well the observed mean extinction curve of our galaxy.

The averaged anisotropy of the scattering phase function g needed in the radiative transfer calculation is obtained in a similar manner to Eqs. (A.1)–(A.4). (A.8)(A.9)where Q_{phase,   i} is the anisotropy efficiency.

Appendix B: The library of attenuations of stellar light for the diffuse stellar components

The second set of simulated data needed to fit the panchromatic SEDs is the library of attenuations in the UV/optical/NIR as a function of and i. As mentioned in Sect. 2.4, the revision of the dust model required a recalculation of the database for the attenuation of stellar light, as presented in Paper III. The overall concept and characteristics of the calculations are the same as in Paper III, but a small change, mainly in the zero point of the calculations, was apparent due to the change in the relative contribution of absorption and scattering to the total extinction. Here we only give an example of a comparison between the attenuation-inclination curves obtained using the new and the old dust model (Fig. B.1).

The attenuation-inclination relations for disks show a systematic change with inclination when changing the dust model. Thus the attenuation for the low inclinations is decreased more than for the high inclinations, with a tendency for the curves to converge at the edge-on inclinations. This means that the shape of the attenuation-inclination is steepened for the present dust model. The curves for bulges show the biggest offset when changing the dust model, but in most of the cases there is no change in the shape of the curves. As one can see the curves run almost parallel, except perhaps for the lowest values of opacity. The smaller change is seen for the thin disk component, where neither the shape nor the zero point are changed significantly.

We also did some tests to quantify the effect of the change in the dust model to the overall energy balance. By integrating the attenuation over all angles we obtained an estimate of the total energy absorbed by dust in a galaxy. This absorbed energy was found to be on average 10% smaller for the attenuations calculated using the new dust model than for those from Paper III.

Appendix C: Formulation of composite attenuation of stellar light

In Sect. 5.1 of Paper III a generalised formula was given, showing how the composite attenuation (that is the overall attenuation of an arbitrary combination of luminosity components from stellar populations in the young stellar disk, the old stellar disk and the bulge) can be derived from the library of attenuation of stellar light from the diffuse component and the attenuation of the clumpy component. At any wavelength, the composite attenuation depends on the relative luminosity of the three stellar components, which we have described in this paper in terms of the parameters SFR, F, old and B/D, which we used to describe the dust emission. Here we re-write the generalised expression for the composite attenuation (Eq. (16) from Paper III) for the specific parameterisation adopted in this paper.

At a given wavelength λ, the composite attenuation Δm_λ in a galaxy is given by:

(C.1)where and L_λ are the intrinsic and the apparent luminosity densities. The quantities and L_λ can be further expressed as a summation of the corresponding quantities for the disk, thin disk and bulge: The apparent and intrinsic luminosity densities for the disk, thin disk and bulge are related as follows:

where , and are the attenuation values expressed in magnitudes for the diffuse component in the disk, thin disk and bulge and , and are the corresponding attenuations expressed in linear form. Using Eqs. (13), (15), (19) and Eqs. (C.4) − (C.6) we can rewrite Eqs. (C.2) and (C.3) as: By making the notation: (C.9)the composite attenuation from Eq. (C.1) becomes: (C.10)where (C.11)For the case that old = 0 in the optical and in the UV Eq. (C.10) becomes: (C.12)The use of the calibration factor F_cal in our procedure means that in practice the equations describing the attenuation due to the diffuse dust illuminated by the young stellar disk need to be rescaled to accommodate different values of F than those used in the calibration. For this we need to use the correction factor for the diffuse component corr^d(F) as defined in Eq. (29), and rescale Eq. (C.5) in a similar way to the formulation of the radiation fields in Sect. 2.5.2: (C.13)In this case Eq. (C.11) becomes: (C.14)and Eq. (C.12) becomes: (C.15)Equation (C.14), together with Eqs. (C.10) and (C.11) are the analog of the original expression for the composite attenuation from Eq. (16) in Paper III. Equation (C.15) together with Eq. (C.12) are the analog of Eq. (17) from Paper III.

Appendix D: How to use the model to fit UV/optical to infrared/submm SEDs

The six physical parameters determining our model prediction of the SED are:

• , SFR, F, old, B/D and i¹⁴

whereby

• , SFR, F, old, and B/D

determine , the model prediction for the dust luminosity density in the IR/submm as given by Eq. (43) of this paper, whereas

• , SFR, F, old, B/D and i

determine the attenuation in the measured UV/optical/NIR emission (whereby the dependence on SFR and old is a weaker dependence due to the dependence of the composite attenuation on the relative amplitudes of the young and old stellar populations, as described in Appendix C). This attenuation is given in magnitude form, Δm_λ, by Eqs. (C.10), (C.11) and (C.14). In the following it is convenient to express it in the linear form: (D.1)We note that two of the six physical parameters, SFR and old, are extrinsic (that is, the quantity scales with the amount of material in an object). This is a consequence of our model galaxy having a fixed size, expressed in terms of the fixed reference scale length of the old stellar population in B-band of  kpc.

The dust emission SED of the diffuse component of a galaxy with a value¹⁵ for differing from will be:(D.2)\arraycolsep1.75ptwhere (D.3)where SFR^model is the star-formation rate of the model galaxy having the reference size , old^model is the normalised luminosity of the old stellar disk population of the model galaxy having the reference size , SFR is the real star formation rate of the galaxy that we want to model, old is the real normalised luminosity of the old stellar disk population of the galaxy that we want to model, and (D.4)Equation (D.2) expresses the fact that radiation field energy density, and hence the colours of the dust emission, varies according to surface density of luminosity. In cases where a galaxy is unresolved ζ is unconstrained by the data, and becomes a further free parameter of the model. Since we may also not know the distance D to the galaxy, it is convenient to express the dust emission SED of the diffuse component from Eq. (D.2) as a flux density by dividing throughout by 4πD², to obtain: (D.5)where θ_gal is the half angle subtended at the Earth by the actual B-band scalelength of the galaxy:

(D.6)We note that, provided the galaxy is sufficiently resolved for θ_gal to be measured, as expressed by Eq. (D.5) depends on the value of ζ (via SFR and old – cf. Eq. (D.3)), but is independent of the distance D.

Our model is constructed such that the total emitted luminosity L_{λ,   star} of UV and optical light powering the dust emission can be directly constrained from available measured apparent UV/optical spatially integrated fluxes, , corrected for attenuation (expressed here in linear form by Eq. (D.1)) as a function of , SFR, F, old, B/D, i and distance D:(D.7) where, as outlined in Sect. 5 of Paper III and in Appendix C of this paper the dependence on SFR, old and B/D is for the optical range only. Depending on the exact range of wavelengths for which is available, L_{λ,   star} can be integrated over wavelength to obtain constraints on the parameters SFR = SFR^con and/or old = old^con as a function of and i. Specifically, if UV data is available we can require that (D.8)Combining Eq. (17) from Sect. 2.3.3 and Eqs. (D.7) and (D.8), this leads to: (D.9)Correspondingly, if optical data is available, we can require that(D.10) After manipulation of Eqs. (12), (13), (15) and (19) from Sects. 2.3.2, 2.3.2 and 2.3.3 and Eq. (D.7), this leads to: (D.11)where (D.12)(D.13)(D.14)Note that old^con is a weak function of SFR as well as , B/D and i due to the need to take into account the contribution of optical photons by the young stellar population, as described in Sect. 2.3.1.

We are now in a position to determine the physical model parameters from a combined set of measured UV/optical flux densities , measured at wavelengths λ^{iobs,   star}, and IR/submm flux densities of the pure dust component (corrected for contamination by direct stellar light at short infrared wavelengths) , measured at wavelengths λ^{iobs,   dust}. To do so, we minimise the function(D.15) subject, if UV data is available, to the constraint from Eq. (D.9) and, if optical data is available, from Eq. (D.11). σ_{iobs,   dust} are the 1σ uncertainties in the measurements and S_{λ,   dust} is given by (D.16)where and are given by Eqs. (D.5) and (41), respectively. In the case that the distance, the optical structure and orientation parameters D, θ_gal, ζ, B/D and i are known the optimization problem posed by Eq. (D.15) is reduced from the parameter set (SFR, old, , F, B/D, ζ) to the parameter set (SFR, old, , F). With just 4 parameters, this might potentially be solvable purely considering the dust emission data, bearing in mind the orthogonal effect of these parameters on the amplitude and colour of the dust emission SEDs, (as described in Sect. 6) if at least 4 data points are available well sampling the whole MIR/FIR/submm range. However, this will often be the exception rather than the rule, and in any case the fit is primarily and more robustly constrained through the optical and UV measurements. This is firstly because, if both UV and optical measurements are available the number of the primary search parameters would be reduced from four to just two – and F. Secondly, and perhaps more significantly, the model would no longer have any degree of freedom in terms of scaling parameters, due to the fact that, as noted above, SFR and old are the only extrinsic parameters in the full parameter set.

Below we illustrate how to use the optical and UV constraints by giving a possible processing path for a galaxy with a known (spectroscopically determined) redshift, for which integrated UV and optical photometry were available, and for which B/D, i and θ_gal are known from optical imaging. The parameter set to be determined is thus , SFR, old, F:

• step i): choose a trial value for each of and F

• step ii): set SFR to SFR^con and old to old^con from Eqs. (D.9) and (D.11) for the trial values of and F.

• step iii): find SFR^model and old^model from Eq. (D.3), substituting for ζ as defined in Eq. (D.4). The value of used in Eq. (D.4) should be derived from the optically measured value using the correction factors tabulated as a function of i and in Tables 1−5 of Paper IV.

• step iv): using the trial values of and F, together with SFR^model and old^model and ζ from step iii) find .

• step v): substitute from step iv) into Eq. (D.5) to compute . Compute (note the use of the extrinsic value of SFR = SFR^con here). Substitute the values for and determined at the wavelengths of the IR/submm observations in Eq. (D.15) to determine χ² for the trial combination of and F.

• step vi): repeat steps i) to v), until the combination of and F that minimises χ2 is found. The values of SFR and old from steps ii) to vi) found for the pair of and F that minimises χ2 are then best fit parameters in SFR and old.

We note that in the case of edge-on galaxies UV data should not be used to constrain SFR, since typically only a few percent of the total UV disk luminosity will be seen, and the solution can be subjected to stochastic variations, because the received photons may only be emitted by a small numbers of star-forming regions. In this case the SFR is better constrained from the FIR emission, as modelled for NGC 891 (Paper I and Sect. 3).

Having found the best fit parameter set, the following two further steps can be made:

• step vii: deredden the UV/optical/NIR spectrum using the fittedvalues of , F, and the measured B/D and i.

• step viii: apply a population synthesis modelling fit to the dereddened UV/optical spectrum from vii to find the star formation history.

Appendix E: Tables describing the fixed parameters of our model

---

Acknowledgments

We would like to thank an anonymous referee for constructive and useful comments that helped improve the quality of the paper. We would also like to acknowledge Aigen Li for providing us with the data for the optical properties of carbonaceous grains used in Draine & Li (2007).

One of us (NDK) acknowledges support from EU Marie Courie grant 39965 and EU REGPOT grant 206469.

References

All Tables

All Figures