Modelling the UV to radio SEDs of nearby star-forming galaxies: new Parsec SSP for Grasil


By means of the updated PARSEC database of evolutionary tracks of massive stars, we compute the integrated stellar light, the ionizing photon budget and the supernova rates of young simple stellar populations (SSPs), for different metallicities and IMF upper mass limits. Using CLOUDY we compute and include in the SSP spectra the nebular emission contribution. We also revisit the thermal and non-thermal radio emission contribution from young stars. Using GRASIL we can thus predict the panchromatic spectrum and the main recombination lines of any type of star-forming galaxy, including the effects of dust absorption and re-emission. We check the new models against the spectral energy distributions (SEDs) of selected well-observed nearby galaxies. From the best-fit models we obtain a consistent set of star formation rate (SFR) calibrations at wavelengths ranging from ultraviolet (UV) to radio. We also provide analytical calibrations that take into account the dependence on metallcity and IMF upper mass limit of the SSPs. We show that the latter limit can be well constrained by combining information from the observed far infrared, 24 , 33 GHz and H luminosities. Another interesting property derived from the fits is that, while in a normal galaxy the attenuation in the lines is significantly higher than that in the nearby continuum, in individual star bursting regions they are similar, supporting the notion that this effect is due to an age selective extinction. Since in these conditions the Balmer decrement method may not be accurate, we provide relations to estimate the attenuation from the observed 24  or 33 GHz fluxes. These relations can be useful for the analysis of young high redshift galaxies.

radio continuum: galaxies–infrared: galaxies–ISM: dust, extinction–stars:formation–galaxies: high-redshift

1 Introduction

Modelling the SEDs of galaxies has proven to be a very powerful tool in our current understanding of the different physical processes that come into play in the formation and evolution of galaxies. Various physical properties of galaxies like stellar, metal and dust content, star formation rate, dust obscuration, etc are estimated by fitting the theoretical SEDs to the observed ones. In the UV to infrared (IR) spectral regions of a SED, dust plays an important role. It absorbs (as well as scatters) the UV-near-infrared (NIR) light and re-radiates it in the IR. This IR emission (3 - 1000 ) may arise from (a) the emission from dust heated by young OB stars (neglecting heating from AGN), (b) the emission from the circumstellar envelopes of evolved stars and (c) the cirrus emission from dust distributed throughout in the ISM and heated by the general interstellar radiation field. Point (b) leads to wrong estimates of the physical properties related to star formation. The radio band which is not sensitive to dust attenuation is usually used to check and complement the interpretations arrived by using the optical/IR bands. Radio emission from normal star-forming galaxies is usually dominated by the non-thermal component (up to 90 per cent of the radio flux) which is believed to be due to the synchrotron emission from relativistic electrons accelerated into the shocked ISM by core-collapse supernovae (CCSN) explosions from massive young stars (above  ) ending their lifes (Condon, 1992). Recent advances in hydrodynamical simulations of CCSN have indicated a range of stellar masses where the stars fail to explode but rather end up directly as a black hole. This raises concern about the widely accepted notion that all stars more massive than end up as a CCSN.

Over the years, great progress has been made both in the development of tools and models use to extract the information encoded in the SEDs and in multi-wavelength surveys that sample the UV to radio SED of local and high redshift galaxies. Furthermore, new observing facilities with unprecedented resolution and sensivity (e.g. SKA, JWST, EVLT etc) will be put in place in the nearest feature. At the same time, semi-analytic models of galaxy formation, in the context of the current cosmological standard model, are appearing and providing realistic predictions of the physical properties of galaxies with which the results of the SED fitting can be compared.

Stellar population synthesis (SPS) still remains the basis of SED modelling. The most common method used in computing the SEDs of SSPs is that of isochrone synthesis (Chiosi, Bertelli & Bressan, 1988; Maeder & Meynet, 1988; Charlot & Bruzual, 1991) which uses the locus of stars in an isochrone to integrate the spectra of all stars along an isochrone to get the total flux. This involves computing first the stellar evolutionary tracks for different masses and metal contents and building the isochrones from the tracks. With the isochrones, the resulting stellar spectrum is computed using stellar atmosphere libraries. In recent years, significant efforts from different research groups have been put into providing homogeneous sets of evolutionary tracks and improving the stellar libraries. As a result of these improvements, SPS models can remarkably reproduce well the UV-NIR SEDs and high-resolution spectra in the optical wavelength band. However, despite these improvements, there are still challenges in the field, especially in the treatment of the phases of stellar evolution that are weakly understood. The most important being the short lived phases: massive stars, thermally pulsing asymptotic giant branch (TP-AGB) stars, blue stragglers and extreme horizontal branch stars. The TP-AGB stars have however increasingly received attention leading to a rapid progress in their modelling (Maraston, 2005; Marigo & Girardi, 2007; Marigo et al., 2017).

The stellar evolution code used in Padova to compute sets of stellar evolutionary tracks that are of wide useage in the astronomical community has recently been thoroughly revised and updated, going by the name, PARSEC(PAdova-TRieste Stellar Evolution Code). More details of this code can be found in Bressan et al. (2012, 2013); Chen et al. (2014) and will be briefly discussed later in the text. In this paper, our main aim is to use the PARSEC evolutionary tracks of massive stars to compute the quantities (the integrated light, ionizing photon budget and the supernova rates) of young SSPs, thereby updating the SSPs used by GRASIL in predicting the panchromatic spectrum of star-forming galaxies. We also revisit the thermal and non-thermal radio emissions. We finally check these new models with selected nearby well-observed galaxies. We perform an analysis of the resulting best-fit SEDs, with the aim of obtaining a new set of SFR calibrations and investigating the dependence of the dust attenuation properties on galaxy types.

The structure of the rest of paper is as follows: In Section 2, we use the new PARSEC code database of stellar evolutionary tracks of massive stars to compute the ionizing photon budget, the integrated light and the supernova rates predicted by young SSP models. Using the integrated spectra of the SPPs in CLOUDY, we computed the nebular emissions. In Section 3, we revise (a) the prescription used in computing the thermal radio emission. (b) the previous non-thermal radio emission model originally described in Bressan, Silva & Granato (2002, hereafter B02) while taking into account recent advances in CCSN explosion models. In Section 4, we check our new radio emission models and SSPs with GRASIL using selected well-observed nearby galaxies. Finally, we discuss the resulting best fit SEDs of all galaxies studied in this paper in Section 5. We draw our conclusions in Section 6.

Throughout this paper we adopt (Asplund et al., 2009) as the solar oxygen abundance. Our models adopt a Kennicutt (1983) IMF (slope of -1.4 for and -2.5 for ). The cosmological parameters we adopt assume

2 Simple Stellar Populations with Parsec

PARSEC is the latest version of the Padova-Trieste stellar evolution code with thorough update of the major input physics including new and accurate homogenised opacity and equation of state tables fully consistent with any adopted chemical composition. More details may be found in Bressan et al. (2012, 2013); Chen et al. (2014). The evolutionary tracks span a wide range in metallicities, , and initial masses, from very low to very massive stars, starting from the pre-main sequence phase and ending at central carbon ignition. Tang et al. (2014) and Chen et al. (2015) computed new evolutionary tracks of massive stars and tables of theoretical bolometric corrections that allow for the conversion from theoretical HR to the observed colour-magnitude diagrams. A preliminary comparison of the new models with color-magnitude diagrams of star-forming regions in nearby low metallicity dwarf irregular galaxies was performed by Tang et al. (2014). The full set of new evolutionary tracks and the corresponding isochrones may be found in and, respectively.

2.1 Integrated SSP spectra

To build the integrated spectra of SSPs, we adopt the spectral library compilation by Chen et al. (2015) who homogenized various sets of existing stellar atmosphere libraries encompassing a wide range of parameters of both cool and hot stars (i.e. masses, evolutionary stages and metallicities). The atmosphere models adopted by Chen et al. (2015) are the ATLAS9 models (Kurucz, 1993; Castelli & Kurucz, 2004), suitable for intermediate and low mass stars, and the Phoenix models (Allard et al., 1997) for the coolest stars. For hot massive stars Chen et al. (2015) generated new wind models with the WM-BASIC code Pauldrach, Puls & Kudritzki (1986) and adopted new models of WR stars from the Potsdam group (see e.g. Sander et al., 2015). Since different sets of libraries have different metallicities, they were merged to obtain homogeneous sets of spectral libraries as described in Girardi et al. (2002) and Chen et al. (2014).

