Combining physical galaxy models with radio observations to constrain the SFRs of high-z dusty star-forming galaxies

Abstract

1 INTRODUCTION

In star-forming galaxies, stars, gas and dust are mixed in a very complicated way, and dust obscuration strongly depends on their relative geometrical distribution. This is particularly related to the ages of stellar populations, as young stars embedded in dense molecular clouds (MCs) are more extinguished than older stars (e.g. Calzetti, Kinney & Storchi-Bergmann 1994; Silva et al. 1998; Charlot & Fall 2000; Poggianti, Bressan & Franceschini 2001).

New data combining deep optical and far-infrared (FIR) photometric imaging, as discussed in Lo Faro et al. (2013, BLF13 hereafter) among many others, then force us to consider relatively realistic models of galaxy synthetic spectra, with radical complications with respect to previous modelling, occurring at two levels. On one side, dust extinction cannot be neglected and should be considered as a function of the age of the stellar populations in the galaxy. In addition, the differential extinction in age and dusty environments entails that geometrical effects in the distribution of stars and dust play a fundamental role in determining the UV to submillimetre (sub-mm) spectral energy distribution (SED), and have to be carefully modelled.

By analysing the far-UV to sub-mm properties of a sample of 31 dusty star-forming luminous and ultraluminous infrared galaxies [(U)LIRGs] at z ∼ 1–2, BLF13 have in fact demonstrated that very idealized approaches, e.g. based on the assumption of a homogeneous foreground screen of dust, single extinction law, and optical-only SED-fitting procedures, may produce highly degenerate model solutions and galaxy SEDs unable to energetically balance the dust reprocessed infrared (IR) emission from the galaxy. For the most dust obscured objects, as observed by BLF13, this can result in severely underestimated stellar masses. Another crucial aspect of the BLF13 analysis concerns the estimates of the star formation rates (SFRs) for high-z dusty star-forming galaxies. The classical Kennicutt (1998, K98 hereafter) calibration, widely used in literature to estimate SFRs from the rest-frame 8–1000 μm total IR (TIR) luminosity, is based on the assumption that the bolometric luminosity produced by a 100 Myr constant star formation (SF) is all emitted in the IR. The BLF13 typical galaxies, instead, appeared to include significant contributions to the dust heating by intermediate-age stellar populations, i.e. older than the typical escape time of young stars from their parent MCs (in the range between ∼3 and ∼90 Myr), and powering the cirrus emission. This brings to a factor of ∼2 lower calibration on average of SFR from the total L_IR, for moderately star-forming galaxies, with respect to the Kennicutt's calibration. Therefore, although the FIR emission is surely related to the recent SFR in a galaxy, it also includes contributions by all populations of stars heating the dust (i.e. the cirrus component; Helou 1986; Lonsdale Persson & Helou 1987; Hirashita, Buat & Inoue 2003; Bendo et al. 2010; Li et al. 2010), as actually confirmed by our analysis. The idea that some or most of the FIR and sub-mm emission of sub-mm galaxies may be due to cirrus was first discussed by Efstathiou & Rowan-Robinson (2003, ERR03 hereafter). ERR03 actually also used the radio data and the FIR/radio correlation in the fitting. Efstathiou & Siebenmorgen (2009) explored this idea further by modelling a sample of sub-mm galaxies with Spitzer spectroscopy and redshifts and confirmed the idea that cirrus is contributing significantly in the FIR. In addition, a fraction of the UV-optical emission escapes the galaxy and is not registered by the IR luminosity.

The estimates of masses and SFRs for high-z galaxies provide fundamental constraints for models of galaxy formation and evolution. Notably, these quantities bear on the position of galaxies on the SFR versus stellar mass plane, with galaxies on the main sequence (MS; e.g. Rodighiero et al. 2011) being interpreted as steadily evolving, in contrast to outliers (off-MS objects) probably undergoing starburst episodes. Therefore, further independent tests of the star formation histories (SFHs) and recent SFR are then required in order to check our interpretation of high-z galaxy SEDs particularly for what concerns the contribution of intermediate-age stars, affecting both the estimated stellar mass and recent SFR. Especially suited to this aim is the radio spectral range.

The radio emission from normal star-forming galaxies is usually dominated by the non-thermal (NT) component (up to ∼90 per cent of the radio flux) which is due to the synchrotron emission from relativistic electrons accelerated into the shocked interstellar medium, by core-collapse supernova (CCSN) explosions (Condon 1992). When considering high-redshift highly star-forming galaxies and starbursts usually the contribution from free–free emission becomes more important and dominates at ∼30 GHz rest frame. By looking for spectral slope variations in a sample of sub-mm galaxies at z ∼ 2–3, Ibar et al. (2010) found the same α ∼ 0.7 between 610 MHz and 1.4 GHz emphasizing that apparently the optically thin synchrotron emission was still dominant at those redshifts.

Only stars more massive than ∼8 M_⊙ produce the galactic cosmic rays responsible for the NT emission. These stars ionize the H ii regions as well contributing to the thermal component. Radio emission is therefore a probe of the recent SF activity in normal galaxies.

The inclusion of radio data into our SED-fitting procedure is therefore crucial as the radio continuum offers a further independent way to estimate and constrain the SFR in galaxies completely unaffected by extinction. Moreover, as the radio emission probes the SFH on a different time-scale with respect to the IR emission, this analysis is also useful for our understanding of the SFH of galaxies. At the same time, including new data directly related to the rate of formation of new stars can help to better constrain the average dust extinction, by comparison with the FIR and UV spectral information.