Figure 1: SSP Integrated spectra per solar mass at ages of 2, 5 and 40 Myr, for metallicities Z = 0.02, 0.004 and 0.0001 and of 40, 120 and 350  indicated by the solid blue, dotted red and dashed green lines respectively. At 2 Myr, the spectra for the various upper mass limits at any of the three metallicities is quite distinguishable. At 5 Myr and Z = 0.02, it is distinguishable only at wavelengths below , the maximum for the lyman ionizing photons. At this age, notice that the spectra for the of 120 and 350 below are superimposed, owing to the fact they must have evolved off the main sequence at this age. At 5 Myr and Z = 0.0001, the reverse is the case, the spectra is distinguishable only at wavelengths above because at this metallicity, stars lifetimes are a bit longer than at Z = 0.02. At 40 Myr, the spectra for all upper mass limits and metalicities are all superimposed and the Lyman ionizing photons are no longer produced.

The properties of SSP are obtained by integration along the corresponding isochrones assuming a two-slopes power law Kennicutt initial mass function. Results are presented for three different values of the upper mass limit, of 40,120, 350 . As an example, Figure 1 shows the synthesized spectra of star clusters at ages of 2, 5 and 40 Myr and for three values of the initial metallicity, Z = 0.02, 0.004 and 0.0001. The effect of assuming different values, at fixed total initial mass, is more pronounced at an age of 2 Myr. At 5 Myr, when stars with masses of 120  and 350  already left the main sequence or already died out, the spectra of SSP with upper mass limits of 120  and 350  are superimposed, while those of 40  remains distinct. At 40 Myr the spectra are almost independent of , though it is evident that the spectrum of the SSP with larger has a lower luminosity, because of the higher mass budget stored in massive stars We already note from this figure that the effects of on the number of ionizing photons disappears at ages larger than about 5 Myr, as discussed below.

2.2 Ionizing Photon Budget

The number of Lyman ionizing photons per sec () and per unit mass emitted by young stellar populations is controlled by hot massive stars, i.e. O-B main sequence stars and Wolf Rayet stars. This number is thus critically dependent on the shape of the initial mass function in the domain of massive stars. In Figure 2, we show the time evolution of for selected SSPs with different upper mass limit () and metallicities. As already said, the adopted IMF is a two-slope power law Kennicutt (1983). The lower limit is  = 0.1 and the upper limits are  = 40 ,  = 120 and  = 350 . The slope of the IMF is X = 1.4 between and  = 1 and X = 2.5 between  = 1 and

Figure 2: Variation of the instantenous number of ionizing photons per second per unit ionizing solar mass of the cluster with age (Myr). Different symbols correspond to different metallicities as illustrated in the plot


From this figure we may appreciate the role of age, metallicity and IMF on the ionizing photon rate . For a given Z and , generally increases to a maximum value and then, once a threshold age is reached, it rapidly decreases to negligible values. At fixed , both the maximum value of and the threshold age decrease, at increasing metallicity. In general this is also true at varying age, i.e. at given , decreases at increasing metallicity. However there are some cases where this is not true, in particular for the SSP of solar metallicity. The effect of is strong. A star cluster with of 350  produces about seven times more ionizing photons than a cluster with a of 40 , at constant total mass. Moreover, the age to attain the maximum becomes lesser at increasing upper mass limit, reflecting the larger relative weight of more massive stars in the ionizing photon budget.

2.3 Nebular Emission with PARSEC’s Ssp

The integrated spectra of the SSPs are used to calculate the nebular emission from the surrounding \ionHII regions which is then added to the original spectrum to obtain the integrated spectra containing both the stellar continuum (with absorption lines and eventually wind emission features) and the nebular features (continuum and lines). For this purpose, star clusters are assumed to be the central ionizing sources of the \ionHII regions that are modelled using the photoionization code CLOUDY (Ferland, 1996). As a further input to the CLOUDY code, we specify that the \ionHII region is modelled as a thin shell of gas with constant density, , placed at a distance  = 15 pc from the central source. The evolution of the ionizing star clusters is followed from 0.1 Myr to 20 Myr, for five different values of the metallicity (Z = 0.0001, 0.0005, 0.004, 0.008, 0.008 and 0.02) and three values of of 40, 120 and 350 M. We note that our goal is not that of providing a detailed dependence of a large ensemble of emission lines on the critical parameters of the \ionHII nebulae. Instead we aim at obtaining a reasonable estimate of the intensities of the main recombination lines and of free-free emission to increase the diagnostic capabilities of our population synthesis codes. Line emission and the corresponding nebular continuum are much less dependent on the characteristic properties of the \ionHII regions (e.g. ionization parameter, individual abundance of heavy elements etc.) than e.g. excitation lines, for which a more detailed set of initial parameters would be more appropriate (see e.g. Panuzzo et al., 2003; Wofford et al., 2016). The main output of this process is a library of emission line intensities, nebular continuum properties and electron temperatures (T) of the \ionHII regions. Then, emission lines and nebular continuum are used to suitably complement the integrated SED of SSPs from the far-UV to radio wavelengths.

3 Radio Continuum Emission

Radio emission associated with the presence of young massive stars comprises essentially two components, the thermal radio emission (also referred to as free-free emission) and the non-thermal radio emission (also referred to as non-thermal radio emission). The former emission comes from the nebular free electrons originating from the ionizing radiation of massive stars. The non-thermal emission is instead believed to be synchrotron radiation that originated from the interaction of relativistic electrons, produced in the ejecta of core-collapsed supernovae (CSSN), with the ambient magnetic field. The radio continuum is thus a tracer of the number of massive stars formed (and exploded) and hence an optimal indicator of the very recent (if not the current) star formation rate.

The fact that the radio emission is a good SFR tracer is supported by the remarkably tight correlation between FIR and non-thermal radio emission. At 1.49 GHz, this correlation is quantified by the q-parameter (Helou, Soifer & Rowan-Robinson, 1985):


where (Young et al., 1989), and are IRAS flux densities in .

3.1 Thermal Radio Emission

At sufficiently high frequency such that free-free self-absorption is negligible, the relation between the specific luminosity of free-free emission and the rate of ionizing photons can be written as (Rubin, 1968; Condon, 1992)


where is the electron temperature, and is the velocity averaged gaunt factor (Draine, 2011):


where is the charge of the ions in the \ionHII region. An approximate velocity averaged gaunt factor were obtained earlier by Oster (1961):


Taking into account that and adopting a common approximation used in literature, we obtain for the free-free emission (Condon, 1992)


The latter equation is often used by assuming a fixed electronic temperature (10 K) to provide useful analytical approximations at radio frequencies especially in the lack of nebular emission (see discussion in B02). Here we directly estimate the intensity of the thermal radio emission of our SSP from Cloudy with the procedure already described in Section 2.3. In Figure 3, we show the evolution of the 1.49 GHz and 33 GHz thermal radio emission for different metallicities, Z = 0.0001, 0.0005, 0.004, 0.008 and 0.02, and different of and . The effects of metallicity and on the thermal emission are easily noted. For example, at Z = 0.0001, the 1.4 GHz thermal emission for is about 3 - 7 times larger than that for In general by increasing the metallicity from Z = 0.0001 to Z = 0.02 the thermal radio emission decreases by a factor of about 3. Instead by increasing from 40 to 350 , thermal emission increases by about one order of magnitude. These important factors, that must be considered in the calibration of the relation between star formation rate and thermal radio emission, are discussed below.

Figure 3: Variation of the 1.49GHz (red) and 33GHz (blue) thermal radio emission with age at different metallicities for different values of . The five different line styles correspond to the five different values of metallicity.

Metallicity effects

Equations 2, 3 and 5 contain an explicit dependence on the electron temperature, that is generally neglected in analytical approximations which assume a constant value of  = 10 K. Since is known to depend on the metallicity of the \ionHII regions, with a variation of more than a factor of two in the range of the observed values, we provide in Appendix A some useful analytical relations. We first list in Table 2 the values of the electronic temperature derived using CLOUDY for our SSP models at various metallicities and upper mass limits. In Figure 17 of Appendix A, we show the relations between oxygen abundance and the electronic temperature. In this figure, crosses, ’X’s and asterisks indicate , and respectively. The average of the empirical fits derived by López-Sánchez et al. (2012) for high-ionization \ionOIII and for low-ionization \ionOIII zones is shown as the dashed black lines. We provide a multiple regression fitting relation (Equation 34) between , and Z that could be easily included in analytical approximations.


In a young stellar system the ionizing photon budget is dominated by massive stars and thus there must be a tight relation between the current star formation rate and the rate of ionizing photons, that ultimately produce the thermal radio emission Condon & Yin (1990). In this section, we provide a calibration of the relation. For this purpose we consider the integrated spectrum of a galaxy up to a time :


where SFR(t’) is the star formation rate at time and is the stellar spectrum of a SSP of age and given metallicity


In the latter equation, is the spectrum of an individual star in a SSP which depends on the fundamental parameters, its mass , age , metal abundance and the IMF, . Using the SSPs already described, a constant SFR of 10 and adopting the Kennicutt (1983) IMF, we obtain the integrated number of ionizing photons , shown in figure 4. Since is dominated by the most massive stars, the number of ionizing photons emitted by a galaxy with constant SFR will initially grow and soon saturate to a constant maximum value, when there is a balance between the newly formed ionizing massive stars and the ones that die. Looking at Figure 2 we see that, almost independently from the metallicity, the characteristic time of the saturation is set by the rapid drop of above about 6 Myr. The effect of on is a direct consequence of the variation of the of SSPs shown in Figure 2. The variation of with metallicity decreases with increasing upper mass limit.

Figure 4: Evolution of the integrated number of ionizing photons per second per unit solar mass initially formed as a function of the age in yr. Different line colors represents different metallciities. The variation of Q(H) with metallicity decreasis at increasing IMF upper mass limit.

After denoting by the calibration coefficient between SFR and Q(H) in the equation below


we collect in Table 2 the values of obtained with our constant SFR models, for different SSPs parameters. To illustrate the significant variations of with metallicity at a given , we show in Figure 18 the plot of against metallicity for  = 40, 120 and 350 Using the values of given in Table 2 for different and , we provide a multiple regression fitting relation between , and metallicity (), Equation 35. As an example, using this fitting relation, for and Z = 0.02 Equation 8 is:


We may compare equation 9 with the one obtained by Murphy et al. (2012) using the Starburst99 stellar population model with a Kroupa (2001) IMF, a metallicity of Z = 0.02 and a constant star formation over 100 Myr.


We see that the Murphy et al. (2012) calibration constant is fairly in agreement with ours.

Sfr - Calibration

By combining equation 8 and equation 2, we derive the relation between the SFR and thermal radio emission:


where . The values of the coefficient for different SSP parameters are provided in Table 2 and the variation of with metallicity and is shown in Figure 19. Using these values of given in Table 2 for different and , we also, as we did for the case , provide a multiple regression fitting relation between , and metallicity (), Equation 36. As an example, using this fitting relation for and Z = 0.02 and assuming the already quoted approximation for the Gaunt factor, Equation 11 becomes


A similar equation is provided by Murphy et al. (2012):


SFR versus H and H calibrations

Using well-known relations between and the intensity of recombination lines (Osterbrock, 1989), we may also obtain the corresponding calibrations for the SFR. For H and H we have


Using the above relations in equation 8, we write in Appendix A analytical equations for the SFR-H (equation 39) and SFR-H(equation 40) calibrations, as a function of Z and M. For the case with and Z = 0.02 we obtain with these analytical equations


Our calibration coefficient in equation  16 is in good agreement with the value of obtained by Bicker & Fritze-v. Alvensleben (2005) using GALEV synthesis code.

3.2 Non-Thermal Radio Emission

We know quite little about the source of non-thermal (NT) radio emission in star-forming galaxies. The mechanism is believed to be synchrotron emission from relativistic electrons that are accelerated by ISM shocks in the outskirts of CCSN explosions (Condon, 1992). B02 derived the following relation between NT radio emission and the CCSN rate

In the above equation is the CCSN rate, is the average non thermal radio luminosity due to a young Supernova Remnant (SNR) and is the average injected energy in relativistic electrons per CCSN event. Moreover, B02 estimated that, at 1.49 GHz,


can be estimated from observations. For our Galaxy B02 adopt (Berkhuijsen, 1984) at 0.4 and . They obtain for a radio slope , a value of = 1.44. We note that the final fate of massive stars is a critical assumption for estimating their contribution to thermal radio emission. The formation of NSs and BHs depends on the amount of mass lost by the massive star through stellar winds and on the hydrodynamics of the explosion. One of the the most troublesome facts is that the variations shown by these pre-supernova stars in their structural properties, like for e.g. the Fe-core and O-core masses, are non-monotonic and pronounced even within small differences in the ZAMS masses (Sukhbold & Woosley, 2014). Some recent works in the attempt to characterize the parameters of successful and failed supernovae have yielded some structural parameters that can be utilized in predicting the fate of SNe (eg. Ertl et al., 2015; Ugliano et al., 2012; Janka, 2012; O’Connor & Ott, 2011). Using these parameters, Spera, Mapelli & Bressan (2015) and Slemer et al. 2017 (to be submitted) were able to characterize the final fate of PARSEC massive stars for the different criteria adopted for successful CCSN explosion. Following Spera, Mapelli & Bressan (2015) and Slemer et al. 2017 we assume that stars with initial mass above about do not contribute to non thermal radio emission, contrary to B02 where all stars of masses above were thought to explode as CCSN. The region between about is critical because the explosion criterion by O’Connor & Ott (2011) produces a much higher number of NSs than that produced by Ertl et al. (2015) criterion. Finally we stress that we assume that progenitors undergoing pair-instability SNe either collapse to a BH or are completely incinerated by a thermonuclear explosion without producing relativistic electrons, similarly to what is assumed for SNIa. This assumption on the CCSN distribution with initial mass modifies the expected non thermal luminosity of star-forming galaxies and requires a re-calibration of the GRASIL model. After accounting for the failed and successful SNe in the evaluation of non thermal emission from Equation LABEL:eq_l_nt_bre, we show in Figure 5 the variation of the 1.49 GHz non-thermal radio emissions with age, for different upper mass limits and metallicities. We note that non-thermal radio emission is not sensitive to the upper mass limit of the IMF, as long as it is larger than the assumed threshold for failed SN, M = 30. In this case non thermal emission begins at an age of 7 Myr and not earlier (3 Myr) as is the case in B02. By coincidence, this critical time corresponds to about the epoch at which thermal emission abruptly fades down, as depicted in Figure 5. Non thermal emission remains effective up to about 35 - 38 Myr, depending on the metallicity. In subsequent sections we will show how thermal and non-thermal radio emission are checked and eventually re-calibrated using the GRASIL code.

Figure 5: Variation of the 1.49 GHz and 33 GHz non-thermal (free-free) emission with age at different metallicities and different IMF upper mass limits. The red lines indicate thermal radio emission and the blue lines non-thermal radio emission. The plot for the thermal radio emission has already been shown in Figure 3 but added here for comparison with the non-thermal radio emission. As easily noticed in the figure, the non-thermal emission is almost unchanged for each of the different IMF upper mass limits, owing to the fact that masses above 30  collapse directly to BH and do not explode to produce non-thermal emission. For this same reason, the non thermal emission for all plots starts at about an age of 7 Myr. For comparison purposes, the plot of the thermal radio emission is included to show it’s pronounced variation with upper mass limits. At 1.4 GHz (upper panel), the non-thermal radio emission dominates that at 33 GHz (lower panel) whereas at 33 GHz, the thermal emission dominates.

4 Calibration of the new SSP suite with Grasil

In the previous sections, we described a new suite of SSPs that will supersede the ones used in the current version of GRASIL. The new suite differs in many aspects from the previous one as we briefly list in the followings. (1) The SSP are based on the most recent PARSEC stellar evolutionary tracks with an updated physics, in particular new mass-loss recipes and finer and wider coverage in initial mass and metallicity; (2) they include the most recent advances in our understanding of the CCSN explosion mechanism and account for the so called ”failed SN″; (3) a more accurate gaunt factor and metallicity-dependent electron temperature were incorporated in the relation between the integrated ionization photon luminosity and thermal radio luminosity (4) in the current suite, we may adopt several values of the IMF upper mass limit Because of all these differences, we need to check and recalibrate some parameters of the SED produced by GRASIL, in particular to check whether we are able to reproduce the canonical value of , observed in a prototype normal star-forming galaxiy, see Eq. 1. For this purpose, in the next section we will use GRASIL with the new SSPs to reproduce the SEDs of a some selected well studied galaxies.

4.1 SED Modelling with Grasil