Despite the strong effort devoted to comparing the different SFR indicators from X-ray to radio (see e.g. Daddi et al. 2007, 2010; Murphy et al. 2011; Kurczynski et al. 2012 among the others) and probing the FIR–radio correlation up to high redshift (e.g. Ivison et al. 2010; Sargent et al. 2010; Mao et al. 2011; Pannella et al. 2014) there are only few works including radio emission modelling in galaxy evolution synthesis with the aim of unveiling the nature of different galaxy populations. Most of them are based on the more classical and semi-empirical multicomponent SED-fitting procedure where the different parts of the SED (optical–near-infrared (NIR, mid-infrared (MIR)–FIR, radio) are modelled individually assuming for the radio the FIR–radio correlation and with the main aim of disentangling the AGN contribution from that of SF dominated sources (see e.g. U et al. 2012).

The work by Bressan, Silva & Granato (2002) was mainly focused on understanding the nature of the FIR–radio correlation, particularly for a local sample of obscured starbursts, by making use of state-of-the-art models of star-forming galaxies. They used grasil to model the IR emission of galaxies and extended its predictions beyond the sub-mm regime, up to 316 MHz. They showed that the delay of the NT emission with respect to the IR and thermal radio can give rise to observable effects in the IR to radio ratio as well as in the radio spectral slope, potentially allowing the analysis of obscured starbursts with a time resolution of a few tens of Myr, unreachable with other SF indicators.

Vega et al. (2008) also included radio modelling (with at least three data points in the radio regime) in their analysis of a sample of 30 local (U)LIRGs, with the aim of studying their starburst nature and disentangling the possible AGN contribution. The well-sampled radio spectra allowed them to put strong constraints on the age of the burst of SF.

This work represents the first attempt to apply this sophisticated physical analysis, including the modelling of radio emission, to well-sampled (from far-UV to radio) SEDs of high-z dust obscured SF galaxies. We investigate the effect of accounting for radio data in the SED-fitting procedure on the derived main physical properties of galaxies, focusing in particular on the constraints on the current SFR, SFH and consequently M_⋆ of galaxies. Our analysis points towards the radio flux as an essential information for interpreting star-forming galaxies at high redshifts and for recovering reliable SFH.

With the above purposes, in this work we have complemented the MIR-selected sample studied in BLF13 by including into our SED analysis new Very Large Array (VLA¹) observations of the integrated 1.4 GHz emission of galaxies in the GOODS-South field. Both because of the small number of (U)LIRGs in our sample with high signal-to-noise (S/N) radio data (4/31), and to extend our analysis also to a different family of galaxies, we have included in our study the six well-sampled NIR-selected BzK galaxies by Daddi et al. (2010) at z ∼ 1.5 in the GOODS-North field. The best-fitting model for each observed galaxy is searched within a grasil-generated library of ∼10⁶ models, covering a huge range of SFHs and dust model parameters, with model SEDs covering self-consistently the UV to radio spectral range. The FIR and radio luminosities as SF indicators are then discussed by making use of the best-fitting SFRs from our thorough SED modelling, and by comparing to the most widely adopted relations by K98, derived from stellar population modelling, and Bell et al. (2003), based on the FIR–radio correlation. Finally, we discuss the consequences of our inferred values of stellar mass and SFR on the MS relation.

The paper is organized as follows. In Section 2, we describe the (U)LIRGs and BzK data sample we use for our analysis, including the new VLA data. The physical model and SED-fitting procedure are presented in Section 3. The results concerning the best-fitting far-UV to radio SEDs are discussed in Section 4, while we compare with the observed FIR–radio relation in Section 5. In Section 6, the radio and FIR as SFR indicators are discussed, and the implications of our models for the SFR–stellar mass relation are presented in Section 7. Our summary and conclusions are presented in Section 8.

The models adopt a 0.1–100 M_⊙ Salpeter (1955) initial mass function (IMF). Where necessary, and indicated in the text, we scale to a Chabrier (2003) IMF by dividing masses by a factor of 1.7. The cosmological parameters we adopt assume H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7 and Ω_M = 0.3.

2 DATA SAMPLE AND OBSERVATIONS

This work is intended to be complementary to the physical analysis performed in our previous paper (BLF13) and aimed at testing BLF13 estimates, in particular those concerning M_⋆ and SFR, against observations through the investigation of the radio emission of high-z star-forming galaxies. For this reason, we first apply our analysis to the same sample of 31 high-z (U)LIRGs whose far-UV to sub-mm properties have been deeply investigated in BLF13 and then, in order to strengthen our conclusions, we extend the analysis by including a small but well representative sample of six BzK-selected star-forming MS galaxies at z ∼ 1.5 for which also high S/N radio observations are available. These two samples benefit from a full multiwavelength coverage from far-UV to radio including Spitzer and Herschel [Photoconductor Array Camera and Spectrometer (PACS) + Spectral and Photometric Imaging Receiver (SPIRE)] data. The two samples are briefly presented in Sections 2.1 and 2.2.

2.1 The (U)LIRG data sample

The sample of high-z (U)LIRGs whose radio properties are investigated here was originally selected from the sample of IR-luminous galaxies in GOODS-S presented by Fadda et al. (2010, F10 hereafter). It includes the faintest 24 μm sources observed with the Spitzer Infrared Spectrograph (IRS; S₂₄ ∼ 0.15–0.45 mJy) in the two redshift bins (0.76–1.05 and 1.75–2.4). These galaxies have been selected by F10 specifically targeting LIRGs at z ∼ 1 and ULIRGs at z ∼ 2. This sample is therefore crudely luminosity selected and F10 did not apply any other selections.

As emphasized by F10, except for having higher dust obscuration, these galaxies do not have extremely deviant properties in the rest-frame UV/optical compared to galaxies selected at observed optical/NIR band. Their observed optical/NIR colours are very similar to those of extremely red galaxy populations selected by large area K-band surveys. Moreover, all the z ∼ 2 ULIRGs are IR-excess BzK [(z − B) versus (K − z); Daddi et al. (2004)] galaxies and most of them have L_FIR/L_1600 Å ratios higher than those of starburst galaxies at a given UV slope. The ‘IR excess’ is mainly due to the strong polycyclic aromatic hydrocarbon (PAH) contribution to the MIR luminosity which leads to overestimation of L_IR and to underestimation of the UV dust extinction.

The sample studied here includes only galaxies powered by SF. All the objects which were classified by F10 as AGN dominated, on the basis of several indicators such as broad and high ionization lines in optical spectra, lack a 1.6 μm stellar bump in the SED, X-ray bright sources, low MIR 6.2 μm Equivalent Width (EW) and optical morphology (see also Pozzi et al. 2012), were excluded by BLF13 from the original sample. With the further requirement of a full complementary optical/NIR photometry provided by the multiwavelength MUSIC catalogue of Santini et al. (2009) the final sample consists of 31 (U)LIRGs among which 10 are at z ∼ 1 (all in the nominal LIRG regime) and 21, mostly ULIRGs, are at z ∼ 2.

All the 31 sources in the GOODS-S field have the advantage of a full FIR coverage from 70 to 500 μm from both Herschel SPIRE and PACS instruments with the SPIRE and PACS data taken, respectively, from the Herschel Multitiered Extragalactic Survey (HerMES; Oliver et al. 2012) and the PACS Evolutionary Probe (PEP; Lutz et al. 2011) programmes. Typical noise levels are ∼1 mJy for PACS 70–160 μm and ∼6 mJy for SPIRE 250–500 μm, including confusion.

2.1.1 Radio observations

In addition to the availability of a spectroscopic redshift and complete FIR coverage we also looked for VLA detections at 1.4 GHz. As already pointed out above, the inclusion of radio emission into our procedure has been revealed to be fundamental to constrain the ‘current’ SFR in galaxies, mostly contributed by very young stars, and their SFH.

The publicly available VLA catalogue in GOODS-S of Miller et al. (2013) includes only 3 of the 31 objects of our (U)LIRG sample. Inclusion in this catalogue required a point-source detection of >5σ, a fairly conservative but well-established limit for radio catalogues. We have extracted radio flux densities for the remaining sources directly from the released radio mosaic image, adopting the flux density at the position of the (U)LIRG coordinates. This measurement, in units of μJy per beam, is a reasonable representation of the flux density of any point source at the location corresponding to the (U)LIRG coordinates. The majority of the (U)LIRGs are associated with positive flux densities at lower S/N than required for the formal published catalogue, as seen in Fig. 1. In fact, we find one additional source for which S/N >4σ and we include this source in our detections presented in Table 1 which contains the coordinates of the (U)LIRGs, the measured flux densities in μJy at those coordinates, and the rms of the radio mosaic at those coordinates.

Figure 1.

Histogram showing the measured radio flux densities at the indicated positions of (U)LIRGs. It clearly shows that the expectation of a Gaussian centred on zero and with a dispersion of the image noise (red histogram) is very different from the actual measured distribution.

Open in new tabDownload slide

As a variation in the histogram shown in Fig. 1, we have also performed a radio stacking analysis on the 31 (U)LIRGs. This involved taking cutouts of the radio mosaic image centred on each of the 31 (U)LIRG coordinates and combining them to evaluate a single statistical representation of the full population. When performing the stack using an average with rejection of the single highest and single lowest value (‘minmax’ rejection), we recovered a significant detected source having a peak flux density of 11.6 μJy and an rms of 1.2 μJy. A stack produced using a straight median instead of the average yielded consistent results. To check the validity of this stack detection, we created a list of 1000 random positions in the vicinity of the (U)LIRG coordinates and evaluated a stack based on these coordinates. At the central pixel of this stack, the measured flux density was 8 pJy beam⁻¹ (i.e. effectively zero) and the rms in the stack image was 0.2 μJy. Thus, we are confident that the (U)LIRG population represents real radio emission with flux densities only occasionally above the threshold used in the formal Miller et al. (2013) catalogue.

When modelling the radio data we have used the following prescriptions to deal with the errors: (i) for the objects having S_1.4 GHz/rms ≳ 3σ we have used the same error as given in the catalogue, usually of the order of 6.1–6.3 μJy; (ii) for those objects having instead 2σ ≲ S_1.4 GHz/rms ≲ 3σ we have assumed a 2σ error bar associated with the measured flux density; (iii) finally we consider as 3σ upper limits all the sources having negative average flux densities measured at the position of the source and/or S_1.4 GHz/rms ≲ 2σ.

We note that the presence of a single data point in the radio part of the spectrum compared to the larger number of data available in the far-UV to FIR regime does not affect significantly the determination of the best-fitting model in terms of its χ² value but it plays a crucial role in constraining the SFH of the galaxy. What we are interested in, in fact, is determining if our best-fitting solutions are able to reproduce, within the uncertainties, the observed radio flux densities. Thus, we intend to verify if the amount of massive young stellar populations dominating the most recent SFR is well constrained by our model as well as the fraction of stars outside the MCs contributing to cirrus heating.

2.2 BzK data sample

The six BzK galaxies analysed here have been selected from the original sample studied in Daddi et al. (2010, D10 hereafter). These are NIR-selected galaxies, in GOODS-N, with K < 20.5 (Vega scale; or K < 22.37 AB), to which the BzK colour selection criterion of Daddi et al. (2004) has been applied, together with the requirement of a detection in deep Spitzer 24 μm imaging, in order to select a sample of star-forming galaxies at 1.4 < z < 2.5. This small sample of BzK galaxies has the advantage of having, in addition to the availability of a spectroscopic redshift, also a very rich photometric data set including full far-UV to NIR coverage, Spitzer (IRAC, MIPS + 16 μm InfraRed Spectrograph peak-up image; Teplitz et al. 2011; Pannella et al. 2014) and Herschel (from both PACS and SPIRE) observations and high S/N radio detections at 1.4 GHz (Morrison et al. 2010). All the BzKs have redshift in the range 1.4 ≲ z ≲ 1.6, that is in the lowest redshift regime probed by BzK selection. This is due to the requirement, by D10, of a radio detection coupled with the choice of observing the CO[2–1] transition. Its rest-frame frequency of 230.538 GHz can be observed with the Plateau de Bure Interferometer (PdBI) only up to z = 1.87.

Among the six targeted galaxies, five redshifts were obtained through the GOODS-N campaigns at Keck using DEIMOS (Stern et al. 2010). The redshift for BzK-12591 was instead derived by Cowie et al. (2004). This galaxy, showing a strong bulge in the Hubble Space Telescope (HST) imaging, has also a possible detection of [Ne v] λ 3426 Å emission line, suggesting the presence of an AGN.

For these objects, thanks to the availability of data from MIR to mm wavelengths, Magdis et al. (2012, M12 hereafter) have recently provided robust estimates of their dust masses (log M_dust in the range 8.52–9.11), based on the realistic models by Draine & Li (2007), independent estimates of the CO-to-H₂ conversion factor (α_CO), molecular gas masses and star formation efficiencies (SFEs) by exploiting the correlation of gas-to-dust mass with metallicity (M_gas/M_dust–Z). The SFRs have been computed using the information coming from the full coverage from MIR to sub-mm offered by Herschel data and therefore do not suffer from the same bias of SFR estimates from 24 μm only. It has been shown, in fact, that the L_IR from 24 μm can be up to a factor of ∼2 higher than the ‘real’ L_IR from Herschel (see e.g. Oliver et al. 2012).

The rich suite of empirical estimates of the main physical parameters of these galaxies, provided by both D10 and M12, allows us to put strong constraints on our best-fitting solutions. In particular, in this paper, we will focus only on the results concerning the radio emission, the SFR and the SFH of these objects. A full physical characterization of their molecular gas and dust properties including a detailed comparison of our predictions to empirical estimates based on observations is matter of a forthcoming paper (Lo Faro et al. 2014).

3 SED MODELLING WITH grasil

The approach used here to physically characterize these galaxy populations and in particular their radio emission, is based on galaxy evolution synthesis technique.

When modelling the SEDs of star-forming galaxies, dust effects become crucial, particularly at high redshift, so we need to include an appropriate dust model accounting for both the absorption and thermal re-emission from dust. Several works dealt with the radiative transfer (RT) in spherical geometries, mainly aimed at modelling starburst galaxies (e.g. Rowan-Robinson 1980, 2012; Rowan-Robinson & Crawford 1989; Efstathiou, Rowan-Robinson & Siebenmorgen 2000; Popescu et al. 2000; Efstathiou & Rowan-Robinson 2003; Takagi, Arimoto & Hanami 2003a; Takagi, Vansevicius & Arimoto 2003b; Siebenmorgen & Krügel 2007). Early models of this kind did not include the evolution of stellar populations. Silva et al. (1998) were the first to couple RT through a dusty ISM and the spectral (and chemical) evolution of stellar populations.

To model the emission from stars and dust consistently in order to get reliable estimates of the main physical parameters of galaxies (stellar mass, average extinction, SFR etc.), we need to solve the RT equation for idealized but realistic geometrical distributions for stars and dust as well as taking advantage of a full multiwavelength coverage from far-UV to radio. The grasil spectrophometric code (Silva 1999; Silva et al. 1998, 2011) satisfies all these requirements. For a detailed description of the code, we defer the reader to the original works. Here we provide a brief summary of its main features relevant to our work.

3.1 grasil main features

grasil is a self-consistent physical model able to predict the SEDs of galaxies from far-UV to radio including a state-of-the-art treatment of dust extinction and reprocessing based on a full RT solution. It computes the RT effects for three different dusty environments: (i) dust in interstellar H i clouds heated by the general interstellar radiation field of the galaxy (the ‘cirrus’ component), (ii) dust associated with star-forming MCs and H ii regions (dense component) and (iii) circumstellar dust shells produced by the windy final stages of stellar evolution.

It accounts for a realistic geometry where stars and dust are distributed in a bulge and/or a disc profiles. In the case of spheroidal systems, a spherical symmetric distribution with a King profile is adopted for both stars and dust, with core radius r_c.

Disc-like systems are modelled by using a double exponential of the distance from the polar axis, with scale radius R_d, and from the equatorial plane with scaleheight z_d.

In the following, we assign the same scalelengths to the stellar and dust distributions.

The clumping of young stars and dust within the diffuse medium together with the accounting for a realistic geometrical distribution for stars and dust, gives rise to an age-dependent dust attenuation which is one of the most important feature of this approach as widely discussed in BLF13 (see also Granato et al. 2000; Panuzzo et al. 2007). The dust model consists of grains in thermal equilibrium with the radiation field, and small grains and PAH molecules fluctuating in temperature.

3.2 Input star formation histories

The input SFHs are computed with che_evo (Silva 1999), a standard chemical evolution code which provides the evolution of the SFR, M_gas and metallicity, assuming an IMF, a SF law SFR(t) = ν_Sch · M_gas(t)^k + f(t), (i.e. a Schmidt-type SF with efficiency ν_Sch and a superimposed analytical term to represent transient bursts), and an exponential infall of gas (dM_inf/dt ∝ exp (−t/τ_inf)).

By varying the two parameters, ν_Sch and τ_inf we are able to recover a wide range of different SFHs, from smooth ones for large values of τ_inf, to ‘monolithic-like’ ones characterized by very short infall time-scales. A very short τ_inf can be used to have the so-called close box chemical evolution model, which ensures that the gas going to form the galaxy is all available at the beginning. Fig. 2 shows some examples of the possible SFHs which can be implemented in the chemical evolution code.

Figure 2.

SFHs from chemical evolution code: from top left to bottom right we show, as examples several SFHs obtained with our chemical evolution code. Top: exponentially declining SFHs at fixed ν_Sch for increasing values of τ_inf (left), then at fixed τ_inf and decreasing values of ν_Sch (right). When we consider longer infall time-scales and lower efficiencies our SFHs assume the form of a classical delay τ model (see the text for details). When, instead a shorter infall time-scale is considered, a more rapidly declining SFH, resembling the so-called τ models, is obtained. The bottom-left panel shows a constant SFR while the bottom-right panel shows a typical ‘starburst’ SFH. The solid black vertical lines highlight an assumed observing time of 2 Gyr. In the bottom-right panel, the time at which the galaxy is observed is highlighted in the figure.

Open in new tabDownload slide

In the figure, it can be seen that by coupling the gas accretion phase with the depletion due to SF we get SFHs which closely resemble, in their functional form, the so-called delay τ models of Lee et al. (2010) (see also the seminal work by Sandage 1986). These are in fact usually characterized by an early phase of rising SFRs with late-time decay [SFR ∝ t^βe^{− t/τ} according to Lee et al. (2010) formalism] and are probably the most suitable to explain the SEDs of high-redshift galaxies and also some local galaxies (Gavazzi et al. 2002).

The average SFHs of the 31 (U)LIRGs whose radio properties are analysed here were first presented in our previous work by BLF13 and are now discussed in more detail in Section 4.1.1 and Fig. 3 where the best-fitting SFH of each object is shown. The SFHs of the six SF BzK galaxies at z ∼ 1.5 are instead discussed in Section 4.2.

Figure 3.

Best-fitting SFHs of z ∼ 1–2 (U)LIRGs. The long dashed red vertical lines indicate the age at which the galaxy is observed. In each panel, the two parameters ruling the SFH are also specified, namely the infall time-scale and SFE. The two objects labelled in red are those for which the inclusion of radio data substantially changed the SFH (see Fig. 6 for details).

Open in new tabDownload slide

In the following, we adopt k = 1, f(t) exponential, and a Salpeter IMF which is the default choice for the chemical evolution code. When required, we convert to the Chabrier IMF by dividing by 1.7.

Our reference library of Simple Stellar Populations (SSPs) is from Bressan, Granato & Silva (1998) and Bressan et al. (2002), which directly includes the effects of dusty envelopes around asymptotic giant branch (AGB) stars and the radio emission, as described below.

3.3 Radio emission in grasil

In the SSP models by Bressan et al. (2002, B02 hereafter) and Vega et al. (2008), which we adopt here, the radio emission is assumed to be the sum of two different contributions: thermal emission from free electrons in H ii regions, and synchrotron emission from relativistic electrons accelerated by CCSN explosions (Condon & Yin 1990; Condon 1992).

The thermal component is assumed to be proportional to the number of H ionizing photons directly derived from the stellar populations and thus dependent on their age, metallicity and IMF, and scales as ≈ν^−0.1 (see equation 1 in B02).

The NT radio emission is assumed to be proportional to the CCSN rate (see section 3 and equations 3, 5 and 17 in B02 for the demonstration of this assumption). It is computed by calibrating the total NT emission/CCSN (Type II SNe, SNeII) rate relation (proportional to the average luminosity per supernova, SN, event) on the observed properties of our Galaxy, thus using the Milky Way (MW) SNII rate and total synchrotron luminosity, and accounting for the ≈6 per cent contribution by SN remnants (equations 5 and 17 in B02). The CCSNe rate is directly provided by the SSPs, since it represents the death rate of stars more massive than 8 M_⊙, i.e. ν_CCSNe = Φ(m_d) · dm_d/dt, with m_d(t) being the initial mass of the dying star in the stellar population of given age t.

As the calibration of B02 depends on the time sampling of the SSPs via dm_d/dt, this was later increased by 30 per cent due to a finer re-sampling of the time steps, in order to maintain the same radio emission under the same conditions (O. Vega, private communication).

The model has been shown to reproduce very well the FIR-to-radio correlation of normal star-forming galaxies, namely q_1.49 GHz = 2.3, as well as the radio emission of local starburst and (U)LIRGs.

We emphasize that, in our models, the SFR depends linearly on the gas fraction (Schmidt law), while the FIR emission depends on the SF, gas fraction and metallicity history, and the SN rate on the recent SFR. The consistency between FIR emission, radio emission and SN rate is thus remarkable and should be considered as a successful test of the model (B02).

4 RESULTS

The models used to interpret the observed SEDs are selected from different libraries self-consistently generated with grasil. Each library includes thousands of spectra corresponding to different combinations of model parameters. The total number of models available is ∼10⁶. Depending on the geometry (spheroid or disc) and on the assumed SFH (with or without burst), the number of free parameters, (described in detail in BLF13), typically ranges between 6 and 9. The number of photometric data points is always larger than 18.

When spectral information is also available, as for example the IRS spectra in our case, the best-fitting solution to the photometry only is then compared to the spectrum and if necessary the specific parameters influencing that region of the spectrum are fine-tuned in order to reproduce both the photometry and the spectrum. This is done for all the solutions having a reliable χ² value.

The best-fitting model provides the physical parameters of the galaxies, in particular their gas, dust and stellar mass, instantaneous SFR and SFH, optical depth and attenuation.

4.1 Interpretation of z ∼ 1 LIRGs and z ∼ 2 (U)LIRGs

As clearly stated in Section 2.1.1 our major interest here is determining whether the best-fitting solutions obtained for the overlapping galaxy sample, (z ∼ 1–2 (U)LIRGs), in our previous paper (BLF13) are able to reproduce, within the uncertainties, the observed radio flux densities. In other words, we want to check if the amount of massive young stellar populations dominating the most recent SFR is well constrained by our model as well as the fraction of stars outside the MCs contributing to cirrus heating. Given the dominant role of cirrus emission to the L_IR of BLF13 (U)LIRG SEDs, it is important to understand if the predictions for the FIR are consistent with those from the radio emission and if the conspicuous cirrus component contributing to the FIR is an ‘effect’ of a poor parameter exploration or if it is real, and required in order to reproduce the NIR-to-FIR properties of our galaxies.

So given the best-fitting SEDs computed in BLF13 using the procedure described in Section 4 and corresponding to the set of model parameters discussed in BLF13 we have simply added to the photometric data points the radio flux densities and extended the fit to the radio.

Below we give a very brief summary of the main results obtained in BLF13.

According to BLF13, all these galaxies appear to require massive populations of old (>1 Gyr) stars and, at the same time, to host a moderate ongoing activity of SF with typical SFRs ≲ 100 M_⊙ yr⁻¹. The bulk of the stars appear to have been formed a few Gyr before the observation in essentially all cases (see also Fig. 3). Their average extinctions and stellar masses are found to be higher with respect to estimates based on optical-only SED-fitting procedures. In particular, the stellar mass difference is found to be larger for the most dust obscured (U)LIRGs at z ∼ 2 for which it reaches a factor of ∼6 [with the median Δ log (M_⋆) ∼ 0.16 and 0.36 dex for the z ∼ 1 and z ∼ 2 (U)LIRGs, respectively]. The predicted SFRs are found to be lower than those computed from L_IR using the Kennicutt relation due to the significant contribution to the dust heating by intermediate-age (t_obs ≳ 10-90 Myr) stellar populations through ‘cirrus’ emission (∼73 and ∼66 per cent of total L_IR for z ∼ 1 and 2 (U)LIRGs, respectively).

4.1.1 Best-fitting far-UV-to-radio SEDs and SFHs

Fig. 4 shows the grasil best fits to the far-UV-to-radio SEDs of our z ∼ 1 LIRGs and z ∼ 2 (U)LIRGs listed according to their ID. The 3σ upper limits are indicated as red arrows while the insets report the fit to their IRS MIR spectra (discussed in detail in BLF13). The fit obtained in our previous work reproduces, well within a factor of 2, the far-UV to radio emission for almost all the (U)LIRGs into our sample without re-fitting. 28/31 (U)LIRGs show modelled radio fluxes within data error bars.

Figure 4.

grasil best fits (solid black line) to the observed SED (red circles) of z ∼ 1–2 (U)LIRGs. IRS spectra appear in the inset window (magenta line). The colour-coded lines represent the unextinguished starlight (orange dot–dashed), extinguished starlight (green dashed), cirrus emission (magenta dotted), MC emission (blue long dashed) and thermal and NT radio components (salvia long line-dashed and brown dot-long dashed, respectively). The latter are shown in the wavelength range where their contribution is significant. We are able to reproduce, well within a factor of 2, the far-UV to radio emission for almost all the (U)LIRGs into our sample. 28/31 (U)LIRGs show modelled radio fluxes within data error bars.

Open in new tabDownload slide

Figure 4

– continued

Open in new tabDownload slide

Figure 4

– continued

Open in new tabDownload slide

Figure 4

– continued

Open in new tabDownload slide

The inclusion of radio data into our models thus seems to confirm our physical solutions. In particular, as discussed in BLF13, given the detailed shape of the broad-band SED our physical analysis appears to be able to give important hints on the main parameters ruling the source's past SFH, i.e. τ_inf and ν_Sch as shown in Fig. 3.

In BLF13, we investigated both SF models with and without a starburst on top of the Schmidt-type part of the SF law (see Fig. 2). For the majority of our (U)LIRGs, a suitable calibration of the τ_inf and ν_Sch allowed us to obtain good fits to the observed SEDs with the continuous models. Fig. 3 shows the best-fitting SFH obtained for each single object in our sample and the relative τ_inf and ν_Sch. The vertical red line highlights the time, in Gyr, at which the galaxy is observed t_gal.

Small values for τ_inf (in the range 0.01–0.1 Gyr) and high values for ν_Sch (in the range 0.8–1.4 Gyr⁻¹), corresponding to an early fast and efficient SF phase, are required for 16/31 objects, 4 LIRGs and 12 z ∼ 2 (U)LIRGs. These are also the objects showing the strongest stellar bump in the rest-frame NIR (see for example L5134 and U4367 in Fig. 4). Smoother SFHs characterized by longer τ_inf (ranging between 0.4 and 1.0 Gyr) are instead required for the six remaining z ∼ 1 LIRGs and 9/21 z ∼ 2 (U)LIRGs. These galaxies present almost ‘flat’ rest-frame NIR bands and higher UV fluxes.

As already discussed in BLF13 and shown here in Fig. 3 all our (U)LIRGs appear to include massive populations of old (>1 Gyr) stars with the bulk of stars formed a few Gyr before the observation in essentially all cases. This seems to correspond to a formation redshift of z ∼ 5–6. Average estimates can be inferred from Fig 3: for the 12 z ∼ 2 (U)LIRGs characterized by very peaked SFHs, (τ_inf ∼ 0.01), about 66–80 per cent of the stellar mass has formed within ∼1 Gyr from the beginning of their SF activity, proportionally higher at higher SFE. For the z ∼ 2 (U)LIRGs presenting, instead, more regular SFH about 30–43 per cent of stellar mass is formed within 1 Gyr. Finally, for the LIRGs, on average, ∼28 per cent of the stellar mass has already been formed within the first 1 Gyr.

4.1.2 Constraints on physical solutions

We recall that all the best fits shown in Fig. 4 have been obtained by assuming the spherically symmetric King profile.

F10 performed a rough analysis of their morphology, based on HST Advanced Camera for Survey (ACS) images in the B, V, i, z filters. Many of these have been found to be extended sources characterized by many clumps. Some of them, instead, have been found to be very compact objects. All of them show red colours indicative of significant dust contribution. However, given their high redshifts it is difficult to draw a complete picture of their morphology. We believe our approximation of spheroidal geometry to be a good choice for these objects, as further confirmed by our results.

We have anyway tested on this sample of high-z (U)LIRGs also a comprehensive library of disc galaxies (∼400 000 models described in Section 4.2), and we have found that for almost all of them a spheroidal geometry is the best choice (see discussion below).

Figure 5.

Logarithmic difference between the rest-frame |$L^{{\rm GR}}_{1.4\,\mathrm{GHz}}$| luminosity computed by grasil and that one estimated from the observed flux density at 1.4 GHz using the relation of equation (2), as a function of the rest-frame observed L_1.4 GHz. 3σ upper limits are highlighted as the open red triangles. Errors on luminosities are not reported here in order to avoid crowding in the figure. They are anyway shown in the best fits. A general trend appears in our scattered solutions with the Δ L_1.4 GHz decreasing at increasing L_1.4 GHz, meaning that at larger rest-frame L_1.4 GHz our model tends, on average, to underpredict the observed radio emission. This indeed pertains only to the ≲ 2σ detections and no trend as a function of luminosity is instead evident when looking to the >4σ detections (filled symbols). We stress here that the values shown in this plot, for the high-z (U)LIRGs, refer to the best-fitting solutions obtained in our Paper I (Lo Faro et al. 2013) and as we will see in Section 4.1.3 the new solutions considered for the objects U4812 and U5152 bring these two objects on the Δ L_1.4 GHz = 0 line.

Open in new tabDownload slide

The first case, U5152, is a ‘genuine’ critical case as it is a 4σ detection in the catalogue. For this object, we measure a |$\Delta L_{1.4\,\mathrm{GHz}} = (L^{{\rm GR}}_{1.4\,\mathrm{GHz}}-L_{1.4\,\mathrm{GHz}})$| of about 0.49 dex. Moreover, for this object we also fail in well reproducing the FIR peak of the spectrum. The FIR modelled emission appears to be much broader and hotter than that suggested by the two IR data points, in particular the PACS 160 μm flux. The UV-to-MIR part of the spectrum is instead well reproduced. We have thus investigated for this object new solutions including a different geometry and also the possibility for a ‘late’ burst of SF on top of its SFH in order to boost the radio contribution from young stars leaving almost unchanged the rest of the SED. These new solutions are discussed in detail in Section 4.1.3.

The three objects U5801, U16526 and U5059 are instead faint detections at ∼2σ so we do not consider these objects as failures of our SED-fitting solutions. Moreover, their overall SED from far-UV to radio is very well reproduced by our model, with radio fluxes well within the error bars. For the object U5801, which shows the largest discrepancy on the radio flux, the FIR peak is defined by only one data point. This can bring to some degeneracy in the solutions being the contribution by young stars less constrained.

A general trend in our scattered solutions appears in Fig. 5 with the Δ L_1.4 GHz decreasing at increasing L_1.4 GHz, meaning that at larger rest-frame L_1.4 GHz our model tends to underpredict the observed radio emission. This holds, however, only for the ≲ 2σ detections. If we consider, in fact, the high S/N radio detections (filled symbols) together with the BzK galaxies not AGN dominated, there is no evidence for a variation of the Δ L_1.4 GHz as a function of luminosity. We emphasize here that all the solutions shown in the plot for the high-z (U)LIRGs refer to the best fits obtained in our previous paper BLF13 and that, considering the new solutions discussed in Section 4.1.3, the two ULIRGs U4812 and U5152 are shifted to the Δ L_1.4 GHz = 0 line. For the low S/N radio detections shown here, the average tendency of our models to underpredict the observed data could be indicative of the need of a little burst of SF on top of our best-fitting SFHs, (as discussed above), able to increase the contribution by young stars to the radio emission leaving the rest of the SED unchanged. Another possibility is to consider a different IMF characterized by a larger fraction of massive young stars, with respect to the Salpeter, as the Chabrier for example. Anyway any strong conclusion in this context is prevented by the high uncertainties in the observed radio flux densities of most of our radio faint sources.

4.1.3 A test-case: U4812 and U5152

U4812, together with U5652 and U5152, are the only three objects at z ∼ 2 to have been detected at ≳ 4σ level and listed in the published radio catalogue by Miller et al. (2013). These are also among the objects showing the highest fluxes in the FIR Herschel bands. For U5652, our previous physical solution reproduces very well the radio data with a Δ L_1.4 GHz ∼ 0. For the remaining two ULIRGs, the modelled radio emission provided by our fits tends to underpredict the observed flux densities by ∼0.2 and 0.49 dex for U4812 and U5152, respectively. Although the discrepancy between the modelled and observed radio fluxes of U4812 is well within a factor of 2, given that this object is one of the few detected at high significance in the radio band we decided to consider new solutions involving a new SED fitting for both the objects.

Based on the BLF13 analysis, U4812 has a TIR luminosity of ∼5.0 × 10¹²L/L_⊙ corresponding to a best-fitting model SFR₁₀ ∼ 160 M_⊙ yr⁻¹. It is also the object showing the largest stellar mass discrepancy (about a factor of 4 to 5), with respect to the estimate based on optical-only SED-fitting procedure, and the largest dust obscuration with an average value of A_V ∼ 3.76. According to our previous analysis, its best-fitting SFH is characterized by a very short infall time-scale and high efficiency of SF (τ_inf = 0.01 Gyr, ν_Sch = 0.88), corresponding to an early fast and intense SF phase with an initial burst followed by a more regular SFR [see Fig. 6 (left: black solid line)]. The galaxy, however, is then observed a few Gyr after the peak. As shown in Fig. 7 (top left) the best-fitting SED based on our previous solution works very well in the far-UV-to-sub-mm range but tends to underpredict the radio flux at 1.4 GHz. This is probably due to the significant contribution (much stronger with respect to the other objects in the sample) in our old best fit by the cirrus component to the rest-frame MIR region. A good compromise is given by the new solution shown in the right-hand panel of Fig. 7. This has been obtained by running the SED-fitting procedure described in Section 4 also on an extended library of disc galaxies including more than ∼400 000 models described in Section 4.2. In the new solution, the contribution of cirrus emission to the FIR remains unchanged, but its role in the MIR part of the SED is overcome by the MC emission whose effect is that of increasing the young component contributing to the radio continuum. In this way, we are able to reproduce the entire SED of the galaxy, included the IRS spectrum, although the depth of the silicate feature and continuum appear to be better reproduced in the first case. This new solution corresponds to the best-fitting SFH shown in Fig. 6 (left) as a red dashed line. This SFH is clearly different with respect to the previous one. Here, the SFH is characterized by a longer infall time-scale, of 4 Gyr and also a higher efficiency of SF. Therefore, it resembles a more gradually evolving SFH typical of BzK galaxies at high redshifts. The best galaxy age is found to be 2.66 Gyr, very close to our previous estimate of 2.3 Gyr. In this new solution the galaxy is observed closer to the peak of SF, which contributes to enhance its SFR to 316 M_⊙ yr⁻¹. Note that this best fit has been actually obtained using the disc geometry. This geometry in combination with the SFH does not seem to affect significantly both the average extinction (⁠|$A^{{\rm disc}}_{\rm V} = 3.92$|⁠) and the FIR luminosity (L_IR = 4.60 × 10¹²L/L_⊙) of this solution with respect to the previous one in BLF13. It affects, instead, the distribution of dust among the dense and diffuse components, enhancing the contribution from MCs and resulting in a lower dust mass (⁠|$\log M^{{\rm disc}}_{\rm dust} = 9.00$| compared to our previous estimate of |$\log M^{{\rm sph}}_{\rm dust} = 9.39$|⁠), i.e. almost a factor of 2.5 difference. Concerning the stellar mass, the longer infall time-scale causes a larger stellar mass by a factor of ≲ 1.5 with respect to our previous best fit. Fig. 6 also reports the best-fitting values, in logarithmic scale, of the M_⋆ relative to the two best fits shown in Fig. 7 (top).

Figure 6.

Comparison between the best-fitting SFH obtained with our previous analysis (black solid line) without including the radio data and the new best-fitting SFH (red dashed line) resulting from the SED fitting in which we considered also the radio data. The best-fitting galaxy age (t_gal), infall time-scale (τ_inf), SFE (ν_Sch) and M_⋆ estimates are also highlighted in the figure.

Open in new tabDownload slide

Figure 7.

Best-fitting SEDs of the ULIRGs U4812 (top) and U5152 (bottom) obtained before (left) and after (right) the inclusion of radio data in the fit, the latter including the new solutions considered. For the object U4812, the two plots shown in top figure refer to a different geometry, spheroidal (left) and disc (right) and different SFHs as reported in Fig. 6. For details see the text.

Open in new tabDownload slide

The other case for which we have explored a different fit is U5152 whose best-fitting SED is shown in Fig. 7 (bottom). For this object our solution (bottom left) well reproduces the far-UV to sub-mm SED (with some uncertainty on the PACS 160 μm flux), but underestimates the observed radio data by a factor of ∼3. Based on the BLF13 analysis and similarly to U4812, the best-fitting SFH of this object (shown in Fig. 6 right) is also characterized by a very short infall time-scale, of the order of 0.01 Gyr, and high SFE (∼0.9) typical of an early intense SF phase with the galaxy being observed few Gyr after the peak (t_gal ∼ 2.35 Gyr).

In order to enhance the radio contribution from young stars leaving almost unchanged the far-UV to FIR SED, we have added a late starburst on top of the quiescent SFH of this galaxy. We slightly increased the SFE of the quiescent SFR in order to limit the amount of gas available for the burst. Anyway the gas mass involved in the burst accounts to only a small fraction, of the order of 2 per cent, of the galactic mass at that epoch. This burst has an e-folding time-scale of ∼10 Myr and takes place when the galaxy is 2.34 Gyr old, with the observational time being 2.35 Gyr (same as in BLF13). The SFH corresponding to this new solution is shown in Fig. 6 (right) as red dashed line and compared to the old prescription (black solid line). The red dashed vertical line in correspondence of t = 2.34 Gyr represents the burst.

Under these new prescriptions, we are able to reproduce fairly well the entire SED of the galaxy, including the IRS spectrum and radio data. The best-fitting SED corresponding to the new solution is shown in Fig. 7 (bottom right). In the new fit, the contribution of MCs emission in the MIR region is enhanced, mainly due to a lower optical depth of MCs as compared to the previous fit.

This new solution does not seem to affect significantly the main physical properties of this galaxy. Both the average extinction (A_V = 3.20) and FIR luminosity (L_IR = 1.27 × 10¹2) well agree with the previous estimates (A_V = 3.27 and L_IR = 1.09 × 10¹2). However, given the presence of a recent burst, its SFR averaged over the last 10 Myr increases by a factor of ∼3 with respect to the old solution (SFR∼ 54 M_⊙ yr⁻¹) up to 171 M_⊙ yr⁻¹. Differently from U4812, where the effect of a different geometry (disc versus spheroids) brings to a greater dust mass by a factor of ∼2.5, here the dust mass corresponding to the new solution (log (M_dust) = 8.34) is a factor of ∼3 lower than the BLF13 one, due to the enhancement of MCs contribution to the MIR (see Fig. 7, bottom) also dominating at sub-mm wavelengths. For the stellar mass, the new solution provides a log M_⋆ = 11.05, a factor of 1.5 lower than that one obtained with the old prescription (see Fig. 6, right) but still larger compared to the one based on optical-only SED-fitting procedure. The lower stellar mass can be explained as due to the decreased amount of evolved stars.

Our analysis thus demonstrates, at least for these two test cases, that radio data are crucial to break model degeneracies and further constrain the SFH and M_⋆ of the galaxy.

Table 2 summarizes all the best-fitting physical properties of the high-z (U)LIRGs including the updated solutions for the two cases discussed in this section. As discussed in our previous paper, when considering in addition to the many combinations of grasil parameters also the different combinations of parameters ruling the SFH, the typical uncertainties on our best-fitting L_IR, M_⋆, A_V and SFR₁₀ are, respectively, 0.13 dex, 0.2 dex, 0.3 mag, and 0.2 dex thus well within the typical uncertainties for this kind of analysis. Fig. 8 gives a more specific idea of the degeneracies in our model solutions, between the average extinction, stellar mass, SFR and IR luminosity colour coded by the value of χ², for a typical z ∼ 2 ULIRG. As we see, among the many solutions considered, acceptable best fits, within |$\chi ^{2}_{\rm min} + 30$|⁠, are clearly identified in the parameter space and not much degeneracy is apparent. The grey dots shown in the panels represent all the solutions having χ² ≲ 100. This large value has been chosen in order to give a better idea of the parameter coverage of our libraries. The degeneracy shown in this figure refers to the case of fixed geometry (King's profile here). Except for the ULIRG U4812 discussed above for which we have been able to obtain a better fit to the entire SED (the radio data here have been crucial) by considering disc geometries for all the other ULIRGs in the sample we did not get any physical solution when exploring a geometry different from a spheroidal one. We cannot therefore provide in this case a reliable measurement of model degeneracies including the assumption of different geometries. Based on the unique case investigated here we can indeed notice that the use of a disc geometry results in a more ‘classical’ and less ‘bursty’ SFH more typical of normal SF galaxies (see discussion below in Section 4.2). The typical uncertainties on the main physical quantities listed above settle, in this case, around a factor of 2.

Figure 8.

χ² values for spectral models of a typical z ∼ 2 ULIRG colour coded (χ² increasing from dark to light colours) as a function of the average extinction and the stellar mass (top), SFR (centre) and TIR luminosity (bottom). We show here all the solutions with χ² ≤ 100 (grey dots) in order to give an idea of the wide range of parameters (and solutions) explored, at fixed geometry, and those with χ² within |$\chi ^{2}_{\rm min}$|+30 (coloured points). In addition to the variation of the input grasil parameters fully discussed in BLF13, we consider here also the variation of the two parameters characterizing the SFH, namely ν_Sch and τ_inf. As evident from the figure not much degeneracy is seen among the model parameters. The typical uncertainty on our best-fitting M_⋆, A_V, SFR₁₀ and L_IR are, respectively, 0.2 dex, 0.3 mag, 0.2 dex and 0.13 dex thus well within the typical uncertainties for this kind of analysis.

Open in new tabDownload slide

4.2 Best-fitting SEDs and SFHs of z ∼ 1.5 BzK star-forming galaxies

We discuss here the results relative to the 6 z ∼ 1.5 BzK SF galaxies presented in Section 2.2.

In addition to the library of star-forming spheroids discussed above, it was particularly important to consider for these objects also the model libraries with disc geometry. The scalelengths R_d and z_d of the double exponential profile have been assumed to be equal for stars and dust in order to limit the number of free parameters of the model. Several values for these scalelengths have been considered in combination with the other model parameters widely discussed in BLF13. Concerning the SFH, we have considered both SFHs typical of normal SF galaxies, namely long infall time-scales (1–5 Gyr) and moderate to high SFEs (0.5–3.0 Gyr⁻¹), and more ‘bursty’ SFHs characterized by shorter infall time-scales. We have built in this way a library including more than 400 000 disc spectra.

For each BzK galaxy of the sample, we have run the SED-fitting procedure described above on the full model library. The best-fits obtained through the |$\chi^{2}_{\nu}$| minimization procedure and physical parameter analysis, are all consistent with a disc geometry, in agreement with D10 morphological analysis. These are reported in Fig. 9.

Figure 9.

Best fits to the observed SEDs (red filled circles) of z ∼ 1.5 BzK star-forming galaxies. The short green and cyan segments at 16 and 24 μm, respectively, represent the modelled fluxes integrated over the filter bands. These are required by the presence of several spectral features ascribed to PAHs falling in these wavelength range.

Open in new tabDownload slide

For all the six BzK galaxies in the sample, we are able to reproduce the observed SEDs from far-UV-to-sub-mm very well. Among these galaxies, 3/6 (BzK-16000, BzK-21000, BzK-17999) have also modelled radio emission in perfect agreement with radio data, 1/6 (BzK-4171) has modelled radio fluxes in agreement with radio data within a factor of ≲ 1.5 and for 2/6 (BzK-12591, BzK-25536) objects our solutions appear to underpredict the radio data by a factor larger than ∼2. All these results are quantitatively summarized in Fig. 5, where the z ∼ 1.5 BzKs are represented by the starred symbols.

The two objects showing a ‘deficit’ of the models with respect to the radio data, namely BzK-12591 and BzK-25536, also show a low value of the observed FIR-radio luminosity ratio q_TIR, both ∼2.06 (see Section 5 for details). The galaxy BzK-12591, as discussed in Section 2.2, shows a strong bulge in HST imaging and has a possible detection of [Ne v] λ 3426 Å emission line, indicative of the presence of an AGN. As our model does not include an AGN component, we tend to underpredict the radio emission for this object. The other galaxy, BzK-25536, has no explicit indication for a presence of an AGN but shows the same low value of q_TIR. Apart from a hidden AGN, it has been argued that also star-forming galaxies in dense environments or particularly star-forming phase may show such low-q values. Miller & Owen (2001) observed cluster galaxies with no sign of AGN and low-q values, which they ascribed to thermal pressure of the intracluster medium (ICM) causing a compression of the galaxy magnetic field and therefore a radio excess. On the other hand, B02 proposed for these same cases that an environment induced fast damping of the SF could give rise to an apparent radio excess since radio emission fades less rapidly than the IR.

However, given the nature of these galaxies and their selection we suggest that an optically obscured AGN more probably contributes to the radio emission of this object. A ‘radio excess’ for these two galaxies has been also independently claimed by Magdis and collaborators (private communication). The predicted and observed q_TIR values of each BzK galaxy in the sample, estimated according to the equation (3), are listed in Table 3 and discussed in Section 5.

From the detailed shape of the broad-band best-fitting SED of these six BzK galaxies, we have derived the SFHs shown in Fig. 10. 3/6 objects (top-left and top and bottom-centre panels) show relatively long infall time-scales (τ_inf ∼ 4 Gyr) typical of gradually evolving star-forming discs at high redshift. Two of them also have high SFEs (ν_Sch = 1.4 Gyr⁻¹) while one of them shows a lower SFE of 0.6 Gyr⁻¹ typical of normal SF spiral galaxies. The remaining three BzK show SFHs characterized by shorter infall time-scales (τ_inf ∼ 0.5 Gyr) and moderate SFEs (ν_Sch = 0.8 Gyr⁻¹), typical of objects in earlier phases of SF. The latter also show smaller ages with respect to the former, typically in the range between 1.78 and 2.66 Gyr. All of them are observed within 2 Gyr from the peak of the SF activity, showing a moderate ongoing activity of SF with SFRs ≲ 107 M_⊙ yr⁻¹ typical of MS star-forming galaxies (D10). Anyway the statistics of the sample is still too low to draw a self-consistent picture of their SFHs.

Figure 10.

Best-fitting SFHs of z ∼ 1.5 SF BzK galaxies. The long dashed red vertical lines indicate the age at which the galaxy is observed.

Open in new tabDownload slide

Table 3 summarizes the main results of our analysis and compare them to the estimates provided by D10 and M12. Differently from D10, whose SFR estimates rely only on 24 μm, the SFRs from M12 have been computed using the information coming from the full coverage from MIR to sub-mm offered by Herschel. As already stated above, the L_IR from 24 μm can be up to a factor of ∼2 higher than the ‘real’ L_IR measured with Herschel (see e.g. Oliver et al. 2012). As shown in Table 3 our SFRs, averaged over the last 10 Myr, well agree with those derived by M12. For comparison, we also report the SFRs based on the Kennicutt calibration. These differ with respect to our estimates by a factor lower than ∼1.4, on average, well within the typical uncertainties for this kind of measure.

The stellar masses derived from the best-fitting SED (SFH) of these objects are also listed in Table 3 and compared to the estimates provided by D10.

Stellar masses in D10 (column 9 of Table 3) are derived by fitting the Maraston (2005, M05 hereafter) SSP models to the UV–optical–NIR (up to 5.8 μm) band of each galaxy. M05 models include a particularly strong contribution of AGB emission, much stronger than e.g. B02 and Bruzual & Charlot (2003, BC03 hereafter), both the latter are in fact based on the same Padova stellar isochrones.

Given that all the galaxies in the sample are star forming, a constant SFR in combination with a large range of ages, a Chabrier (2003) IMF and different metallicities from half solar to twice solar, is adopted by D10. Dust effects are accounted for by assuming a homogeneous foreground screen of dust and the Calzetti et al. (2000) reddening law. The 1σ error on D10 stellar masses is of the order of 0.10–0.15 dex. As evident from Table 3, D10,M_⋆ are much lower than our estimates by a factor ranging between 2.4 and ∼7 (for BzK-17999). As shown in M05 and Maraston et al. (2006) the differences in the (1) stellar evolutionary models used to construct the isochrones, (2) the treatment of the Thermally Pulsing (TP)-AGB phase and (3) the specific procedure used for computing the integrated spectra can lead to large differences in terms of stellar ages and masses when comparing the BC03 stellar models to the M05 ones. In particular, M05 SSPs are typically brighter and redder than BC03 ones, this results in lower ages and lower stellar masses with respect to BC03 by a factor of ∼60 per cent. We believe that, once corrected for the use of different stellar population models, the combination of a different SFH and dust attenuation treatment is again the major source of discrepancy as demonstrated by BLF13.

5 FIR–RADIO CORRELATION: COMPARING MODELS TO OBSERVATIONS

Despite the FIR/radio correlation is now well established up to high redshifts (e.g. Ivison et al. 2010; Sargent et al. 2010; Mao et al. 2011; Pannella et al. 2014), its physical origin is still debated.

B02 found the tightness of the FIR/radio correlation to be natural when the synchrotron mechanism dominates over the inverse Compton, and the electron cooling time is shorter than the fading time of the SN rate. Both these conditions are met in star-forming galaxies, from normal spirals to obscured starbursts. However, since the radio NT emission is delayed, deviations in the correlation are expected both in the early phases of a starburst, when the radio thermal component dominates, and in the post-starburst phase, when the bulk of the NT component originates from less massive stars.

Fig. 11 shows the distribution of the predicted q_TIR values for our z ∼ 1 LIRGs (blue horizontal lines) z ∼ 1.5 BzKs (green wide spaced horizontal lines) and z ∼ 2 (U)LIRGs (black 45° angled solid lines). The two distributions of LIRGs and (U)LIRGs are quite similar with the lower redshift LIRGs reaching higher values of q_TIR. The q_TIR distribution of BzK is shifted instead towards lower values. Mean values and standard deviations are also highlighted in the figure and compared to the observational estimates provided by Sargent et al. (2010; in red).

Figure 11.

Distributions of the predicted q_TIR values for the z ∼ 1 LIRGs (blue horizontal lines) z ∼ 1.5 BzK (green wide spaced horizontal lines) and z ∼ 2 (U)LIRGs (black 45° angled solid lines). Mean values and standard deviations are also highlighted in the figure. The red arrow and horizontal line indicate the typical q_TIR of z ∼ 1–2 (U)LIRGs, as computed by Sargent et al. (2010), and its range of values.

Open in new tabDownload slide

Sargent et al. (2010) studied the evolution of the IR/radio relation out to z ∼ 2 for a statistically significant volume-limited sample of IR-luminous galaxies selected in the COSMOS field. Their sample includes ∼1692 star-forming ULIRGs and ∼3000 SF ‘IR bright’ (L_TIR ≥ L^knee_TIR(z)) COSMOS sources up to z ∼ 2. They found no evolution of the median TIR/radio ratio among the ULIRG sample, with a median value at z ∼ 2 of |$2.672^{+\infty }_{-0.121} (2.892)$|⁠. The value within brackets represents the median in high-redshift bins before the correction (Δ q_TIR ∼ 0.22), needed in order to compensate for the relative offset between medians at high and low redshift that arises artificially due to the increased scatter (σ_TIR) in the data at z ≳ 1.4 (Sargent et al. 2010).

Our predicted q_TIR for z ∼ 2 (U)LIRGs (2.76 ± 0.10), well agrees, within the errors, with the results from Sargent et al. (2010) discussed above and also with the median q_TIR measured by the authors for local (U)LIRGs (⁠|$\left\langle \,q_{\mathrm{TIR}}\,\right\rangle \,2.703_{-0.050}^{+0.050}$|⁠). Our results thus confirm, from a modellistic point of view and in line with theoretical and numerical simulation expectations, that ULIRGs should follow the local IR–radio relation until at least z ∼ 2 (Murphy et al. 2009; Lacki & Thompson 2010). This implies that magnetic fields are sufficiently strong to ensure cosmic rays electrons to predominantly lose energy through synchrotron radiation rather than inverse Compton scattering off the CMB.

Good agreement, within the errors, is found also between our predicted q_TIR for z ∼ 1 LIRGs and the median q_TIR derived by Sargent et al. (2010) for z ∼ 1 SF IR-bright sources (⁠|$2.672^{+0.069}_{-0.061}$|⁠). These are actually very similar to the median q_TIR measured for z ∼ 2 SF ULIRGs. In Fig. 11 this is emphasized by the two samples showing very similar q_TIR distributions. An even better agreement with the observational estimates provided by Sargent is found for the six BzKs which show an average q_TIR of 2.65±0.05. Of course given the low statistics of this sample we cannot draw strong conclusion from this comparison. It gives anyway important hints about our solutions. No strong evolution of the FIR/radio correlation is thus observed also between our z ∼1 and ∼2 L ≳ 10¹¹ L_⊙ objects.

Our physical model thus seems to be able to reproduce the radio properties of high-z (U)LIRGs including the FIR–radio correlation up to z ∼ 1–2. This provides a further important constraint to model the SFR and SFHs underlying the observed SED.

6 RADIO CONSTRAINTS ON THE CURRENT SFR OF GALAXIES

Recent works by Daddi et al. (2007) and D10 have found very good agreement, within a factor of ∼2, between different SFR indicators (UV dust-corrected, MIR and radio 1.4 GHz) for a GOODS sample of BzK-selected galaxies including both individual sources and stacked sources. All these SFR estimates, however, rely on calibrations based on similar assumptions, namely K98. This calibration assumes that the L_bol of a constant SF lasting 100 Myr is totally emitted in the IR (Leitherer & Heckman 1995, LH95 hereafter; K98). For a constant SF, the L_bol after the first 10 Myr evolves relatively slowly because the rate of birth and death of the most massive stars (with lifetimes ≲ 10 Myr and dominating the L_bol) reaches a steady state (see figs 2 and 8 of LH95). The K98 SFR/L_IR calibration adopts the mean bolometric luminosity for a 10–100 Myr continuous SF, solar abundance, Salpeter IMF of the starburst synthesis models of LH95, and assumes that L_IR = L_bol.

In BLF13, we have shown that, due to the significant contribution of cirrus emission to the total L_IR whose heating source includes already evolved stellar populations (ages older than the typical time-scale for young stars to escape from the parent MCs, (t_esc) typically ranging between ∼3 and ∼90 Myr), our inferred SFR₁₀ are systematically lower than those based on the K98 calibration, by a factor of ∼2 to 2.5.

The inclusion of radio emission into our procedure is therefore extremely important to understand if this conspicuous cirrus component contributing to the FIR is an ‘effect’ of a poor parameter exploration or if it is real, and required in order to reproduce the NIR-to-FIR properties of our galaxies. One possibility in fact is that, forcing the contribution by intermediate-age ( ≳ 50 Myr) stellar populations to reproduce the NIR peaked fluxes of our galaxies, we get the right FIR emission, as the integrated luminosity (mass) of these stars is large, but we lack in SNII production. This increases the L_FIR/L_1.4 GHz ratio of the model but fitting the FIR we underpredict the radio. It is important, therefore, to check any possible systematics within the model.

Fig. 12 compares our model SFRs averaged over the last 10 Myr, thus sampling the most recent SF activity, to the SFR obtained from the rest-frame radio luminosities at 1.4 GHz, as provided by the best-fitting physical model, using the Bell et al. (2003) calibration of equation (4). For the comparison, we have re-scaled this relation, originally computed for a Salpeter IMF to a Chabrier one. Red filled circles and blue squares represent, respectively, our z ∼ 2 and ∼1 (U)LIRGs while the starred symbols mark the z ∼ 1.5 BzKs. The green solid line represents the SFR estimates based on the K98 calibration plotted as a function of the SFR_1.4 GHz^Bell derived from the modelled rest-frame radio luminosity. Its slope is ∼0.90. The dashed blue line is the 1–1 correlation line.

Figure 12.

Comparison of our grasil-estimated SFR₁₀, averaged over the last 10 Myr, with the SFR derived from the modelled (filled points) rest-frame L_1.4 GHz luminosity using the Bell et al. (2003) calibration (equation 4). The z ∼1.5 BzK galaxies are shown as green open stars. The 1–1 correlation is shown as blue dashed line. The green solid line represents the SFR estimates provided by the K98 calibration plotted as a function of the radio SFRs derived from the modelled L_1.4 GHz. The slope of the green line is ∼0.90. The yellow shaded area highlights the region of the plot where the SFR estimates are within a factor of ∼2 from the 1–1 correlation line. While a very good agreement is clearly visible between our SFR₁₀ and the best-fitting radio based |${\rm SFR}^{{\rm Bell}}_{1.4\,\mathrm{GHz}}$|⁠, the SFRs based on the Kennicutt relation appear to be systematically off-set with respect to those provided by radio luminosity. In particular, we observe that the discrepancy between the FIR- and radio-based SFRs seems to be larger at lower redshifts (see the text for details).

Open in new tabDownload slide

Our SFR₁₀ (from grasil) appear to be in perfect agreement with the |${\rm SFR}^{{\rm Bell}}_{1.4\,\mathrm{GHz}}$| (from Bell et al. 2003) derived from the best-fitting |$L^{{\rm GR}}_{1.4\,\mathrm{GHz}}$|⁠. For our model predictions, we have derived the calibration factor between modelled SFR and radio rest-frame luminosity, corresponding to the ratio |${\rm SFR}_{10}^{{\rm Salp.}}$|/|$L^{{\rm GR}}_{1.4\,\mathrm{GHz}}$|⁠, and compared it to the same factor as derived by Bell et al. (2003) (i.e. 5.52 × 10⁻²²). For the filled points and starred symbols of Fig. 12 we have measured a factor of (5.53 ± 0.38) × 10⁻²².

The estimates from the Kennicutt relation are, instead, systematically off-set with respect to those provided by radio luminosity, by a factor ranging between 1.5 and ≳ 2. In particular, we observe that the discrepancy between the FIR- and radio-based SFRs is ≳ 2 up to z ∼ 1.8 then larger than 1.5 up to z ∼ 2.5 and lower than 1.5 at higher redshifts. In other words, the discrepancy seems to be larger at lower redshifts.

In Fig. 13, we compare our model SFR₁₀ (left) and those based on the Kennicutt calibration, SFR_K (right), to the radio estimates where the rest-frame radio luminosity is computed directly from the observed flux density using the relation specified in equation (2) (left: open blue squares and red circles; right: filled green circles and squares). The blue dashed line represents, in both panels, the 1–1 correlation line, while the green solid line is the linear regression for the filled green points.

Figure 13.

Left: comparison of our grasil-estimated SFR₁₀ with the SFR derived from the observed (open points) rest-frame radio luminosity, L_1.4 GHz, using the Bell et al. (2003) calibration (equation 4). The 1–1 correlation is shown as blue dashed line. In both panels, the yellow shaded area identifies the region of the plot, where the SFRs are within a factor of 2 from the 1–1 correlation line. Right: comparison between the SFR estimates provided by the K98 calibration and the SFR_1.4 GHz^Bell derived directly from the observed rest-frame radio luminosity. The green filled circles and squares represent, respectively, the z ∼ 2 and 1 (U)LIRGs while the green open stars the z ∼1.5 BzKs of D10. The green solid line is the linear regression of the green points. Its slope is ∼0.73. A larger scatter around the 1–1 correlation line is present here. It reflects the scatter of our best-fitting solutions already discussed in Fig. 5 and pertains mostly to the faint detections. Anyway most of our predictions are well within a factor of ∼2 (yellow shaded area) from the 1–1 correlation line.

Open in new tabDownload slide

As shown in Fig. 13 (left) a larger scatter around the 1–1 correlation line is evident when considering the |${\rm SFR}^{{\rm Bell}}_{1.4\,\mathrm{GHz}}$| estimated from the observed rest-frame radio luminosity. This scatter reflects the scatter of our best-fitting solutions already discussed in Fig. 5 and pertains mostly to the z ∼ 2 2σ– ≲ 3σ detections. It is, in fact, lower for the z ∼ 1 LIRGs. As discussed above, only one ULIRG, U5152, and the two BzK dominated by an AGN show a significant (larger than a factor of 2) discrepancy between our model and the observed radio fluxes and most of our predictions are well within a factor of ∼2 (yellow shaded area) from the 1–1 correlation line. Averaging out our results, we have measured a calibration factor between our model SFRs and the observed rest-frame radio luminosity of ∼(6.06 ± 3.00) × 10⁻²², with the 6/31 (19 per cent) 3σ upper limits excluded from this computation. Also in this case our solutions are in agreement, within the errors, with the empirical calibration provided by Bell et al. (2003).

The important point here is that there is no systematic effect in our solutions which tend to scatter both above and below the empirical calibration. On the contrary, compared to the SFRs estimated from the FIR using the Kennicutt calibration, we were systematically lower by a factor of ∼2. Therefore, this seems to go in the direction of confirming our physical predictions and it gives an indication that the radio emission more than the FIR is able to accurately predict the current SFR mostly contributed by young massive stars. It does not seem to suffer from contamination by intermediate-age stellar populations, as it can happen for the FIR. Of course, we have to take into account the large errors associated with the radio data, but even with this in mind our results are in good agreement.

In the right-hand panel of Fig. 13, we compare the SFR_K based on the Kennicutt calibration and FIR luminosity (filled green circles and squares), to the radio estimates derived from the observed rest-frame radio luminosity. Also here there is a large scatter of the data points mostly above the 1–1 correlation line. The slope of the green line is ∼0.73. What is evident here is that the FIR-based SFRs tend to be higher with respect to the radio ones by a factor larger than 2 up to z ∼ 1.5. Moreover, up to z ∼ 1.8 most of the points tend to lie above the 1–1 correlation line. Beyond redshift 1.8 the discrepancy between the FIR- and radio-based SFRs decreases up to a factor below 1.5. So there seems to be, in the trend shown in this figure, a redshift dependence of the SFR based on Kennicutt calibration and on the radio flux, already highlighted in Fig. 12. As we have discussed above, the FIR, differently from the radio emission, depends on the dominant population of stars heating the dust. We have shown that when cirrus emission powered mostly by intermediate-age stars dominates at FIR wavelengths the Kennicutt calibration can overestimate the current SFR of galaxies. We have thus investigated the correlation between the logarithmic difference of the FIR- and radio-based SFRs and the fractional contribution of cirrus emission to the TIR luminosity computed as L_IR^cirr/L_IR^TOT. It is worth stressing again that in our model the sources responsible for the cirrus heating are all the stars outside the MCs, i.e. all the stars with ages older than the typical escaping time-scale. The results are shown in Fig. 14 as filled green circles (z ∼ 2) and squares (z ∼ 1). There is not a clear correlation in the first (left) plot, due to the large scatter, rather a ‘trend’ showing that the discrepancy between the two different SFR estimates becomes larger for higher fractional contributions of cirrus emission to the IR luminosity. A correlation appears to be more evident in the second (right) plot where the scatter is strongly reduced as all the quantities are derived from the model. There seems to be also a redshift dependence in the sense that being z ∼ 1 LIRGs (filled squares) characterized by higher cirrus fractions, on average, they also show the larger discrepancies.

Figure 14.

Logarithmic difference of the SFR estimated from the rest-frame IR luminosity using the K98 calibration, SFR_K, and that one derived from the rest-fame radio luminosity using the Bell et al. (2003) empirical relation, SFR_1.4 GHz^Bell, as a function of the fractional contribution of cirrus emission to the TIR luminosity as provided by the best-fitting model. In the left-hand panel, L_1.4 GHz is directly derived from the observed data using the relation 4, while in the right-hand panel, it comes from the best-fitting model. There is not a clear correlation in the left-hand panel, due to the large scatter, rather a ‘trend’ showing that the discrepancy between the two different SFR estimates becomes larger, on average, for higher fractional contributions of cirrus emission to the IR luminosity. The correlation becomes, instead, clearer in the right-hand panel, where the scatter is strongly reduced as all the quantities are derived from the model. There seems to be also a redshift dependence in the sense that being z ∼ 1 LIRGs (filled squares) characterized by higher cirrus fractions, on average, they also show the larger discrepancies.

Open in new tabDownload slide

The extension of our analysis also to low-redshift MS galaxies would be necessary in order to draw a more complete picture of the variation of the two different SFR estimators as a function also of redshift. This will be part of our future works.

7 HOW DOES OUR PHYSICAL ANALYSIS AFFECT THE SFR-M_⋆ DIAGRAM?

Many recent studies have found evidence that the SFR in galaxies correlates with their stellar mass along a MS relation which evolves with redshift and representing a ‘steady’ mode of SF (e.g. Guzman et al. 1997; Brinchmann & Ellis 2000; Bauer et al. 2005; Bell et al. 2005; Papovich et al. 2006; Reddy et al. 2006; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007; Pannella et al. 2009; Rodighiero et al. 2010, 2011; Karim et al. 2011). Above the MS with higher specific Star Formation Rates (sSFRs) (by a factor ranging between 4 and 10), there are the so-called outlier galaxies characterized by a ‘starburst’ mode of SF generally interpreted as driven by mergers. These off-MS (or outliers) galaxies have been found to contribute only 10 per cent of the cosmic SFR density at z ∼ 2 (Rodighiero et al. 2011). This has been interpreted as a further indication that most of the stellar mass forms in continuous mode of SF. Under this picture, high-z LIRGs and (U)LIRGs seem to mostly reflect the high SFR typical for massive galaxies at that epoch, so they are not brief stochastic starbursts as their local counterpart. They simply represent the early gas-rich phase of smoothly declining SFH of ≥L_⋆ galaxies, as demonstrated by our analysis.

The availability of radio data has allowed us to put stronger constraints on our solutions, particularly in terms of the SFH and current SFR of these objects. Moreover, being the SFR and M_⋆ of a galaxy related one to each other and dependent on the age distribution of stellar populations and giving the radio constraints also on the fraction of very young stars contributing to the NT radio luminosity through CCSNe explosions, the analysis performed in this work goes in the direction of strengthening our predictions also in terms of the stellar content of these galaxies. Of course the results presented here rely on the specific choice of the IMF being a Salpeter then converted to a Chabrier IMF by assuming the conservative value specified in Section 2.

We investigate here the effects of our physical analysis on the SFR–M_⋆ relation described above. While our derived galaxy SFRs have been discussed thoroughly in this paper, we refer the reader to BLF13 for the details of stellar mass determination. However, as the entire discussion about the placement of galaxies in the MS depends on how our SFR and M_⋆ have been estimated, it is worth summarizing here the key results of our previous work on M_⋆ measurements.

The stellar masses of all the galaxies in our sample have been self-consistently derived by BLF13 by coupling the chemical evolution models with the RT based spectral synthesis code grasil. The stellar mass is computed by accounting for the contribution of both still surviving stars and dead remnants. The fit is performed taking into account the information coming from the entire SED from far-UV to radio and, as the dust extinction and reprocessing are computed by solving the differential equations for RT, the energy balance is always naturally preserved in the code. As already highlighted in Section 4.1, when comparing our stellar mass estimates to those provided by more classical SED-fitting procedures we found systematically larger values by a factor up to ∼6 for the most dust obscured (U)LIRGs at z ∼ 2. We also found for the same objects higher extinctions by ΔA_V ∼ 0.81 (z ∼ 1) and 1.14 (z ∼ 2).

We investigated in BLF13 the possible origin of the larger differences in our M_⋆ estimates by considering different configurations for the factors which are supposed to be the main drivers of these discrepancies, that is to say, the different SFHs (continuous versus starburst), evolutionary models, wavelength coverage and dust extinction and reprocessing parametrization. Micha lowski et al. (2012) found in the adopted SFH, (double component versus τ model), and evolutionary models the major factor (∼2.5 in M_⋆) affecting the stellar mass estimates of high-z dust obscured sub-mm galaxies when a τ model is considered. We tested all these hypotheses (see discussion in section 5.1 of BLF13) and concluded that the different SFH and SSPs cannot account for the larger discrepancies we find among the most dust obscured ULIRGs at z ∼ 2; for these the dominant factor is the dust extinction. We found in fact, the stellar mass discrepancy to be larger for higher values of A_V concluding that the stellar mass which is missed by fitting the optical data alone is hidden in dust. With our analysis, we showed how extremely idealized approaches, not accounting for the information coming from the entire observed SED and thus based on optical-only SED-fitting procedures, may produce highly degenerate model solutions and galaxy SEDs unable to energetically balance the dust reprocessed IR emission from the galaxy (see e.g. fig. 7 in BLF13). For the most dust obscured objects this can result in severely underestimated stellar masses. Similar conclusions have been obtained recently also by Mitchell et al. (2013).

Our results, in terms of larger stellar masses and lower SFRs with respect to classical estimates, seem to have a strong impact on the SFR–M_⋆ diagram.

Fig. 15 shows the stellar mass–SFR relation for different galaxies at 1.5 < z < 2.5. Most of the points are the same ones discussed in fig. 1 of Rodighiero et al. (2011) (see left-hand panel). The main thick black solid line represents the MS relation from Rodighiero et al. (2011). The light yellow region (right-hand panel) defines the locus occupied by objects having sSFR within a factor of 2 from the MS while the two solid lines above the yellow region mark, respectively, the loci 4 and 10 times above the MS (along the SFR axis). The small grey circles are BzK galaxies from both COSMOS and GOODS catalogues. The small red and cyan circles represent, respectively, PACS-COSMOS and PACS-GOODS-S detected sources. The open blue circles and green stars shown in the left-hand panel represent our sample of z ∼ 2 (U)LIRGs and z ∼ 1.5 BzKs. For all these data points, the stellar masses have been measured by using BC03 templates in combination with the hyperz code of Bolzonella, Miralles & Pelló (2000). The SFRs are derived from the TIR luminosity using the K98 calibration.

Figure 15.

SFR–M_⋆ diagram: grasil model predictions versus empirical estimates based on observations. In both panels, small grey circles are z ∼1.5–2.0 BzK galaxies from GOODS and COSMOS survey while small red and cyan circles are, respectively, PACS-COSMOS and PACS-GOODS-S detected sources. The black solid line represents the MS relation computed by Rodighiero et al. (2011), while the yellow region defines the locus occupied by objects with sSFRs within a factor of 2 from the MS. The other two black lines highlight the x4 and x10 off-MS regions. The big open blue circles and green stars shown in the left-hand panel represent, respectively, our z ∼2 (U)LIRGs and z ∼1.5 BzKs whose stellar masses have been computed with the hyperz code (thus based on optical-only SED-fitting procedure), and SFRs have been derived by using the empirical K98 calibration. The big blue filled circles and green thick stars shown in the right-hand panel represent, instead, the z ∼2 (U)LIRGs and the six z ∼1.5 BzKs whose M_⋆ and SFRs have been self-consistently derived here through our physical analysis. The magenta stars identify, a subsample of the x10-off-MS galaxies of Rodighiero et al. (2011). The big magenta star shown in the left-hand panel highlights the region where the outliers have been selected with their SFR and M_⋆ derived by using classical approaches. For the small magenta stars shown in the right-hand panel the stellar masses and SFR come from our physical analysis.

Open in new tabDownload slide

In BLF13, we have applied our physical analysis also to a sample including some of the most extreme outliers of Rodighiero et al. (2011) (those having sSFR above a factor of 8–10 with respect to MS) in order to explore its effects also on starburst classified galaxies. The big open magenta star shown in the left-hand panel of Fig. 15 highlights the region of the diagram where the outliers have been selected with their M_⋆ and SFRs computed using classical approach.

The big blue filled circles and green thick stars shown in the right-hand panel are, respectively, the z ∼ 2 (U)LIRGs and the six BzKs analysed here. Their stellar masses and SFRs have been self-consistently computed using our physical model. The same holds for the small magenta stars representing the high-z and highly SF outliers analysed in BLF13.

Our physical analysis applied to high-z dusty star-forming (U)LIRGs and BzKs thus has the effect of shifting downwards, below the MS, these objects with the exception of the two (U)LIRGs showing lower A_V. The effect strongly depends on the level of dust obscuration of the galaxy and on the nature of SF (starburst versus quiescent SF). We have seen in BLF13 that for the most dust obscured objects the stellar mass can be significantly larger than that derived from ‘optical-only’ SED-fitting approaches. As for the SFRs, we have found that they can be systematically lower, by a factor of ∼2 to 2.5, with respect to those based on K98 calibration due to significant cirrus contribution to IR luminosity, powered by intermediate-age (t_⋆ ≳ 10–90 Myr) stellar populations. The combination of these two effects on the SFR–M_⋆ relation is clearly visible when we look at the so-called outliers (magenta stars). The starting point is represented by the big open magenta star lying above the x10-off MS black line in the left-hand panel of Fig. 15. After applying our analysis we get the magenta stars shown in the right-hand panel. We can clearly distinguish in the figure two different situations, corresponding to the two well distinct groups of magenta stars. For the most dust obscured objects whose IR SEDs are dominated by an ongoing burst of SF in the MC (fig. 5 of BLF13), the SFRs derived from the model are consistent with the K98 calibration, within the errors, while the stellar masses are still larger than classical predictions and for some objects reach a factor of ∼7. The effect here is that of shifting these most extreme cases to the region between the x4- and x10-MS solid lines, still above the MS. For the dust obscured objects whose IR SED is dominated by cirrus component, thus characterized by gradual evolving SFHs, we measure SFRs lower by a factor of ∼2 to 2.5 than Kennicutt calibration, in addition to large stellar masses. Here the effect of our analysis is that of bringing back to the MS these peculiar galaxies. Similar results, for a sample of high-z sub-mm galaxies, have been recently obtained also by Micha lowski et al. (2012).

The question we would like to address now is: ‘are these’ ‘outliers’ a real population of off-MS galaxies or due to a severe mass underestimate related to their high dust obscuration, are they ‘simply’ highly star-forming MS galaxies?’. Moreover combining the two results discussed above [relative to normal star-forming (U)LIRGs and starburst galaxies], it seems that the global effect of our analysis is that of simply changing the normalization of the relation. Of course to investigate more quantitatively this possibility we need to apply the same physical analysis to a larger statistical sample of both MS and starburst galaxies to see whether there is a shift downwards (to larger stellar masses) of the entire sequence or only a compression due to the fact that only the most dust obscured objects have stellar masses which have been underestimated. Moreover, another important aspect that should be taken into consideration is the ‘real’ effect of assuming a different IMF on the estimates of stellar mass and SFR. We have adopted here a conservative value to scale from one IMF to another but assuming a different stellar mass distribution for the stars in the chemical evolution code could affect in different way, according to their SFH, the contribution by old and young stars to the UV-to-FIR SED. And this could be different from the simple assumption of a fixed scaling factor for all the stellar populations. All these aspects will be dealt in forthcoming works.

8 SUMMARY AND CONCLUSIONS

In this paper, we have extended the thorough SED analysis performed by BLF13, by including also the modelling of radio emission for a sample of high-z (U)LIRGs and BzK-selected galaxies for which a full multiwavelength data set from far-UV to sub-mm was already available.

In BLF13, we have found that the values of M_⋆ and SFR estimated through our approach can be, in some cases, significantly different with respect to previous works based on more classical methods (e.g. hyperz-mass). In particular, we have found higher stellar masses for the most dust obscured galaxies and lower SFR by a factor of ∼2 for the best fits requiring a significant contribution by intermediate-age stars (with ages older than 10–90 Myr on average) to the NIR and the FIR (through cirrus emission).

We have investigated here the origin and relevance of the cirrus component contributing to the FIR emission of our z ∼ 1–2 (U)LIRGs by complementing their far-UV to sub-mm data with VLA radio observations at 1.4 GHz carried out in the GOODS-S field. The addition of radio luminosity to the spectral multiband fitting is crucial to further constrain our inferred SFR and history considering that the radio emission is strictly connected to the most recent SF activity.

The main points of our work are summarized below.

We have collected for the sample of z ∼ 1–2 (U)LIRGs of BLF13 both published and unpublished 1.4 GHz VLA observations in the GOODS-S field. Given the low number of (U)LIRGs with high S/N radio data, and also to extend our analysis to a different population of galaxies, we have included in our study six well-sampled NIR-selected BzK galaxies at z ∼ 1.5 by D10 in GOODS-N.

For each galaxy, we have searched the best fit to their far-UV to radio SEDs, within a grasil-generated model library, with ∼10⁶ models, covering a large SF and dust-related parameter space. The modelling of the radio emission includes both the thermal and NT components, as described in B02.

We are able to reproduce, well within a factor of 2, the far-UV to radio emission for all the (U)LIRGs with rms ≳ 3–4σ, and 4/6 BzKs in our sample. The inclusion of the radio data in the SED fitting has required a different best fit for the sources U4812 and U5152 with respect to BLF13, with a different geometrical distribution for stars and dust for the former, and the addition of a late burst of SF for the latter, to enhance the radio emission.

 The extensive spectrophotometric coverage of our high-z galaxy sample up to radio allows us to set important constraints on the SFH of individual objects. For essentially all (U)LIRG and BzK objects, we find evidence for a rather continuous SFR and a peak epoch of SF preceding that of the observation by a few Gyr. This seems to correspond to a formation redshift of z ∼ 5–6.

 We underestimate the radio data of the two BzKs, BzK-12591 and -25536. The first one has confirmed AGN (Cowie et al. 2004) and shows low value of q_TIR = 2.06. The latter does not have confirmed AGN contributing to its SED but also shows the same q value. The low q observed for this galaxy may be indicative of a hidden AGN.

We have compared our predicted SFR averaged over the last 10 Myr, SFR₁₀, to the values derived from radio luminosity using the Bell et al. (2003) calibration, and from the IR luminosity using the K98 calibration. We find very good agreement between our SFR₁₀ and the |${\rm SFR}^{{\rm Bell}}_{1.4\,\mathrm{GHz}}$| derived from the rest-frame radio luminosity. The estimated SFR_K are instead systematically off-set by a factor of 1.5 to 2, with respect to those provided by the radio luminosity. We find that the discrepancy between the FIR- and radio-based SFRs is larger at lower redshifts. We also find that the agreement between the two different indicators seems to be, on average, better for the BzK-selected galaxies, as also found by D10. We do not observe any systematic effect in our solutions (in terms of SFRs), which tend to scatter both above and below the empirical calibration. On the contrary, compared to the SFRs estimated from the IR calibration, we are systematically lower by a factor of 2.

The discrepancy between the IR-and radio-based SFRs can be explained by the fact that differently from the radio emission, the FIR depends on the dominant population of stars heating the dust. We have shown that when an important contribution to the dust heating is provided by relatively evolved stars, usually via the cirrus emission, the Kennicutt calibration can overestimate the SFR. Since z ∼ 1 LIRGs are characterized by higher cirrus fractions, they also show the larger discrepancies. The extension of our analysis also to low redshift MS galaxies is required in order to draw a more complete picture of the variation of the two different SFR estimators as a function of both redshift and cirrus content.

 This comparison of the SFRs derived from different indicators (and methods) thus seem to go in the direction of confirming the BLF13 results, and provide an indication that radio emission more than FIR is able to accurately predict the SFR. It does not seem to suffer from contamination by intermediate-age stellar populations, as it happens for the FIR where, in our solution, cirrus dominates. Of course we have to take into account the large errors associated with the radio data, but even with this in mind our results are in good agreement.

We finally show how our results, in terms of larger stellar masses (BLF13) and lower SFRs with respect to previous estimates, can significantly affect the SFR–M_⋆ diagram. The amount of the shift in mass or SFR strongly depends on the level of dust obscuration of the galaxy and on the nature of SF (starburst versus gradually evolving SF).

 At high dust obscurations we estimate larger masses, but we find two different situations for the SFR: (a) for the most dust obscured objects whose IR SEDs are dominated by an ongoing burst of SF in MCs, the SFRs derived from the model are consistent with the K98 calibration. These objects are shifted in the region between the x4- and x10-MS, still above the MS; (b) for the dust obscured objects whose SED fit does not require a starburst, our estimated SFRs are a factor of ∼2 lower than the K98 calibration. We shift these objects to or below the MS. Similar results, for a sample of high-z sub-mm galaxies, have been recently obtained also by Micha lowski et al. (2012).

 To quantitatively investigate if the net effect of our analysis is that of simply changing the normalization of the relation, or a compression due to the fact that only the most dust obscured objects have stellar masses which may have been underestimated, we will need to analyse a larger statistical sample of both MS and starburst galaxies. Another important aspect now under investigation is the ‘real’ effect of assuming a different IMF on the estimates of stellar mass and SFR. In fact, assuming a different stellar mass distribution for the stars in the chemical evolution code could affect in different way, according to their SFH, the contribution by old and young stars to the UV-to-FIR SED. And this could be different from the simple assumption of a same scaling factor for all the stellar populations.

The results discussed in this paper mostly rely on VLA observations that were not completely matched in sensitivity to our Herschel and Spitzer data, and consequently only a few (U)LIRGs had reliable radio detections. Our analysis points towards the radio flux as an essential information for interpreting star-forming galaxies at high redshifts and for recovering reliable SFH. The currently operational and much improved sensitivity of the Jansky Very Large Array (J-VLA) promises to add important data as soon as the relevant cosmological fields will be surveyed at faint enough limits. And, of course, forthcoming facilities like Meerkat, or future ones like Square Kilometre Array (SKA), will provide the full extensive radio coverage needed for all the relevant use.

This paper uses data from Herschel's photometers SPIRE and PACS. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy) and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy) and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). We acknowledge support from ASI (Herschel Science Contract I/005/07/0). We thank Denis Burgarella, Veronique Buat, Alessandro Boselli, Giulia Rodighiero and Eduardo Ibar for the useful discussions which contributed to improve the paper and Maurilio Pannella to have provided us the optical-to-MIR photometry for the six BzK galaxies analysed here. We also thank the anonymous referee for his/her helpful comments which have contributed to increase the clarity of the paper.

REFERENCES