GRASIL is a spectro-phometric code able to predict the SED of galaxies from the FUV to the radio, including state-of-the-art treatment of dust reprocessing (Silva et al., 1998, 2011; Granato et al., 2000), production of radio photons by thermal and non-thermal processes (BO2) and an updated treatment of PAH emission (Vega et al. (2005)). The reader is referred to these original papers for a fully detailed description of the code. It is also worth noting that nebular emission was already included in GRASIL by Panuzzo et al. (2003) but with a completely different procedure which also accounted for different electron densities. In this respect our procedure is more simple because nebular emission is added directly to the SSPs, but it allows a more versatile use of GRASIL.

For sake of convenience, we briefly summarise here the main features of GRASIL. The first step is to compute a chemical evolution model that provides the star formation history (SFH) and the metallicity and mass of the gas. Other quantities are computed by the code CHEEVO, such as mass in stars, SN rates, detailed elemental abundances, but they are not used for the spectro-photometric synthesis. GRASIL’s main task is to compute the SED resulting from the interaction between the stellar radiation from CHEEVO and dust, using a relatively flexible and realistic geometry where stars and dust are distributed in a spheroidal and/or a disk profiles. A spherical symmetric distribution with a King profile is adopted in the case of spheroidal systems while, for disk-like systems, a double exponential of the distance from the polar axis and from the equatorial plane is adopted. The dusty environments are made up of dust (i) in interstellar HI clouds heated by the general interstellar radiation of the galaxy (refered to as the cirrus component), (ii) associated with star-forming molecular clouds (MC) and (iii) in the circumstellar shells of Asymptotic Giant Branch (AGB) stars. GRASIL performs the radiative transfer of starlight through these environments. We remind the reader that the reprocessing of starlight by dust in envelopes of AGB stars is already taken into account in our SSPs. For the intrinsic dust properties, GRASIL adopts grains (graphite and silicate) in thermal equilibrium and a mixture of smaller grains and PAHs fluctuating in temperature. The dust parameters are tuned in order to match with the observed emissivity and extinction properties of the local ISM (Vega et al., 2005).

M100: a close to solar metallicity galaxy

The galaxy M100 was selected from the Spitzer Nearby Galaxies Survey (SINGS) (Kennicutt et al., 2003), a collection of 75 nearby galaxies with a rich data coverage from the far-UV to the radio, because it is one of the best sampled objects, including the intensities of main emission lines. This will provide the opportunity to check, for the first time in GRASIL, the consistency of the SED continuum and emission lines fitting. To obtain the chemical evolution model we adopt the parameters of the chemical evolution code CHE-EVO given by Silva et al. (1998) who was able to well-reproduce M100. They are infall time scale  Gyr, efficiency of the Schmidt-law star formation law and the total infall mass = . Figure 7 shows the SFH resulting from the adopted chemical evolution parameters. The SFH is indicated by the solid blue line and the gas mass history by the dotted purple line. To perform the SED fit, we consider several time steps along the chemical evolution of the galaxy for which we built a library of SED models with different GRASIL parameters. We note that since the observed SED refers to the whole galaxy, the derived parameters correspond to a luminosity average of the physical properties, as commonly obtained by this kind of fitting process. As a particular check, we also investigate the effect of changing the upper mass limit of the IMF, ,120 and , on some physical quantities derived from our best-fit SED, in particular the SFR. The main GRASIL parameters varied are the molecular gas fraction , the escape time and the optical depth of molecular clouds (MC) at 1 . The best fits were obtained by minimizing the merit function which is given by


where , and are the model flux values, the observed fluxes and observational errors respectively. is the number of photometric bands used to obtain the best fit. Our GRASIL best-fit SEDs of M100 for and , are shown in the left and right panels of Figure 6 respectively, while the corresponding best fit parameters are summarized in Table  1. We do not show the case with , because its thermal-radio emission was significantly exceeding the observed one. This discrepancy could not be cured by varying any other parameter in the fit. This result is, by itself, quite interesting and shows that for a normal star-forming galaxy it is difficult to have an average IMF extending up to such high initial masses. The different dust components in diffuse medium (cirrus) and molecular clouds are shown and indicated by the blue and red dashed lines, respectively. The thermal and non-thermal radio emission components are indicated by the cyan and purple dashed dotted lines, respectively. The thermal component can be seen to be negligible for the case of . We include, for the first time in GRASIL, selected emission lines in the SED. The labels in the plots refer to some important quantities derived from our best fit: the bolometric luminosity from UV to radio (BOL), the total IR luminosity from 3 - 1000 (FIR), the predicted extinction uncorrected luminosity (Ha), the observed luminosity (Ha), attenuation (A(H)), attenuation in the V() band (A(V)), the SFR averaged over the last 100 Myr (), the q-parameter as defined by equation 1 (q(1.4GHz)), the coefficients of the SFR(H) and SFR(IR) calibrations ( and respectively). A summary of other important quantities derived from our best fit is provided in Tables 3 and  4. For both values of the adopted , the best fits match very well the overall observed UV-radio SEDs. However, for = 120, the model over-predicts by factor of 3 the observed (i.e. non corrected by extinction, hereafter ’transmitted’) H luminosity, (taken from Kennicutt et al. (2009)). This upper mass limit also produces a higher UV emission and a slight excess of thermal emission. On the other hand, the observed H luminosity is well reproduced in the case of . The over-predicted H luminosity resulting from adopting = 120  could be lowered by increasing the escape time of young stars in their parent clouds, leading to a larger absorption in the MCs. This will also lower the predicted far-UV emission which, in the current best fit, is slightly larger than the observed value. However, a larger absorption from MCs will increase the predicted flux above the observed one. Note that at the cirrus component has a pronounced minimum and its contribution to the overall MIR emission is only a small fraction of the total. Thus, the flux, being dominated by the MC emission, is indeed a strong diagnostic for the amount of light reprocessed by the MC component. From our best fit, we obtain a CCSN rate to NT radio luminosity ratio (Equation LABEL:eq_l_nt_bre) that is about a factor of 1.35 larger than the value obtained by B02 using the previous radio model. That is,


We anticipate here that this value is confirmed by all subsequent SED fits. We also note that this value is almost independent from the adoption of for the reasons already discussed previously. The predicted average SFR resulting from the panchromatic fit is only slightly affected by the adoption of a different value of . By increasing from 40 to 120 the average SFR decreases by about 16 per cent.

Figure 6: M100 GRASIL best-fit SED for = 40  (right panel) and (left panel). The different dust components, diffuse medium and molecular clouds are represented by the dashed blue and red lines respectively. The thermal and non thermal radio components are represented by the cyan and purple dashed dotted lines. The thermal component can be seen to be negligible for the case of . For the different , note the differences in the predicted quantities labelled in the plots (see text for more details), in particular the H luminosities, attenuations (A(H)) and SFR(H) calibration coefficient. We estimated , a factor of 1.35 larger than that obtained in B02. Our average value of q = 2.4 which is in line with the observed value of 2.35 (Equation 1) is an indication of our good calibration.
Figure 7: Star formation history of M100 adopted in our CHE-EVO model. The model parametrs are labelled: galaxy’s age (solid red line) , infall time (dashed green line) and schmidt-law efficiency. These parameters were adopted by Silva et al. (1998) based on constraints to reproduce the observed final mass of M100. The SFH is indicated by the solid blue lines while the gas mass history by the dotted purple line.

NGC 6946: individual star-bursting regions

In the previous section we were able to reproduce fairly well the observed SEDs of the prototype galaxy M100, in particular, we used the best fit model to calibrate the constant of the non-thermal radio emission model. In this section we wish to check the thermal component of our radio emission model. For this purpose, we compare our synthetic SEDs with those of selected extranuclear regions of NGC 6946, a well studied starburst galaxy at a distance of 6.8 Mpc. This galaxy is dominated by very young starburst regions and shows relatively large thermal over non-thermal ratios (eg Israel & Kennicutt, 1980; Heckman et al., 1983), and thus is particularly well suited for a test of the free-free emission originating from star formation. Murphy et al. (2010) observed these regions at 1.4, 1.5, 1.7, 4.9, 8.5 and 33 GHz and complemented their full SEDs with existing data from the UV to the submm range. For eight observed regions Murphy et al. (2010) estimated an average 33 GHz thermal component of 85 per cent of the total, while for one region (named extra-nuclear region 4), this percentage was found to be significantly lower, 42  per cent, likely for the presence of a so-called anomalous dust emission component which suppresses the thermal component. Particulary interesting for our purpose is that Murphy et al. (2010) also provide the intensity of the H emission for each region, that can be directly compared to the predictions of our model.

To model the SEDs of the extra-nuclear star-forming regions of NGC 6946 we adopt a spherical symmetric distribution with a King profile for the stars and dust. Most important, we use the data not corrected for the local background emission of the galaxy, though (Murphy et al., 2011) provide also data corrected for this emission. We thus adopt, for each region, a chemical evolution model composed by a star-burst superimposed to a quiescent star formating component. Our choice is dictated by the fact that the subtraction by (Murphy et al., 2011) may bias our results and we prefer to eventually split the contribution of the starburst and the old disk components directly from our fits.

Figure 8: SFH of one of the 8 extranuclear regions, NGC 6946_enuc_3 for case. A starburst age of  = 14 Myr (green solid line) produced the best fit but we added = 5 Myr (blue dotted line) case only to illustrate that we do not expect any contribution to the non-thermal emission from the starburst at this age according to our radio model in Figure 5. At burst ages between 7 and 30 Myr, we expect a contribution to the non-thermal radio emission also from the starburst. The red dotted-dashed line indicates the age (11.000 Gyr) at which the SED of this region was observed.

For purposes of clarity we show in Figure 8 a plot of the star formation history of the extra-nuclear region 3 for . The corresponding chemical evolution model parameters are labelled in the left side of the figure. For these models, the galaxy is observed at a given age  Gyr. The ongoing starburst has an exponentially declining SFR that begins at , where t is the age of the starburst. Two SFR examples are given for a starburst age of t = 5 Myr and t = 14 Myr respectively. In the first case we are not expecting any contribution to the non-thermal emission from the starburst while, in the second case since the age is larger than 7 Myr (see Figure 5) we expect a contribution to the non-thermal radio emission also from the starburst. Thus any age between t =  7 and 30 Myr can provide SEDs with different percentages of thermal vs non-thermal radio emission. We note that this can be used to mimic possible different contributions of non thermal emission from the underlying disk. To obtain the GRASIL best fit models we mainly varied the following parameters: the escape time of young stars from their birth clouds , the optical depth of molecular clouds at 1 , molecular gas fraction , the submm dust emissivity index and the start time of the burst (represented by the starburst age ). The parameters of the best fit SEDs are shown in Table 1 while, the corresponding plots, are shown in Figure 9. In this figure, left panels show the results obtained adopting and right panels show those for .

Figure 9: GRASIL best fit SEDs for the extra-nuclear regions of NGC 6946 (red diamonds) for upper mass limit of 40 (left panels) and 120 (right panels). To test our radio decomposition, we over-plotted on top our thermal radio component (dotted-dashed cyan lines) the local background subtracted radio data flux (cyan diamonds) by Murphy et al. (2011).
Figure 10: continued

For all regions the observed IR SED can be fairly well reproduced by using either of the IMF upper mass limits. In general, the value of is found to lie between 2.5 and 2.6. This is about 0.2 dex larger than the value observed in the normal star-forming galaxy M100 (  =  2.35) implying that, in these star-forming regions, the ratio between radio and FIR luminosity is about a factor 1.6 lower than in the normal star-forming galaxy M100. The estimated young starburst ages support the notion that this is due to a lack of non-thermal radio emission as predicted by the original models by B02. This is also evident from the radio slope which is found to be flatter than the value observed in the normal star-forming galaxy (- 0.8). It is worth reminding here that we adopted the calibration () of non-thermal radio luminosity obtained from the fits to M100. We note also that since the radio emission is dominated by the thermal component, which has no free parameters, the SFR in these regions is very well determined by the 33 GHz point, modulo the IMF. Another characteristics of the fits is the difficulty of reproducing the UV data, in spite of a fairly well fit in all the other bands from mid-IR to radio bands. In particular the runs with tend to show larger UV fluxes than the corresponding cases made with , which is against what could be expected from the dependence of the UV luminosity on the IMF. This can be noted in Figure 9 where we see that adopting we either reproduce or underpredic the far-UV data while, adopting , we either reproduce or overpredict them. At the same time, the average SFR in the models with is about two to three times larger than that obtained with the models that adopt . This is in evident contrast with what we have found in the fits of normal galaxies where we see that the SFR, obtained with different values of , differ by 20 per cent at maximum and it is a consequence of the young age of such star bursting regions, where the bolometric luminosity is dominated by the most massive stars. Finally, the attenuation in the models with is always lower than that obtained with the models adopting . As the non-thermal radio emission starts about 7 Myr after the beginning of the burst of star formation (see Figure 5), we have the opportunity to perform an accurate decomposition of the radio flux into thermal and non thermal components. We present in Table 4 the thermal fraction resulting from this decomposition and we show the two radio components in Figure 9. To show the quality of our decomposition, we also show in Figure 9 and only for the case with , the background subtracted radio data fluxes (cyan diamonds) by Murphy et al. (2011). In the next section, we will discuss the SFR calibrations derived using our best fit model of the galaxies studied in this work.

(Gyr/Myr) (Myr)
(1) (2) (3) (4) (5) (6)
= 40 M☉
M100 12.0 2.5 2.0 12.0 0.10
NGC 6946_1 8.0 1.0 2.1 9.07 0.30
NGC 6946_2 9.0 0.5 1.7 12.00 0.05
NGC 6946_3 8.0 1.0 2.0 24.48 0.20
NGC 6946_5 8.0 0.3 2.0 4.15 0.25
NGC 6946_6 12.0 1.2 2.2 5.33 0.40
NGC 6946_7 7.0 1.0 2.1 14.81 0.30
NGC 6946_8 7.0 0.6 2.1 6.12 0.22
NGC 6946_9 8.0 1.0 2.1 5.33 0.30
= 120 M☉
M100 12.0 0.9 2.0 14.81 0.10
NGC 6946_1 11.0 0.6 2.2 5.33 0.12
NGC 6946_2 11.0 1.0 1.7 16.60 0.15
NGC 6946_3 8.0 0.9 2.2 12.00 0.14
NGC 6946_5 18.0 1.0 1.8 18.75 0.50
NGC 6946_6 11.0 0.7 2.2 5.33 0.12
NGC 6946_7 8.0 1.0 2.2 5.33 0.16
NGC 6946_8 8.0 0.9 2.2 5.33 0.12
NGC 6946_9 9.0 0.6 2.2 5.33 0.12

Column (2) gives the age of the galaxy () in Gyr for the normal star-forming galaxy M100 and the age of the burst () in Myr for the NGC 6946 star-bursting regions. Parameters in other columns are as described in text.

Table 1: GRASIL best fit parameters (for upper mass limits of 40 and 120) for the modelled SEDs of M100 and NGC 6956 SF regions.

5 Discussion

In this section we discuss the results obtained from the best fits of the SEDs of M100 and the extra-nuclear regions of NGC 6946. We begin with the SFR calibrations obtained with the new SSPs and then discuss the impact of the newly added emission line prediction on the resulting attenuation.

5.1 Star formation rate calibrations from Grasil best fits

We first show in Table 3 the luminosities of the best fit models in the selected photometric bands, from the UV to radio wavelengths. The upper panel refers to the results obtained adopting while, the lower panel refers to those obtained with . The FIR luminosity in column 11 was obtained by integrating the IR specific luminosity from 3 to 1000 . We also show the predicted intrinsic (suffix ) and transmitted (suffix ) intensities of the H and H emission lines, the value of (Equation 1), the relative contribution to the 3 - 1000  FIR luminosity by the molecular clouds component and by the cirrus component, and the ratio between the 3 - 1000  FIR and the bolometric luminosity. In the columns of radio luminosities, we enclose in parenthesis the fractional contribution of the thermal radio component to the total radio emission, derived from the models.

Effects of Mup
By construction, i.e. since we are discussing best fits SEDs, the predicted luminosities obtained with the two different are very similar, but for those in the emission lines and in the UV bands. These are the most difficult to model because they are the most sensitive to the attenuation. This is not the only reason however, because even the intrinsic H and H luminosity differs by a significant factor in the two cases. Indeed, in the two cases the ionizing photon flux (and the far and near-UV) differ much more than the corresponding bolometric luminosities. Said in another way, one may be able to obtain the same bolometric flux with two different IMF upper limits but the amount of ionizing photons may be significantly different, and this mainly affect the intensity of the recombination lines and the free-free radio emission. Note however that the 33 GHz luminosity does not show the same strong dependence on the adopted Mup as that shown by the recombination lines, and this is likely due to the fact that a non negligible contribution of the non-thermal luminosity is present even at this high radio frequency, which should be larger in the case of lower . The best fit transmitted intensities of the H emission lines are compared with the observed ones in Figure 11, for both cases of . We do not include M100 in this plot to better show the case of the star-forming regions. The solid line indicates the one-to-one correlation. We see that, in the case of extra-nuclear regions, the predicted values of the case with are larger than those of the case with .

Figure 11: Comparisons between our model’s transmitted intensities of the H emission lines with the observed ones, for (red daimonds) and (blue triangles) cases. The solid line indicates the one-to-one correlation whereas the horizontal lines indicate the errors in the observed H line. The green daimond indicates the case with for the extra-nuclear region 8.

In general our best fit models are able to reproduce the observed values with an accuracy of about 50 per cent. Looking to individual regions, we see that regions 1 ,2, 3, 5, 6 are compatible with while, for regions 9 and 7, a larger value of seems suggested, up to . From the figure it also appears that, in order to fit region 8, an even higher IMF upper mass limit (see the green daimond in Figure 11 for the case of ) should be used. Furhermore, we note that the in all cases where we are able to reproduce the observed H emission we are also able to reproduce fairly well the observed FUV luminosity, which is also strictly related to the intrinsic ionizing photon flux. On the contrary the NUV flux is not well reproduced by our models, likely indicating that its origin is not strictly related to the star formation process. Since all model reproduce the observed bolometric luminosity and, in particular, the 24  and the 33 GHz data, we conclude that, by combining in a consistent way this relevant information, it is possible to obtain a fairly robust estimate of the upper mass limit of the IMF . Indeed we stress that the FIR luminosity and those at H, 33 GHz and 24  bracket, on one side, the intrinsic ionizing photon flux and, on the other, the attenuation and re-emission by molecular clouds where the stars reside during their initial phase, which are the most sensitive to the fractional contribution of the most massive stars.

As far as the prototype normal galaxy M100 is concerned, we remind that the observed H luminosity was matched by the model with . It is important, at this stage, to remind that these conclusions should be taken with care because the current SSP adopted do not account for binary evolution. Indeed it is known that binary evolution can produce more ionizing photons at later stages than single star evolution, because of the loss of the envelope by binary interaction (Wofford et al., 2016). In particular we may say that regions 1 ,2, 3, 5, 6 may be reproduced fairly well without invoking a high and/or significant effects by binary evolution. The latter effect could instead explain the location of region number 8 that cannot be reproduced even with , as indicated by the green diamond in Figure 12). This effect, which add a new dimension to the problem, is under study in PARSEC and it will be included in the nearest future.

The q ratio.
The second thing that we note by inspecting Table 3 is that the value of the q does not significantly depend on the upper limit of the IMF. This is because a significant fraction of both the bolometric luminosity and the radio luminosity are contributed by stars less massive than M = . Indeed after a few Myr, the upper mass of the SSP is already significantly lower than M = . Furthermore, the q ratio is generally higher in the extranuclear regions of NGC 6946 than in the normal SF galaxy M100, by a factor of 1.4. We remind that in our model, CCSNs, responsible for the non thermal emission, are produced only by the mass range M   and the individual star-bursting regions have an age which is younger than the threshold for CCSN production. Thus they did not yet reach a stationary equilibrium for non-thermal radio emission and not only their radio slopes appear flatter but also their q are higher than that of normal star-forming galaxies B02.

The SFR calibration.
Using the luminosity in the different bands, as presented in Table 3, and the average SFR obtained from the GRASIL fits, we show in Table 4, the corresponding SFR calibrations C(b) = SFR/L(b), where b indicates the band. The average SFR is derived by considering the last 100 Myr time interval for the normal star-forming galaxy M100, and the age of the burst in the star-forming regions of NGC 6946, which is shown in column 16 of Table 4 for the NGC 6946 extranuclear regions. This choice is suggested by the fact that within an entire galaxy, the star-forming regions distribute continuously with time while, this is obviously not the case for individual regions. The average SFR () so computed is shown in the second column of Table 4. These calibrations are directly obtained by GRASIL. For sake of comparison, we also add, in Table 4, some calibrations that can be obtained by using the un-attenuated fluxes derived from GRASIL best fit models, either directly or by using the analytical relations obtained by means of the SSPs. These are useful to compare the differences between the SFR derived from a panchromatic fit that realistically account for ISM processes to analytical fits made with SSP fluxes. For example, the FUV, H and H values shown in column 3, 4 and 6, refer to the calibrations obtained by dividing the average SFR by the corresponding un-attenuated intrinsic fluxes of the models, listed in Table 3. For H and H we also show, in columns 5 and 7, the values obtained by inserting their corresponding intrinsic fluxes in the analytical calibrations obtained from simple SSPs models, Equations 40 and 39 respectively. For the 33 GHz radio emission, we add in column 15, the value of the calibration C(33) = SFR/L(33 GHz), directly derived from the SSP models (Equation 38). The last rows in the first and second sections, indicated by , show the median values of all SFR calibrations derived from the star-bursting regions of NGC 6946. In the third section we collect some common SFR calibrations from the literature, with the corresponding bibliographic sources given in the caption.

Figure 12: Left panel: SFR calibration constants C(), in units of , at  = 1.4, 4.9, 8.5 and 33.0 (Table 4), for the star-bursting regions and using the total radio emission. Triangles and daimonds indicate cases of and respectively. Right panel: The same as in the left panel but for M 100 and using both the total radio (crosses) and non-thermal radio emission (xs). In both panels, all red symbols and lines indicate cases of and blue symbols and lines cases of . Also all lines represent least-square fits.

In Figure 12 we plot the calibration constants C() at  = 1.4, 4.9, 8.5 and 33.0 , in units of , against the corresponding radio frequencies, for both the star-bursting regions (left panel) and M100 (right panel) and for both cases of . The relations (in logarithmic units) are almost linear above 1.4 GHz, and we show the best fit regressions as solid lines. Concerning the dependence on M, there is a clear difference between the starburst and the normal star formation regime in the sense that, in the former case the constant is about a factor of three larger independently from the frequency. Again this effect is due to the young age of the starburst regions coupled with the dependence of the thermal radio emission on M. Indeed for the normal SFR regime, the dependence on M becomes significant only at high frequencies, where the thermal emission start dominating. From the fits we derive the following general relations between SFR (in ) and radio luminosity (in ). For Star-bursting regions, with average metallicity 0.006, we have at any frequency between 1.4  and 33 , to better than 10 per cent,


while, for the normal star-forming galaxy, M100, with slightly more than solar metallicity, we have for the total radio emission


and for the non-thermal radio emission, we have


Relations similar to Equations 22,  23 and  24 above were provided by Murphy et al. (2011), assuming Kroupa IMF (with ), solar metallicity and a constant SFR over a timescale of 100 Myr. For total radio emission


and for non-thermal radio emission


These relations are in between our relations in 22 and 23, running only  20 per cent lower than our one for normal star-forming galaxies, below 33 GHz. However at frequencies higher than 33 GHz the difference grows rapidly and at 100 GHz our values are 40 per cent and 200 per cent higher than that of Murphy et al. (2011) for starburst and normal galaxies, respectively. Interestingly, the relation derived by Schmitt et al. (2006) with Starburst99, assuming a Salpeter IMF (with ), solar metallicity and a continuous SFR of 1


agrees well with our relations in 22 and 23 at frequencies from about 33 GHz to 100 GHz. Above 100 GHz the dust emission contribution becomes dominant.

It is interesting to compare the GRASIL calibration at 33 GHz with that derived from the SSP models (Equation 38), i.e. columns 14 and column 15 in Table 4. The values show some discrepancies, especially for normal galaxies. However we remind that the values in column 15 refer to free-free emission alone while the 33 GHz flux in column 14 may include a non negligible contribution by synchrotron emission. If we take this contamination into account, using the percentage contribution of the thermal radio component enclosed in parenthesis, the corrected values of column 14 are in fair agreement with those predicted by Equation 38. This is expected because, at these wavelengths, attenuation, including also the free-free one (Vega et al., 2008), should not be important. The above overall agreement indicate a surprising robustness of the radio calibration, irrespective of the differences in the underlying models. This will be particularly important for the high redshift galaxies especially in their early evolutionary phases.

As a further comparison, we provide at the bottom of Table 4 other SFR calibrations taken from the literature. We remind the reader that, in this work, we adopt a Kennicutt (1983) IMF and a typical time-scale for averaging the SFR rate of 100 Myr for normal galaxies and equal to the age of the starburst in the extra-nuclear star bursting regions. The details on the parameters adopted by different authors may be found in the quoted papers.

5.2 Dust attenuation properties

It has been already shown in (Silva et al., 1998; Granato et al., 2000) that the attenuation is the result of the interplay between the extinction properties of dust grains and the distribution between dust and stars in space and in time. Further discussion can also be found in (Charlot & Fall, 2000; Panuzzo et al., 2003). We now investigate the properties of the attenuation curves of the galaxies studied here, with particular interest to understand its dependence on the galaxy type, i.e. starburst vs. normal regime, and also the effects of using different parameters in the analysis, such as M in the IMF.

We first show in Figure 13, the attenuation curves, A/A, of M100 (green solid line) and of the NGC 6946 extra-nuclear regions studied in this work (black solid lines). They are simply derived by comparing the transmitted galaxy flux to the intrinsic one, as a function of the wavelength. The top panel refers to the case of while the case with is shown in the bottom panel.

Figure 13: Attenuation curves, A/A, for M100 and NGC 6946 extra-nuclear regions, for the cases of (top panel) and (bottom panel). The black solid lines are the attenuations curve of the eight regions of NGC 6946 while the dotted blue line represents the median of these curves. The median values of and for the star-forming regions are also given in the labels. The green line is the attenuations curve of M100. For comparison purposes, we added the attenuation curves of Calzetti et al. (2000) (red dashed line) and Cardelli, Clayton & Mathis (1989) (purple dotted-dashed line) The later attenuation seems to agree very well with that of M100 at wavelengths above 4000 Å.

The dotted blue lines represent the median curve of the NGC 6946 extra-nuclear regions obtained by sampling the curves of individual regions at selected wavelengths. The median values of the quantities and R = A/E(B-V) for the star-bursting regions are also given in the labels. We have also labelled some prominent emission lines ( line is in brown). For comparison purposes, we added the attenuation curves of Calzetti et al. (2000) (red dashed lines) and Cardelli, Clayton & Mathis (1989) (purple dotted-dashed lines). We also collect in Columns (2 - 9) of Table 5 the attenuation in H, H, in the FUV (0.16 ) and NUV (0.20 ), in the B-band(0.45 ), and in the interpolated continuum at H (named A48), in the V-band(0.55 ) and in the interpolated continuum at H (named A65). Columns 10 and 11 give the optical depth of the molecular cloud at 1  derived from GRASIL and the ratio, respectively. In Column 12 we show the difference A-A which corresponds to . Using the Balmer decrement method for the and lines, we also obtained the same quantity, :


and show its value in Column 13. We assume an intrinsic line intensity ratio,  = 2.86, which is valid for =10,000 K (Hummer & Storey, 1987). Column 14 gives the escape time of young stars from their birth cloud in Myr. Column 15 gives the FIR flux in in the 40 - 120  interval derived using the 60 and 120  fluxes (Helou et al., 1988).

Several points can be noted by inspection of Figure 13 and Table 5. First of all, we see that, while in normal galaxies the attenuation in the lines is significantly higher than that in the surrounding continuum, (a factor of 4 on average for M100), this is not true for the star bursting regions where this factor is not more than 1.3. The fact that the attenuation in the lines is significantly higher than that in the surrounding continuum is well known and it is a manifestation that young stars are more dust enshrouded than older stars. However, in the case of the star-bursting regions the underlying continuum is essentially produced by the same stellar populations, with a small contribution from the older populations, so that the attenuation in the lines is not much different from that in the continuum. To get more insights from the models we compare the characteristic R values in the lines and in the nearby continuum. We see that for star-forming galaxies, R = AA-A) is much larger than the value obtained in the corresponding interpolated underlying continuum regions (R). For example, for M100, R = 16(10) while R = 4(4), for M = 40(120) M. In the starburst regions, R and R have similar values between 6 and 10. The high values so determined for R would suggest a high neutral absorption in these objects but this is not the case because we have not modified the grain size distribution which, by default in GRASIL is the one of the Galactic ISM. We also see that in M100, R, in the continuum, has a value of 3, compatible with the common extinction laws provided in literature (Calzetti et al., 2000; Cardelli, Clayton & Mathis, 1989). Instead this effect is due to the fact that the attenuation within molecular clouds is so high (see the values of in Table 5) that the emission of young stars is almost completely absorbed, until they escape from the clouds. This effect mimics the neutral absorption. This effect is enhanced for the emission lines because they are generated only by stars younger than about 6 Myr, so that the fraction of line flux absorbed by molecular clouds with respect to the total one, is even higher. Adopting a larger than 7 Myr will erase line emission in our model, irrespective of the galaxy types considered here. A major consequence of the age selective attenuation described above is that it leads to wrong estimates of the attenuation in the line emission, if one uses methods involving attenuation-affected observables, like the Balmer decrement method. A more accurate method to estimate the intrinsic H luminosity (hence the H attenuation) is by using the 33 GHz and 24  luminosities, since they are optimal tracers of the radiation emitted by the most massive. Adopting a median value of 78.5 per cent as the percentage contribution of the thermal radio component in young star-bursts, we first obtain the thermal radio emission component of our model’s 33 GHz total radio emission luminosity. Using the analytical relation given in Equation 2 (T and gaunt factors used here are as given in Equations 34 and  3 respectively), we derive the number of ionizing photons corresponding to this luminosity. With Equation 35, we then estimate the corresponding intrinsic H luminosity. A plot of the estimated intrinsic H luminosity (in ) against our model’s 33 GHz thermal radio luminosity (L in ) is presented in Figure 14. We added in this figure also the normal SF galaxy M100, with an average thermal fraction of 16.0 per cent(42.0 per cent) for M = 40 (M = 120 ), respectively. The fitting relation (solid line) considering only the star bursting regions for both cases of M, is given by


Figure 14: Plot of our model’s 33 GHz thermal radio luminosity against the intrinsic H luminosity. Solid line represents linear fit for the star-bursting regions. The H luminosity here is derived from the thermal radio luminosity using Equations 2 and 35. We adopted a median thermal fraction of 78.5 per cent in estimating the thermal radio component of the 33 GHz total radio luminosity. We also added in this plot the M100 which has an average thermal fraction of 16.0 per cent(42.0 per cent) for M  =  40 (M  =  120 ).

luminosity (L33 in erg s?1Hz?1) is presented in Figure 12.

The corresponding attenuation at H can be obtained by using the intrinsic value provided by equation 29


As mentioned earlier, this attenuation is larger than the one that can be obtained with the Balmer decrement method because the latter might not accounts for the strong obscuration in the early life of massive stars (see Table 5). Thus, in presence of an estimate of the radio flux at 33 GHz, Equation 30 should be preferred for young starbursts. A similar relation, though with somewhat larger scatter, can be obtained for the 24  flux, which is shown in Figure 15. The relation between the model’s intrinsic H (in ) and 24  luminositiy, L = L in , is


almost independently from M. Correspondingly, the derived attenuation in H is

Figure 15: Plot of the 24  luminosity against the intrinsic H luminosity, both extracted directly from our best-fit models. Solid line represents linear fit. We excluded M100 in this plot so as to better represent the star-bursting regions.

It is also interesting to provide a relation between the attenuation at 1600 Å (Column 4 of Table 5) and the ratio of FIR (Column 11 of Table 3) to F (Column 3 of Table 3). The plot is shown in Figure 16, along with the linear fitting relation obtained for the star-bursting regions (solid line), for both M cases.


We note that the normal galaxy M100 is out from the relation. This may likely be due to the non negligible contribution to the FIR intermediate age stellar populations that are not contributing to the far UV.

Figure 16: Plot of attenuation at 1600 Å against the ratio, , both extracted from our best-fit models. The solid black line represents our linear fit as given in Equation 33. The red crosses and blue ’X’s indicate M and 120  cases respectively.

6 Conclusions

In this paper we have used the recent PARSEC tracks to compute the integrated stellar light, the ionizing photon budget and the supernova rates predicted by young SSP models for different IMF upper mass limits. The SSPs spans a wide range in metallicities, from 0.0001 to 0.04 and initial masses ranging from 0.1 to 350M. In building the integrated spectra, we have adopted the spectral library compiled by Chen et al. (2015). This library is a result of homogenising different sets of libraries through a process fully described in Girardi et al. (2002) and Chen et al. (2014). Using this integrated spectra in the photoionization code CLOUDY, we have also calculated and included in the integrated spectra the nebular continum and the intensities of some selected emission lines. With the new SSPs we can then be able to predict the panchromatic spectrum and the main recombination lines of star-forming galaxies of any type, with GRASIL.

We made also two major revisions in the radio emission model. First, we revisited the relation between free-free radio emission, number of ionizing photons and radio frequency and came up with a relation that (a) takes into account the full dependence of the electron temperature on the metallicity and the effect of considering different IMF upper mass limits and (b) incorporated a more accurate gaunt factor term. The free-free radio emission shows a significant variation with metallicity and IMF upper mass limit, decreasing by a factor of about 3 when Z increases from 0.0001 to 0.02 and inceasing by more than an order of of magnitude when the mass limit is increased from 40 to 350 , at constant total mass. These differences cannot be negelected, especially when calibrating SFRs using tracers that depend strongly on the ionizing photon budget. Hence, the corresponding SFR calibrations (H, H and thermal radio emission) take into account this dependence on metallicity and IMF upper limits.

Second, we revised the non-thermal radio emission model originally described by B02, taking into account recent advances in the CCSN explosion mechanisms which indicates a range of stellar masses where the stars fail to explode, the so called Failed SNe. We adopt a threshold value of 30 M, above which stars do not produce non-thermal radio emission. This has two immediate effects: (a) the beginning of the non-thermal radio emission is delayed by about 7 Myr, the lifetime of a 30 M star; (b) the non-thermal radio emission is independent from the IMF upper mass limit as long as it is above 30 M. It is also assumed that pair-instability SNe, discussed by Slemer et al. (2017) do not produce synchrotron radiation.

With these revisions in the new SSPs and in the radio emission model, there is the need to check some parameters of the GRASIL code. We first re-calibrate the proportionality constant between the supernova rate and the corresponding non-thermal radio luminosity, ), using one of the best-sampled nearby normal star-forming galaxy, M100. We are able to reproduce very well the far-UV to radio SED of this galaxy, for IMF upper mass limits of 40  and 120 , with an average of 2.42. We obtain a value of = 1.94 which is larger than that obtained previously by B02 by a factor of 1.35, but it is similar to that obtained by Vega et al. (2005). We find that using M produces a thermal radio emission that significantly exceeds the observed one. This excess thermal radio emission cannot be cured by varying any other parameter in the fit. This result suggests that, for a normal star-forming galaxy, it is difficult to have an average IMF extending up to such high initial masses. What depends strongly on the upper mass limit of the IMF is the number of ionizing photons, which beqars on the intensity of the recombination lines and on the thermal radio emission. Indeed, in spite of being able to reproduce the FIR and 1.4 GHz radio emission, the model that adopts M = 120  over predicts the H emission by a factor of about three while, in the case of M = 40 , the discrepancy is only of a few percent.

We then check the new thermal radio emission model against the well studied thermal radio dominated star-forming regions in NGC 6946. Even in this case we are able to reproduce very well the observed SEDs from NIR to radio wavelengths for both cases of the IMF upper mass limits. In fitting the SEDs of these thermal radio-dominanted star-bursting regions, we adopt the value of resulting from the non-thermal radio calibration with the normal star-forming galaxy M100. The estimated values of the ratio in these regions lie between 2.5 and 2.6, implying a relatively lower non-thermal emission than in normal star-forming galaxies. The additional evidence of flatter radio slopes supports the notion that there is lack of the non-thermal emission as predicted by the original non-thermal radio emission models by B02. The resulting ages of the bursts, range from 7 to 12 Myr, confirming that these observations can be used to determine the star-burst ages. The fit obtained with the two IMF upper mass limits exhibit interesting differences, in particular in the predicted H and UV luminosities which are the most sensitive to the IMF and to dust attenuations. We show that, by combining information from the FIR, 24 , 33 GHz and H, we can determine the preferred value for M. This is the first time, to the best of our knowledge, that the IMF upper mass limit can be consistently determined. However, for region 8 we cannot reproduce the observed H flux even with M. This could be a case where binary evolution, which is known to produce an excess of ionizing photons due to the loss of the envelope of massive stars by binary interaction, may be required.

With luminosities and averaged SFR extracted from our best fit models (Table 3), we derive SFR calibrations in different bands for both normal galaxies and star-burst regions (Table 4). Since thermal radio emission and that in the recombination lines depend strongly on the upper mass limit of the IMF, we provide multiple regression fitting relations with M and metallicity (e.g. Equations 38, 39 40).

Finally, exploiting the realistic treatment of dust performed by GRASIL, we investigate the properties of the attenuation curves of the galaxies studied in this work with the aim of understanding their dependence on the galaxy type. We find that, while in the normal SF galaxy M100 the attenuation in the lines is significantly higher than that in the surrounding continuum, for the star bursting regions of NGC 6946 the two attenuations are similar. A common property is that the predicted R values obtained for the line is large, mimicking a significant neutral absorption. For star bursting regions, large values are found also for the continuum while, in the case of M100, the R value of the continuum is normal. Large R values have been found in high redshift star-forming galaxies by Fan et al. (2014), using a completely different population synthesis tool. A major disturbing consequence of the age selective attenuation, which could not be revealed in a foreground screen model, is that it leads to wrong estimates of attenuation when using methods involving observed line fluxes, like the Balmer decrement. Of course, more accurate methods are those that combine UV or line fluxes with fluxes that are not affected by attenuation, such as the well known relation between the FUV (at 1600 Å) attenuation and the flux ratio ratio, FIR/F. For the star-forming regions we revisit the former relation (Equation 33) and we provide new relations between the H attenuation and the observed Hand 24  (Equation 30) or 33 GHz (Equation 32) fluxes. The above mentioned relations, which we show to be almost independent from M, can be extremely useful in estimating the attenuations in young high redshift galaxies. For this purpose we are working on a larger galaxy sample to increase the statistical significance of our results.


We thank ?? for helpful discussions. A. Bressan and L. Girardi acknowledge financial support from INAF through grants PRIN-INAF-2014-14. P. Marigo and Y. Chen acknowledge support from ERC Consolidator Grant “STARKEY”, G.A. n. 615604. F. Perrotta was supported by the RADIOFOREGROUNDS grant of the European Union’s Horizon 2020 research and innovation programme COMPET-05-2015, grant agreement number 687312.

Appendix A Variation of electron temperature and the constants, and , with IMF Upper mass limits and metallicity

(1) (2) (3) (4) (5)
0.0200 0.554 18.03 0.64 1.47
0.0080 0.738 13.56 0.89 1.96
0.0040 0.943 10.61 0.11 2.51
0.0005 1.376 7.27 0.16 3.66
0.0001 1.823 5.48 0.17 4.85
0.0200 1.730 5.78 0.65 4.60
0.0080 2.153 4.64 0.94 5.73
0.0040 2.520 3.97 0.12 6.70
0.0005 3.065 3.26 0.17 8.15
0.0001 3.761 2.66 0.19 10.00
0.0200 2.819 3.54 0.65 7.50
0.0080 3.469 2.88 0.95 9.22
0.0040 3.739 2.67 0.12 9.94
0.0005 4.269 2.34 0.17 11.35
0.0001 5.144 1.94 0.20 13.68

Column(1): metallicty. Column(2): production rate of ionizing photons per unit solar mass formed. Column(3): calibration constant () of the SFR-IQ(H) relation as given in equation 8. Column(4): electron temperature in K computed using CLOUDY. Column(5): calibration constant () of the - relation as given in equation 11.

Table 2: Calibration constants, and , for different metallicities and upper mass limits.
Figure 17: Plot of electron temperature T (obtained from CLOUDY) against metallicity Z for  = 40, 120 and 350 indicated by the red ’X’s, blue crosses and green xterisks respectively. The lower axis is the metallicity while the upper one is the corresponding oxygen abundance, . The blue dashed-dotted line indicates the metallicity at Z = 0.0134 which we adopted in this work as the solar metallicity. The average of the empirical fits derived by López-Sánchez et al. (2012) for high-ionization \ionOIII and for low-ionization \ionOII zones is shown by the black dashed line.
Figure 18: Variation of the constant with age at different metallicities (Z = 0.0001, 0.0005, 0.004, 0.008 and 0.02) and different upper mass limits. The symbols used to indicate the different