© ESO, 2013

1. Introduction

Star formation activity is a key observable for understanding the physical processes in the build up of galaxies. The SFR in galaxies depends on the physics of star-forming regions, merger history, gas infall and outflows and on stellar and AGN feedbacks. The SFR distribution and its evolution therefore offer a crucial test for any model of galaxy evolution. There are different indicators of the ongoing star formation (e.g., Kennicutt 1998), such as ultraviolet continuum (λ ~ 912−3000 Å) produced by massive, young stars (t ~ 10⁸ yr, Martin et al. 2005); nebular recombination lines from gas ionized by the hot radiation from early-type stars (λ ≤ 912 Å, t ~ 10⁷ yr); far-infrared (FIR) emission from dust heated by UV light; and non-thermal radio emission, such as synchrotron radiation from supernova remnants (see Kennicutt 1983, 1998 for their respective calibrations with SFR).

Recently, measurements of SFRs in large samples of galaxies based on the above indicators have provided interesting insights in the dominant mode of baryon accretion for star-forming galaxies: the tight correlation, with small scatter, between stellar mass and SFR observed up to z ~ 2 (e.g., Salim et al. 2007; Noeske et al. 2007b; Daddi et al. 2007; Wuyts et al. 2012) suggests that a secular, smooth process, such as gas accretion, rather than merger-induced starbursts, may be the dominant mechanism governing star formation in galaxies. Also, the relation between specific SFR and stellar mass (e.g., Brinchmann et al. 2004; Salim et al. 2007; Noeske et al. 2007a) reveals that less massive galaxies have their onset of star formation occurring later than more massive ones (e.g., Cowie et al. 1996), and that a simple mass-dependent gas exhaustion model can reproduce the observed decline of the cosmic SFR since z ~ 1−1.5 (e.g., Lilly et al. 1996; Schiminovich et al. 2005; Le Floc’h et al. 2005). Such a scenario is also consistent with recent evidence of a higher fraction of molecular gas in massive star-forming galaxies at z ≥ 1 with respect to nearby galaxies (e.g., Erb et al. 2006; Daddi et al. 2008, 2010; Tacconi et al. 2010).

Extending such measurements to high redshift for large samples of galaxies poses several challenges: optical recombination lines are often too weak and are shifted to near-IR wavelengths, where current spectroscopic capabilities are limited. Far-IR indicators are also of limited use at high redshift, since the modest sensitivity and resolution of infrared telescopes make these observations sparse and restricted to the most massive galaxies, unless stacking techniques are used (e.g., Zheng et al. 2007; Pannella et al. 2009; Karim et al. 2011; Lee et al. 2010; Heinis et al. 2013). On the other hand, with increasing redshift, the ultraviolet continuum emission is progressively shifted to optical and near-IR bands, making it easily accessible with the largest ground-based telescopes (e.g., Burgarella et al. 2007; Daddi et al. 2007; Reddy et al. 2006; Steidel et al. 1996; Bouwens et al. 2011).

In the absence of dust, the ultraviolet luminosity of a galaxy is proportional to the mean SFR on a timescale t ~ 10⁸ yr (Donas & Deharveng 1984; Leitherer & Heckman 1995). The presence of dust in the interstellar medium of galaxies hampers such a direct estimate of the SFR from UV observations. Starlight, especially in the UV, is efficiently absorbed and scattered by dust grains, which heat up and re-emit the absorbed energy at FIR wavelengths. This means that a fraction of UV photons will not escape the galaxy, and that neglecting this effect will lead to severe underestimates of the “true” SFR (e.g., Salim et al. 2007). Calzetti (1997) and Meurer et al. (1999) proposed a way to correct the UV continuum for dust attenuation that relies on the existence of a tight correlation between the slope of the UV continuum in the region 1300 ≤ λ/ Å ≤ 2600 (i.e., β-slope) and the ratio between the total IR (8 ≤ λ/  μm ≤ 1000) and ultraviolet luminosities (i.e., infrared excess, IRX ~ L_IR/L_UV). This relation allows the estimation of the total IR luminosity from the shape of the ultraviolet spectral region, providing a way to quantify the amount of attenuation of UV starlight (i.e., A_UV) reprocessed by dust. It has been abundantly used to derive the SFR of high redshift galaxies, for which the UV slope is the only accessible quantity (Reddy et al. 2006; Smit et al. 2012).

To account for departures of quiescent star-forming galaxies from the star-burst IRX-β relation, several authors (e.g., Kong et al. 2004; Seibert et al. 2005; Salim et al. 2007) proposed a modified version of this relation, expressed in terms of A_FUV − β, which leads to a smaller dust correction at a given slope β, though with a significant scatter (σ(A_FUV) ~ 0.9). It is worth mentioning that such a large scatter is not surprising in view of the steepness of the A_FUV − β relation, where a small uncertainty in UV color (e.g., Δ(FUV − NUV) ~ 0.1) produces a large variation in A_FUV (ΔA_FUV ~ 0.8). Johnson et al. (2007a) proposed another method calibrated for quiescent star-forming galaxies which relies on the combination of the D₄₀₀₀ break, a spectral feature sensitive to the age of a galaxy, and a long baseline color. They explored different colors and functions and obtain the smallest residuals by using the (NUV − 3.6  μm) color (σ(IRX) ~ 0.3 for their star-forming population). Treyer et al. (2007) compared the SDSS SFR estimates of Brinchmann et al. (2004), based on nebular recombination lines, with those obtained from the UV continuum corrected for the effect of dust attenuation with the methods of Seibert et al. (2005); Salim et al. (2007), and of Johnson et al. (2007a). They found that the method of Johnson et al. (2007a) leads to the smallest scatter [σ(SFR) ~ 0.22 vs. 0.33]. However, the method of Johnson et al. (2007a) depends on spectral diagnostics, such as the D₄₀₀₀ break, which are difficult to obtain for high redshift galaxies.

In this work, we analyze the behavior of the IRX within the rest-frame color–color diagram (NUV − r) vs. (r − K). We describe a new relationship between IRX and a single vector, called NRK , defined as the combination of the colors (NUV − r) and (r − K). This new diagnostic provides an effective way to assess the total IR luminosity L_IR and then the SFR for individual galaxies with photometric information widely available in large surveys, and does not require complex modeling such as SED fitting techniques. The (NUV − r) vs. (r − K) color diagram adopted in this work is similar to the (U − V) vs. (V − J) diagram proposed by Williams et al. (2009) to separate passive or quiescent from star-forming galaxies, but better leverages the role of SFH and dust by extending into the extreme wavelengths of the SED. Total infrared luminosities were estimated from the deep MIPS-24 μm observations of the COSMOS field, with standard technics extrapolating mid-IR flux densities with luminosity-dependent IR galaxy SED templates as described in e.g., Le Floc’h et al. (2005). Our approach is motivated by the tight correlations that exist between mid-IR and total IR luminosities, both at low redshifts (Chary & Elbaz 2001; Takeuchi et al. 2003) and in the more distant Universe (Bavouzet et al. 2008). In fact, stacking analysis and direct individual identifications of distant sources in the far-IR with facilities like Spitzer and Herschel (e.g. Lee et al. 2010; Elbaz et al. 2010) have shown that these extrapolations from mid-IR photometry provide reliable estimates of total IR luminosities up to z ~ 1.3 (dispersion of ~0.25 dex, no systematic offset), at least for star-forming galaxies initially selected at 24 μm. Also, given the relative depths of the different IR and submillimeter observations carried out in the COSMOS field with e.g., Spitzer, Herschel or JCMT, we note that the deep 24 μm data in COSMOS provide the largest sample of star-forming galaxies with measurable IR luminosities. We thus decide to limit our analysis to massive (M ≥ 10^9.5 M_⊙) star-forming sources first selected at 24 μm and lying at z ≤ 1.3 so as to obtain reliable L_IR estimates. In a companion paper Le Floc’h et al. (in prep.), we will extend our analysis to higher redshifts and lower masses by stacking along the NRK vector the MIPS and Herschel data available in the COSMOS field.

The paper is organized as follows. In Sect. 2 we describe the dataset, the sample selection and the estimates of physical parameters. In Sect. 3 we discuss the behavior of the mean IRX (i.e., ⟨IRX⟩) in the (NUV − r) vs. (r − K) diagram and the calibration of the ⟨IRX⟩ vs. NRK relation. In Sect. 4, we compare our predicted IR luminosities with the ones based on the 24 μm observations and SED fitting technique. We also investigate the dependence with galaxy physical parameters and apply our method to the well established SFR-mass relation. In Sect. 5 we discuss the method, the possible origin of the relation with a complete library of model galaxy SEDs and dust models. We draw the conclusions in Sect. 6.

Throughout the paper we adopt the following cosmology: H₀ = 70 km s^-1 Mpc^-1 and Ω_M = 0.3, Ω_Λ = 0.7. We adopt the initial mass function of Chabrier (2003), truncated at 0.1 and 100 M_⊙. All magnitudes are given in the AB system (Oke 1974).

2. The dataset

2.1. Observations

2.1.1. The Spitzer-MIPS 24 μm observations and COSMOS catalog

The region of the sky covered by the Cosmic Evolution Survey (i.e., COSMOS Scoville et al. 2007) is observed by Spitzer-MIPS (Multi Imaging Photometer) at 24, 70 and 160  μm over a 2 deg² area (Sanders et al. 2007). In this work, we use the deep 24  μm observations and the catalogue of sources extracted by Le Floc’h et al. (2009). This catalog is 90% complete down to the flux limit density S₂₄ ~ 80  μJy. The sources are cross-matched with the multi-wavelength photometric COSMOS catalog. This includes deep ultraviolet GALEX imaging (Zamojski et al. 2007), ground based optical observations with intermediate and broad band filters, near IR photometry (Capak et al. 2007; McCracken et al. 2010; Taniguchi et al. 2007) and deep IRAC photometry (Sanders et al. 2007; Ilbert et al. 2010) for a total of 31 filter passbands. The combined photometry, as well as accurate photometric redshifts, are available from the catalog of Ilbert et al. (2010, version 1.8). The typical depths at 5σ are 25.2, 26.5, 26.0, 26.0, 25.0, 23.5 and 24.0 in the NUV, u, g, r, i, K and 3.6  μm pass-bands, respectively. A vast majority of 24  μm sources is also detected in the optical bands (~95% with u ≤ 26.5 and i_AB ≤ 24.5).

We adopt spectroscopic redshifts from the bright and faint zCOSMOS sample (Lilly et al. 2009), when available, otherwise we use the photometric redshifts computed by Ilbert et al. (2009). The photo-z accuracy for the 24  μm sample is σ = 0.009 at i_AB ≤ 22.5 (with 5700 z_spec) and 0.022 for the fainter sample (with 470 z_spec)¹.

We limit our analysis to a maximum redshift z ~ 1.3 in order to derive reliable IR luminosities from the 24  μm flux density (see next section). We reject AGN dominated sources according to their mid-IR properties based on the diagnostic of Stern et al. (2007). The selected sample consist of ~16 500 star-forming galaxies and ~400 evolved/passive galaxies, with S₂₄ ≥ 80  μJy and z ≤ 1.3. The separation between star-forming and passive galaxies is based on the position of the galaxies in the rest-frame (NUV − r) versus (r − K) diagram as discussed in Sect. 2.3 and Appendix B.

2.1.2. A low-z sample

We complement the COSMOS field with a lower redshift sample selected from Johnson et al. (2007a), which includes SWIRE observations (Lonsdale et al. 2003) in the Lockman Hole and the FLS regions. The multi-wavelength observations combine the GALEX, SDSS and 2Mass photometry. We restrict the sample to the sources detected in the three MIPS passbands (24, 70 and 160  μm). The final sample consists of ~730 galaxies with z ≤ 0.2 (mean redshift: ) of which ~560 are star-forming galaxies (with the same separation criterion as above).

2.2. Rest-frame quantities and physical parameters

2.2.1. The infrared luminosity

The total Infrared luminosity (L_IR) refers to the luminosity integrated from 8 to 1000  μm and is derived by using the code “Le Phare”² (Arnouts et al. 1999; Ilbert et al. 2006) combined with infrared SED templates of Dale & Helou (2002). For the low-z sample, the IR luminosity is estimated by fitting the 8, 24, 70 and 160  μm bands with a free scaling of the SEDs (see Goto et al. 2011, for details). For the COSMOS sample, the L_IR is derived by extrapolating the 24  μm observed flux density with the star-forming galaxy templates of Dale & Helou (2002). While such SEDs are provided as a function of 60/100  μm luminosities ratio, we rescale the templates following the locally-observed dust temperature-luminosity relationship, so as to mimic a luminosity dependence similar to the SEDs of the libraries proposed by Chary & Elbaz (2001) and Lagache et al. (2004). In this scheme, the monochromatic luminosity at any given IR wavelength is linked to L_IR by a monotonic relation, similar to the correlations that have been observed between the mid-IR emission and the total IR luminosity of galaxies in the local Universe (e.g., Spinoglio et al. 1995; Chary & Elbaz 2001; Takeuchi et al. 2003; Treyer et al. 2010; Goto et al. 2011). In this way, a 24  μm flux density observed at a given redshift is associated to a unique L_IR. Note that we did not consider any AGN SEDs in this work, since AGN-dominated sources were removed from our initial sample (see Sect. 2.1.1).

Such extrapolations from mid-IR photometry have been widely used in the past, especially for the interpretation of the mid-IR deep field observations carried out with the infrared Space Observatory and the Spitzer Space Telescope. Given the large spectral energy distribution (SED) variations between the rest-frame emission probed at 24  μm and the peak of the IR SED where the bulk of the galaxy luminosity is produced, the estimate of L_IR with this method yet depends on the assumed SED library. To quantify this effect, we use different IR star-forming galaxy templates of the literature, finding systematic differences of only ≲0.2 dex up to z ~ 1 (e.g., see Fig. 7 of Le Floc’h et al. 2005). Furthermore, stacking analysis of mid-IR selected sources with MIPS-70 μm and MIPS-160 μm data in COSMOS (Lee et al. 2010) and other fields (e.g., Bavouzet et al. 2008) has allowed accurate determinations of average IR luminosities for galaxies stacked by bin of 24 μm flux. These studies have shown very good agreement with the 24 μm extrapolated luminosities of star-forming galaxies up to z ~ 1.5, thus confirming the robustness of the technic. Even more convincingly, direct far-IR measurements of the total IR luminosities for a 24 μm flux-limited sample of star-forming galaxies at z ≤ 1.5 with the Herschel Space Observatory have recently revealed a tight correlation with the luminosities extrapolated from 24  μm photometry and star-forming SED templates (Elbaz et al. 2010). The comparison between the two estimates shows a dispersion less than 0.15 dex at z ≤ 1, demonstrating again the relevance of the method at least for mid-IR selected sources at low and intermediate redshifts, such as in our current study. For these reasons we restrict our analysis to redshifts lower than z ~ 1.3, where the method discussed above to derive the L_IR is widely tested and robust.

2.2.2. The other physical parameters and luminosities

We use the code “Le Phare” combined with the population synthesis code of Bruzual & Charlot (2003, hereafter BC03) to derive the physical parameters for each galaxy. Details regarding the SED fitting are given in Appendix A. We perform a maximum likelihood analysis, assuming independent Gaussian distributed errors, including all the available photometric bands from 0.15 to 4.5  μm. The physical parameters are derived by considering the median value of the likelihood marginalized over each parameter, while errors correspond to the 68% credible region.

We adopt the same approach as Ilbert et al. (2006) to derive the rest-frame luminosities (or absolute magnitudes). We use the photometry in the nearest rest-frame broadband filter to minimize the dependency to the k-correction. These luminosities are consistent within 10% with those derived according to the best fit template (smallest χ²) or those provided by Ilbert et al. (2009) based on a smaller set of empirical SEDs. Our results are not significantly affected by the adopted choices of luminosities³.

We note also that the great accuracy of the COSMOS photometric redshifts in our redshift domain allows us to neglect the impact of the photo-z errors in the quantities discussed here. However, we adopt spectroscopic redshifts when available.

Throughout this paper, the stellar masses always refer to the estimates from the SED fitting technique, assuming a Chabrier IMF truncated at 0.1 and 100 M_⊙, while the total SFR is defined as the sum of the unobscured ultraviolet and total IR luminosities. For the latter, we adopt a relation similar to that proposed by Bell et al. (2005) and adjusted for a Chabrier (2003) IMF: (1)where the total IR luminosity is defined as and ℒ_NUV is the monochromatic NUV luminosity: ℒ_NUV/L_⊙ = νL_ν(2300  Å). We adopt a factor of ~2.3 instead of 1.9 as proposed by Bell et al. (2005) who uses the FUV luminosity. According to stellar population models, our factor is more appropriated to correct the L(2300 Å) in total UV luminosity. The use of the FUV luminosity better traces the emission of short-lived, massive stars, while in this work we use the NUV luminosity. The reason is that we can obtain a more accurate rest-frame NUV than FUV luminosity, thanks to the GALEX NUV and CFHT u-band observations at different redshifts. As shown by Hao et al. (2011), this choice does not impact the reliability of the SFR estimates.

It is worth noting that, as shown in Appendix A, the instantaneous SFR derived from the SED fitting is consistent within less than a factor of two with the SFR estimated with Eq. (1). Considering that the SED fitting relies only on the 0.15 − 4.5  μm bands, this agreement between the different SFR estimates over at least 2 orders of magnitudes is remarkable.

Figure 1 shows the SFR (top panel; as estimated with Eq. (1)) and stellar mass (bottom panel) as a function of redshift for the 24  μm star-forming galaxies. As already shown by Le Floc’h et al. (2005), below z ≤ 0.5 the population is composed of moderately star-forming galaxies with SFR ≤ 10 M_⊙ yr^-1. The fraction of luminous IR galaxies (LIRGs) gradually increases from z = 0.5 to z = 1 and becomes dominant at z ≳ 1 in our sample. At all redshift, the fraction of ultra luminous IR galaxies (ULIRGs) is negligible.

Our sample is dominated by galaxies with log   (M/M_⊙) ≳ 9.5. To characterize how representative the 24  μm population is with respect to the entire star-forming (hereafter SF) population at a given mass and redshift, we define a 50% completeness mass limit (M_50%(z)) as the stellar mass where the ratio / ~ 0.5; where and are the V_max weighted comoving volume densities of 24  μm- and K-selected (K ≤ 23.5) samples of star-forming galaxies respectively. Above this limit (~10¹⁰ [10^10.5]  M_⊙ at z ~ 1.1 [1.3]), we consider the physical properties of the 24  μm population to be representative of the whole SF sample.

2.3. The (NUV – r) versus (r – K) color–color diagram

Recently, Williams et al. (2009) have shown that quiescent and star-forming galaxies occupy two distinct regions in the rest-frame (U − V) versus (V − J) color–color diagram (i.e., UVJ) and validated their separation scheme with a morphological analysis (Patel et al. 2011). In the present work, we adopt the (NUV − r) versus (r − K) color–color (i.e., NUVrK ) diagram to increase the wavelength leverage between current star formation activity and dust reddening. The NUVrK diagram provides also an efficient way to separate passive and star-forming galaxies. In Fig. 2, we show the mean sSFR derived from the SED fitting for the entire COSMOS spectroscopic sample with 0.2 ≤ z ≤ 1.3. The density contour levels (black solid lines) reveal the presence in the top left part of the diagram of a population with low specific SFR [log   (sSFR^SED/yr^-1) ≲ − 10.5], well separated from the rest of the sample. We define this region with the following criteria: (NUV − r) > 3.75 and (NUV − r) > 1.37 × (r − K) + 3.2. In Appendix B we discuss further the star-forming and passive galaxy separation based on the morphological information from HST imaging of the COSMOS field.

In Fig. 3, we show the density distribution of the 24 μm sample in this NUVrK diagram and the color tracks for five models with exponentially declining star formation histories. The colors indicate models with different e-folding times τ = 0.1, 1, 3, 5, 30 Gyr (red, purple, orange, green, blue lines respectively), while the filled and open squares mark different ages (t = 0.1, 0.5, 1, 3, 5 and t = 6.5, 9.0, 12 Gyr, respectively).

The (NUV − r) is a good tracer of the specific SFR, since the NUV band is sensitive to recent (i.e., t ≤ 10^8.5 yr) star formation and the r-band to old stellar populations (e.g., Salim et al. 2005). Models with short star formation timescales quickly move to the top side of the diagram, since star formation ceases early and the integrated light becomes dominated by old, red stars. On the other hand, models with longer e-folding times show bluer colors at all ages, since they experience a more continuous star formation which replenishes the galaxy with young, hot stars emitting in the UV. This behavior has been widely used to separate the active and passive galaxy populations based on (NUV − r) vs. stellar mass (or luminosity) diagrams (e.g., Martin et al. 2007; Salim et al. 2005). However, dust attenuation in star-forming galaxies can strongly alter the (NUV − r) color, producing variations up to several magnitudes [e.g., Δ(NUV − R) ~ 2 for E(B − V) ~ 0.4 and SMC law]. Following Williams et al. (2009) we use a second color (r − K), which does not vary significantly with the underlying stellar population, even for passive galaxies [while (NUV − r) does], to partially break this degeneracy. As shown in Fig. 3, the models of spectral evolution span a much smaller range of (r − K) color than that observed in the 24  μm sample, unless dust attenuation is included. We can qualitatively reproduce the observed ranges of (NUV − r) and (r − K) colors by applying a reddening excess E(B − V) ≤ 0.4 − 0.6, assuming a starburst attenuation law or an SMC-like extinction curve.

3. Infrared excess of star-forming galaxies in the NUVrK diagram

3.1. Definition of the IRX and NRK parameters

We now explore the behavior of the infrared excess in the NUVrK diagram. We define the infrared excess as IRX = L_IR/ℒ_NUV, where L_IR and ℒ_NUV are defined in Sect. 2.2.2. This ratio is weakly dependent on the age of the stellar population, dust geometry and nature of the extinction law (Witt & Gordon 2000).

Figure 4 shows on the left, the volume weighted mean IRX (⟨IRX⟩, color-coded in a logarithmic scale shown in the top right panel), and, on the right, the dispersion around the mean (σ(IRX)) in the NUVrK diagram in four redshift bins. At any redshift, ⟨IRX⟩ increases by ~1.5–2 dex from the blue side (bottom-left) to the red side (top-right) of the diagram, while the dispersion around the mean remains small and constant (≤0.3 dex). We note also the presence of stripes of constant ⟨IRX⟩ values, which we discuss in Sect. 5. The presence of such stripes allows us to describe the variation of ⟨IRX⟩ in the NUVrK diagram with a single vector perpendicular to those stripes. This vector, hereafter called NRK , can be defined as a linear combination of rest-frame colors NRK(φ) = sin(φ) × (NUV − r) + cos(φ) × (r − K), where φ is an adjustable parameter which we require to be perpendicular to the ⟨IRX⟩ stripes. We empirically estimate φ in each redshift bin (Δz = 0.1) by performing a linear least square fit: IRX = f [NRK(φ)] (the linear approximation is justified in the next section). Since the dispersion σ [IRX(φ)] reaches a minimum when the vector NRK(φ) is perpendicular to the stripes, we minimize σ [IRX(φ)] to find the best-fit angle, φ_b. We find 15° ≤ φ_b ≤ 25° in all redshift bins, with a median and semi-quartile range φ_b = 18° ± 4°. Given the narrow distribution of φ_b and the fact that a change up to Δφ ± 7° does not affect our results, we adopt a unique definition for the vector NRK at all redshifts by fixing φ_b = 18°, (2)In Fig. 4 we show the vector NRK (black arrows) and a number of lines perpendicular to it (gray lines) corresponding to constant value of NRK in the range 0 ≤ NRK ≤ 4. We note that NRK has a different orientation in the NUVrK diagram with respect to the starburst and SMC dust attenuation vectors, as we discuss in Sect. 5.

3.2. The calibration: ⟨IRX⟩ versus NRK

With the definition of NRK , we can now derive the relationship between ⟨IRX⟩ and NRK. In Fig. 5 we measure the mean IRX values (⟨IRX⟩) and the associated dispersions per bin of NRK for the 24  μm sample (solid black circles) in different redshift bins. At high redshift, a tight correlation is observed with a small dispersion (σ ≤ 0.3) compared to the evolution of ⟨IRX⟩ (Δ ⟨  IRX  ⟩  ~ 2 dex). At z ≥ 0.4, we obtain almost the same results for the spectroscopic sample (gray squares with yellow shaded region). Although this sample is 1/10 of the 24  μm sample, we observe the same dispersions, suggesting that we are dominated by an intrinsic, physical scatter rather than poisson noise.

In the first panel (0.2 ≤ z ≤ 0.4), we also include the lower redshift sample from Johnson et al. (2007a) (, gray squares with yellow shaded region). While the ⟨IRX⟩ vs. NRK relation remains the same as at higher redshift, the dispersion increases in both of the low-z samples. This is related to the increasing fraction of less active galaxies with a higher contribution of evolved stars affecting the IR luminosity and contributing to blur the correlation between IRX and NRK (e.g., Cortese et al. 2008).

Finally we measure the ⟨IRX⟩ vs. NRK relation for three stellar mass bins. We observe no significant difference in the relation derived in each bin and the global one. This, along with the SED reconstruction analysis shown in Appendix C, supports our assumption of neglecting the dependence on stellar mass in the calibration of ⟨IRX⟩ vs. NRK in the mass range considered here.

In the parametrization of ⟨IRX⟩ as a function of NRK and redshift, we assume that the two quantities can be separated: (3)where f(z) is a third-order polynomial function describing the redshift evolution, and a_N a constant which describes the evolution with NRK. We bin the data in redshift (Δz = 0.1) and NRK (ΔNRK = 0.5) and perform a linear least square fit to derive the five free parameters. We obtain for the redshift evolution f(z) = a₀ + a₁.z + a₂.z² + a₃.z³ with a₀ = −0.69 ± 0.06; a₁ = 3.43 ± 0.33; a₂ = −3.49 ± 0.55; a₃ = 1.22 ± 0.28, and for the dependence on NRK a_N = 0.63 ± 0.01.

The resulting fits from this calibration are shown in Fig. 5 (solid red lines). We show, in the top panel of Fig. 6, the linear fit of ⟨IRX⟩ as a function of NRK after rescaling all the ⟨IRX⟩ values at z = 0 and in the bottom panel, the polynomial fit of ⟨IRX⟩ vs. z after rescaling all the ⟨IRX⟩ values at NRK = 0 The small uncertainty in the slope parameter a_N supports our initial choice of a linear function to describe the relation between ⟨IRX⟩ and NRK. We note that due to adoption of a polynomial function, the redshift evolution is only valid in the range well constrained by the data (0.1 ≤ z ≤ 1.3) and should not be extrapolated outside this range.

A simple interpretation of the redshift evolution of ⟨IRX⟩, at fixed NRK, is the aging of the stellar populations. In fact, Fig. 3 shows that, at a fixed position in the NUVrK diagram (or NRK value), a galaxy at higher redshift, which has a younger stellar population, needs a larger dust reddening than a galaxy at lower redshift, which hosts older, intrinsically redder stars. This effect being more pronounced between 0 ≤ z ≤ 0.5 where the universe ages by ΔT ~ 5 Gyr, compared to 4 Gyr in between 0.5 ≤ z ≤ 1.3, and which is also enhanced by the global decline of the SF activity in galaxies with cosmic time.

4. The infrared luminosity and SFR estimated from NRK vector

The relationship between IRX and the vector NRK allows us to predict the IR luminosity for each galaxy according to its NUV, r, K luminosities and redshift as follows , where NRK and ⟨IRX⟩ are estimated from Eqs. (2) and (3), respectively.

4.1. Comparison with the reference LIR

In Fig. 7 (left and middle panels) we compare the IR luminosity predicted with the NRK method () with our reference IR luminosity derived from the mid/Far-IR bands (L_IR, see Sect. 2.2.1). The IR luminosities estimated with the two methods for the ~16 500 star-forming galaxies in the COSMOS 24  μm sample (left panel) agrees with almost no bias and a dispersion of ~0.2 dex over the entire luminosity range. Less than 1% of the galaxies shows a difference larger than a factor of 3. The prediction of the IR luminosity for the LIRG population, which dominates the COSMOS sample (the 24  μm population peaks at L_IR ~ 2 × 10¹¹ L_⊙), is excellent, considering the small number of parameters involved in the NRK method.

At low redshift, the prediction of the IR luminosity for the SWIRE galaxies (middle panel), which are 10 times less luminous than the 24  μm COSMOS population (L_IR ~ 2 × 10¹⁰ L_⊙), shows a larger scatter (σ ~ 0.27 dex). This reflects the larger dispersion observed in Fig. 5 which is likely related to the wider range of galaxy properties at low z. Johnson et al. (2007b) have also modeled the IRX as a function of different rest-frame colors and D_n(4000) break for the SWIRE galaxies (see their Table 2 and Eq. (2)). The fit residuals for their predictions of IRX vary from σ(IRX) ~ 0.27 to 0.36, depending on the adopted colors. Even with the use of the D_n(4000) break as a dust-free indicator of star formation history, they do not achieve a better accuracy than that obtained with the method presented in this work.

Finally, our method provides a smaller dispersion in the prediction of the IR luminosity than the dust luminosity obtained from the SED fitting (Fig. 7 right panel), where we obtain a dispersion of ~0.26 dex and a bias of Δ ~ − 0.13. We note that the L_IR derived from SED fitting is computed by integrating all the stellar photons absorbed by dust according to the adopted attenuation law and reddening excess (see Appendix A), while the reference L_IR is based on the extrapolation of the 24  μm luminosity. Both measurements suffer from independent source of uncertainties, thus the bias in the IR luminosity may be related to either methods. However, the small scatter (less than a factor 2) and bias show that the SED fitting is a robust and reliable method to estimate the total IR luminosity, when accurate observations on a wide enough wavelength are available (e.g., Salim et al. 2009). This is indeed the case for the COSMOS dataset used here (i.e., 31 passbands from Far-UV to Mid-IR: 0.15 ≤ λ ≤ 4.5  μm).

4.2. Dependence of the predicted IR luminosity with physical galaxy properties

Despite the use of a single vector, NRK , and the absence of mass dependency, our recipe provides a reasonable estimate of the total IR luminosity L_IR. However, galaxies with physical properties that deviate from the bulk of the 24  μm population may be less accurately described in this framework. To test this issue and determine the range of validity of the method, we show in Fig. 8 the residuals (i.e., the difference) between the reference (i.e., ) and predicted IR luminosity as a function of redshift (top-left panel), total SFR (top-right), stellar mass (bottom-left) and specific SFR (bottom-right). The total SFR refers to the (UV+IR) SFR (see Eq. (1)), while the stellar mass is obtained with the SED fitting (see Sect. 2.2.2). We perform this comparison for the L_IR predicted with both the NRK and SED fitting methods.

Top-left panel of Fig. 8 does not reveal any bias with redshift in the IR luminosity predicted with both methods. This is expected by construction for the NRK method, while it confirms the good performance of the SED fitting technique despite the larger scatter than the one obtained with the NRK method.

Top-right panel of Fig. 8 shows that the residuals as a function of the total SFR are relatively stable, with an almost zero bias and a dispersion of ~0.2 dex for SFR ≤ 100 M_⊙ yr^-1. At higher SFR, where galaxies approach the ULIRG regime, the residuals start to deviate from zero, indicating that both the NRK and SED fitting methods under-predict the reference IR luminosity. The difference with the NRK predictions is likely to be caused by the rarity of these highly star-forming objects in our redshift range (0 ≤ z ≤ 1.3, see also Le Floc’h et al. 2005), which makes these objects under-represented in the NRK calibration. The disagreement with the prediction of the SED fitting may have a different origin. As shown in Appendix A, the starburst attenuation law is favored at high SFR, however ULIRGs are found to deviate from Meurer’s relation (Reddy 2009; Reddy et al. 2012) and exhibit a higher IRX at a given UV-slope. This effect has been observed at higher redshift (z ≥ 1.5), but it could also be the reason of the under-estimate of L_IR by the SED fitting technique at lower redshift. We also note that our e-folding SF history models used in the SED fitting may be too simplistic and models including burst episodes would be more appropriated for these actively star-forming galaxies. Finally, we also note that the reference L_IR, based on the extrapolation of the 24  μm luminosity, may reach its limit of validity in this regime. Indeed, the conversion of the 24  μm luminosity to total IR luminiosity becomes uncertain for ULIRGs (i.e., for L_IR ≥ 10¹² L_⊙, see Bavouzet et al. 2008; Goto et al. 2011). Also, the merger nature of ULIRGs at low z (Kartaltepe et al. 2010) makes predictions inaccurate in absence of a complete set of far-IR observations.

Bottom-left panel of Fig. 8 does not reveal any bias with stellar mass for the IR luminosity predicted with the SED fitting technique, while the dispersion increases at the extreme sides. On the other hand, the NRK method under-estimates L_IR for stellar mass log   (M/M_⊙) ≲ 9.3. As shown in Le Floc’h et al. (in prep.) from the analysis of a mass-complete sample obtained with stacking techniques of Herschel/SPIRE data at 250/350/500  μm, the NRK calibration presented in this work needs to be modified for galaxies with M ≲ 10^9.5 M_⊙. At high stellar mass (i.e., log   (M/M_⊙) ≳ 11) the NRK method over-estimates L_IR. This reflects the correlation between stellar mass and sSFR, where quiescent galaxies become the dominant population at high stellar masses (see below).

Bottom-right panel of Fig. 8 shows that the NRK method systematically (as reflected by the small scatter) over-estimates L_IR at low sSFR (i.e., log   (sSFR/yr^-1) ≤ − 10). We investigate this bias in more details in Appendix D and find that the sSFR derived with NRK saturates at such low sSFR, while the reference sSFR keeps decreasing. We therefore propose a simple analytical correction for this bias in order to reconcile the two sSFR estimates, which can then be used to correct . The results for this sSFR-corrected are shown in Fig. 8 as filled green squares. The correction provides a better match to the 24  μm parent distribution across the entire range of stellar mass and sSFR, but at the cost of a larger dispersion for the global sample (i.e., from σ = 0.21 to σ = 0.24). It is worth noting that the reference sSFR for quiescent galaxies might also suffer from systematic error. In fact, recent analysis with Herschel data have reported evidences for a warmer dust temperature in early- than in late-type galaxies (Skibba et al. 2011; Smith et al. 2012). This will result in an over-estimation of L_IR as derived from the 24  μm luminosity (by adopting a too high L_IR/L_24 μ ratio). Combined with the UV luminosity partially produced by evolved stars, these two effects can lead to an over-estimate of our reference SFR (based on UV+IR). As seen in Fig. 8, the L_IR comparison with the SED fitting method shows a large dispersion at low sSFR. The reason is illustrated in Appendix A (Fig. A.1), where a significative fraction of the galaxies with log   (sSFR_tot/yr^-1) ≤ − 10 shows a much lower sSFR with the SED fitting method, leading to a lower L_dust.

In conclusion, the NRK and SED fitting methods provide two independent and reliable estimates of L_IR over a large range of redshift and galaxy physical parameters for the vast majority (~90%) of the 24  μm star-forming sample. However, differences arise when considering the extreme sides of the population, such as highly star-forming (i.e., SFR_tot ≥ 100 M_⊙ yr^-1) or quiescent galaxies (i.e., log   (sSFR_tot/yr^-1) ≤ − 10), which represent ≲1% and 7% of the whole sample, respectively. While the NRK method suffers from a bias at low sSFR, it is not clear the amplitude of this effect, since there is no robust indicators of SFR in the low activity regime. In the next section, we obtain an independent estimate of this bias for the quiescent population by directly comparing the SFRs estimated with the NRK and SED fitting methods.

4.3. Comparison of SFR estimates with NRK and SED fitting techniques

Once calibrated with a Far-IR sample, the NRK technique can be applied to any sample of galaxies. In this section we compare the SFRs derived with the SED fitting and NRK methods for two samples of galaxies, the 24  μm and a K-selected (down to K ≤ 23.5) samples. We restrict the analysis to galaxies with stellar masses M ≳ 2 × 10⁹ M_⊙ and redshift 0.05 ≤ z ≤ 1.3. The SFR^NRK is estimated by means of Eqs. (1)–(3). As discussed above, we found that this technique is better suited for active star-forming galaxies than the less active population (low sSFR). As shown in Fig. 2, we exploit the capability of the NUVrK diagram to separate galaxies with different sSFRs and we define a new vector NRK_ssfr, perpendicular to the starburst attenuation vector (NRK_ssfr = cos(54°) × (NUV − r) − sin(54°) × (r − K)). In Fig. 2, we show constant values of NRK_ssfr (i.e., the norm of the NRK_ssfr vector), corresponding to the range − 2 ≤ NRK_ssfr ≤ 3 as gray lines, which is a good proxy to follow the variation of the mean sSFR. The limit adopted to define the passive population corresponds to NRK_ssfr ≳ 1.9.

The ratios (SFR^NRK/SFR^SED) for the 24  μm (top panel) and the K (bottom panel) selected samples are shown in Fig. 9 for different intervals of NRK_ssfr. In each panel, we report for each bin of NRK_ssfr the percentage of objects in this bin, the percentage of catastrophic objects (estimated as ABS(Δ) > 3σ), the median (Δ) and scatter (σ). For ~85% of the sample, with NRK_ssfr  ≤  1.3, the two SFRs are in excellent agreement, with small biases, scatters and catastrophic fraction for both the 24  μm and K-selected samples. For 1.3 ≤ NRK_ssfr ≤ 1.5 (~6–7% of the two samples), the fraction of catastrophic objects sharply increases to ~10%, the dispersion also increases and the median starts to deviate from zero, with SFR^SED being higher than SFR^NRK. This trend becomes more severe for 1.5 ≤ NRK_ssfr ≤ 1.9 (the remaining ~7 to 9% of the samples). The two SFRs disagree, with a large fraction of catastrophic objects (≳25%), a significant non-zero median and large dispersion, with some slight differences between the two samples. This region with inconsistent SFR estimates is shown in Fig. 2 as a shaded area.

This comparison shows that the NRK method can be successfully applied to a sample of galaxies larger than the 24  μm sample, since it provides SFR estimates in good agreement with the SED fitting for the vast majority (i.e., ~85–90%) of the star-forming population. However, the methods diverge when including low-activity (low-ssFR) galaxies (see also Johnson et al. 2007a; Treyer et al. 2007). To alleviate this limitation, we show that the population with discording results (≤10% of the whole sample) can be easily isolated in the NUVrK diagram.

4.4. Application to the SFR vs. stellar mass relationship

A correlation between SFR and stellar mass in star-forming galaxies is observed both at low (Salim et al. 2007) and high redshift (Noeske et al. 2007b; Elbaz et al. 2007). This relation is surprisingly tight, with an intrinsic scatter of ~0.3 dex. We can therefore exploit the presence of this correlation and test the ability of the NRK method to reproduce the slope and amplitude of this relation at different redshifts.

In Fig. 10 we show the V_max weighted mean of SFR in bins of stellar mass for different samples of star-forming galaxies: the COSMOS 24  μm sample, a flux limited K-selected sample (K ≤ 23.5) and, in the lowest redshift bin, the local GALEX-SWIRE sample. Symbols indicate SFR estimated with different methods: (L_IR + L_UV) from Eq. (1) (black stars), NRK method (red filled circles and open blue squares), NRK method corrected for the sSFR bias (green-filled squares).

At low-z, the different SFR estimates provide similar results and the SFR-M_∗ relation is in excellent agreement with the one measured by Salim et al. (2007) for the GALEX-SDSS sample and derived with a SED fitting method (solid red line). At higher redshifts, we compare our finding with the results of Karim et al. (2011), which are based on the 1.4 GHz radio continuum emission of stacked star-forming galaxies in different stellar mass and redshift bins (solid blue lines). We consider their results from the two-parameter fits given in their Table 4. The SFR-M_∗ relation derived from individual galaxies with both the NRK sSFR bias-uncorrected and corrected methods agrees well with the radio stacking from Karim et al. (2011). The major difference between the NRK and the NRK-sSFR corrected method is in the dispersions around the mean SFRs. The original NRK method shows a small dispersion in the high mass end of the relation, due to an over-estimate of the SFR for massive, evolved galaxies, artificially moving them towards the SFR-M_∗ sequence. This effect varies with redshift, as described in Appendix D, thus we indicate with the shaded area in Fig. 10 the regions in which the SFR is over-estimated by a factor greater than two. It can be seen that, at all redshifts, the bin corresponding to the most massive galaxies is affected by this bias, hence producing a smaller dispersion on the SFR-M_∗ relation with respect to that observed with the 24  μm and NRK-sSFR corrected methods. We note that this effect becomes irrelevant at lower masses.

Overall, even if the stellar mass does not enter in the calibration of the NRK method, the SFR estimated with this method for individual galaxies can reproduce the slope and normalisation of the SFR-M_∗ relation, along with its redshift evolution.

5. Modeling the ⟨IRX⟩ in the NUVrK diagram

In the previous section we have shown that galaxies with different ultraviolet-to-infrared luminosity ratios (IRX) are well separated in the NUVrK diagram. To confirm the validity of this approach, and to explain the physical origin of the observed trends, we appeal to a library of galaxy SEDs computed with the BC03 spectral evolution model. We follow the approach of Pacifici et al. (2012) and extract a set of star formation and chemical enrichment histories from the semi-analytic post-treatment of De Lucia & Blaizot (2007) of the Millennium cosmological simulation (Springel et al. 2005). The star formation and chemical enrichment histories computed in this way reproduce the mean properties of nearby SDSS galaxies. Thus, they span only limited ranges in sSFR, around 10^-10 yr^-1, and in the fraction of the current stellar mass formed in the last 2.5 Gyr at z = 0. Following Pacifici et al. (2012), to account for the broader range of spectral properties of the galaxies in our sample with respect to the SDSS, we re-draw the evolutionary stage at which a galaxy is looked at in the library of star formation and chemical enrichment histories (we do this uniformly in redshift between 0.2 and 1.5). We also resample the current (i.e., averaged over the last 30 Myr) SFR from a Gaussian distribution centered on log   (sSFR) = −9.1, with a dispersion of 0.6. This choice of parameters reproduces the observed global distribution of sSFR (i.e., summed over all redshift bins) and the distribution of galaxies in the NUVrK diagram, after accounting for dust attenuation as described in the next section. We adopt the Chabrier (2003) IMF.

5.1. Dust attenuation model

To include the effect of dust attenuation, we adopt the dust prescription of Chevallard et al. (2013). This extends the two-components, angle-average dust model of Charlot & Fall (2000) to include the effect of galaxy inclination and different spatial distributions of dust and stars on the observed SEDs of galaxies. Chevallard et al. (2013) combine the radiative transfer model of Tuffs et al. (2004, hereafter T04) with the BC03 spectral evolution model. To accomplish this, they relate the different geometric components of the T04 model (a thick and thin stellar disks, and a bulge, attenuated by two dust disks) to stars in different age ranges. Here, for the sake of simplicity and to limit the number of adjustable parameters, we describe attenuation in the diffuse ISM using only the thin stellar disk model of Tuffs et al. (2004). This is supported by the finding by Chevallard et al. (2013) that, in a large sample of nearby star-forming galaxies, the thin stellar disk component of the T04 model accounts for ≈80 percent of the attenuation in the diffuse ISM. Also, we note that adding the T04 thick disk component has a weak effect on the results. The dust content of the diffuse ISM in the T04 model is parametrized by means of the B-band central face-on optical depth of the dust disks τ_B⊥ . This determines, at fixed geometry, the attenuation of starlight by dust at any galaxy inclination θ, which measures the angle between the observer line-of-sight and the normal to the equatorial plane of a galaxy. At fixed τ_B⊥ , the integration over the solid angle of the attenuation curve in the T04 model yields the angle-average attenuation curve (see Sect. 2 of Chevallard et al. 2013). As in Charlot & Fall (2000), we couple the attenuation in the diffuse ISM described by the T04 model with a component describing the enhanced attenuation of newly born stars (t < 10 Myr) in their parent molecular clouds. Following Charlot & Fall (2000), we parametrize this enhanced attenuation by means of the fraction 1-μ of the total attenuation that arises from dust in stellar birth clouds, in the angle-average case.

The attenuation of the radiation emitted by a stellar generation of age t at inclination θ can therefore be written as (4)where the superscripts “BC” and “ISM” refer to attenuation in the birth clouds (assumed isotropic) and the diffuse ISM, respectively. In this expression the attenuation curve for the diffuse ISM is taken from the T04 thin stellar disk model, and following Wild et al. (2007, see also da Cunha et al. 2008, we compute the attenuation curve in the birth clouds as (5)where the V-band optical depth of the birth clouds is related to the angle-average optical depth of the diffuse ISM as (6)To compute the ratio of the infrared-to-ultraviolet luminosities IRX in this model, we take the IR luminosity to be equal to the fraction of all photons emitted in the range 912  Å ≤ λ ≤ 3  μm in any direction that are absorbed by dust (dust is almost transparent at λ > 3  μm). Assuming that photons at IR wavelengths emerge isotropically from a galaxy, we write the IR luminosity L_IR as (7)where is the unattenuated luminosity emitted by stars (assumed isotropic), Ω is the solid angle, and is the integral of Eq. (4) over the star formation history of the galaxy. We compute the monochromatic ultraviolet luminosity at the frequency ν corresponding to λ = 2300 Å as ℒ_NUV(θ) = νL_ν(θ), where L_ν(θ) is given by (8)and if the luminosity emitted by all stars at λ = 2300 Å in the direction θ, and c is the speed of light.

5.2. Model library

We use this model to compute a library of 20 000 SEDs of dusty star-forming galaxies, which we divide in bins of constant (NUV − r) and (r − K). We compute the mean sSFR and ⟨IRX⟩ in each bin and explore their distribution in the NUVrK diagram. After some experimentation, we find that a Gaussian distribution of τ_B⊥ centered at 7, with a dispersion of 3, truncated at the maximum value of available models τ_B⊥ = 8, and a Gaussian distribution of μ centered at 0.3, with a dispersion of 0.2, and truncated at μ = 0 and μ = 1, allow us to well reproduce the data, as shown in the top-left panel of Fig. 11. We note that a uniform distribution of galaxy inclinations would produce a large tail of highly attenuated galaxies, which is not observed in the data (at (r − K) > 2.5). Hence, in Fig. 11 we have adopted the observed distribution of axis ratios, converted to inclinations using the standard formula for an oblate spheroid (e.g., Guthrie 1992) (9)where q₀ is the intrinsic axis ratio of the ellipsoid representing the galaxy, which we fix to q₀ = 0.15.

A comparison between top-left and bottom-right panels of Fig. 11 shows that the models reproduce, at least qualitatively, the distribution of ⟨IRX⟩ in the color–color plane. The reddest near-IR colors, (r − K) > 1.4, correspond to galaxies seen at large inclination. This is consistent with Fig. B.2, which shows that the galaxies with reddest (r − K) colors have the largest measured ellipticities, i.e., they are more inclined. A large inclination makes the disk appear more opaque, since photons have to cross a larger section of the dust disk before they escape toward the observer.

The location and shape of the IRX stripes in the theoretical NUVrK diagram depend on several galaxy physical parameters, namely evolutionary stage, current SFR, dust content and distribution. The evolutionary stage and current SFR determine the relative amount of young and old stars in the galaxy, which controls the ratio of unattenuated ultraviolet to optical and near-IR luminosity of the galaxy. The global dust content, + , and the distribution of dust between ambient ISM and birth clouds affect the (NUV − r) and the (r − K) colors in different ways. For galaxies with a non-negligible fraction of young stars (i.e., log   (sSFR) ≳ − 9), the (NUV − r) color is mainly driven by the stellar birth clouds optical depth , and the (r − K) color by the optical depth of the diffuse ISM .

We test the effect of varying the optical depth of stellar birth clouds by computing the same library of 20 000 galaxy SEDs as described above, but fixing . Top-right panel of Fig. 11 shows that neglecting birth clouds attenuation prevents us from reproducing the reddest (NUV − r) color observed in the data. The stripes appear almost perpendicular to the (r − K) color driven primarily by the diffuse ISM. Also, this model predicts smaller values of ⟨IRX⟩ at a fixed position in the (NUV − r) vs. (r − K) diagram, since the UV photons emitted by young stars do not suffer enhanced attenuation by the dusty birth cloud environment, which would be re-emitted at IR wavelengths increasing the overall IR luminosity.

We also study the effect of neglecting the dependance of dust attenuation on galaxy inclination. To achieve this, we compute the same library of 20 000 SEDs as above, but we fix the attenuation curve to the angle-averaged curve . Bottom-left panel of Fig. 11 shows that this prevents us from reproducing the reddest (r − K) colors of the observed galaxies, which correspond to highly inclined objects (see top-left panel of Figs. 12 and B.2).

We have shown in Fig. 11 that to reproduce qualitatively the observed distribution of galaxies and the value and orientation of the ⟨IRX⟩ stripes in the (NUV − r) vs. (r − K) plane we need a prescription for dust attenuation which includes both a two-component medium (i.e., ISM + birth clouds) and the effect of galaxy inclination. We can now consider the “Full model” shown in the top-left panel of Fig. 11 and study how different dust properties vary in the (NUV − r) vs. (r − K) plane. Figure 12 shows the same library of SEDs as in the top-left panel of Fig. 11. As for Fig. 11, we divide the galaxy SEDs in bins of constant (NUV − r) and (r − K), and compute in each bin the mean value of galaxy inclination 1 − cosθ, V-band attenuation optical depth seen by stars younger [older] than 10⁷ yr [], and slope of the optical attenuation curve in the diffuse ISM .

The top-left panel of Fig. 12 shows that the galaxy inclination systematically increases as (r − K) increases, with a weaker dependence on the (NUV − r) color. This can be understood in terms of the attenuation optical depth in the diffuse ISM, which increases with increasing galaxy inclination, hence making (r − K) redder. The (NUV − r) color is also influenced by the variation of galaxy inclination, but to a much less extent since it also depends on the birth clouds attenuation optical depth (see Eq. (4)). The variation of the mean galaxy inclination shown in the top-left panel of Fig. 12 is also in qualitative agreement with Fig. B.2, which shows that the mean observed ellipticity of the galaxies in our sample increases from the bottom-left to the top-right side of the (NUV − r) vs. (r − K) plane.

The top-right panel of Fig. 12 shows the variation of the V-band attenuation optical depth suffered by stars younger than 10⁷ yr (i.e., , see Eq. (4)). At small (r − K) the stripes of constant are more parallel to the (NUV − r) color, while they become more and more inclined as (NUV − r) and (r − K) increase. This behavior can be understood in terms of the different fraction of light emitted by young stars attenuated by dust in the birth clouds and in the diffuse ISM. At small (i.e., blue) (r − K), young stars are mostly attenuated by the birth clouds component, which makes the (NUV − r) larger (i.e., redder, by decreasing NUV, at fixed r) without affecting (r − K). As (r − K) increases, the attenuation in the diffuse ISM increases, because galaxies are more inclined or have a larger dust content, and so the fraction of light emitted by young stars attenuated by this component raises too. As a consequence, the (NUV − r) color is determined by attenuation in both components, and the stripes of constant change orientation.

The variation of the V-band attenuation optical depth suffered by stars older than 10⁷ yr (i.e., ) shown in the bottom-left panel of Fig. 12 follow that of the (r − K) color, as indicated by the orientation of the stripes of constant almost perpendicular to (r  −  K). This is not surprising, since (r  −  K) traces stars older than 10⁷ yr, which are attenuated by dust in diffuse ISM. When moving from left to right on the (r  −  K) axis, increases from 0 to 2.5, indicating that the amount of attenuation suffered by stars in galaxies with very large (i.e., red) (r − K) is substantial.

The bottom-right panel of Fig. 12 shows the slope of the optical attenuation curve in the diffuse ISM , obtained by fitting a power law to the model attenuation curves in the range 0.4 ≤ λ ≤ 0.7  μm. The slope becomes smaller (i.e., the attenuation curve becomes flatter) when moving from the bottom-left to the top-right side of the diagram. This effect, as described in Chevallard et al. (2013, see their Fig 4), is a general prediction of radiative transfer models which consider disk galaxies with a mixed distribution of dust and stars. The variation of the slope of the attenuation curve in the (NUV − r) vs. (r − K) plane can account for the fact that the observed stripes of constant IRX are not perpendicular to the SMC and Calzetti attenuation vectors (see Fig. 4) and that the SED fitting predicts a systematic variation of the slope of the attenuation curve as a function of sSFR (see Fig. A.2).

With this analysis we have shown that we must account for the effect of geometry (i.e., of the spatial distribution of dust and stars) and galaxy inclination to reproduce the attenuation of starlight by the diffuse ISM, which mainly affects the (r − K) color. We have also shown the importance of accounting for the enhanced attenuation of newly born stars by their birth clouds to reproduce the reddest (NUV − r) colors and to match the observed values of ⟨IRX⟩. In the end, the good qualitative agreement between the “fully model” galaxies and the data (i.e., top-left and bottom-right panels of Fig. 11) confirms that the NUVrK diagram encodes valuable information about the global energy transfer between starlight and dust and galaxy inclination. We defer to a future work a more detailed and quantitative analysis of the data presented here, which would help us to better constrain the global amount, distribution and redshift evolution of the dust in star-forming galaxies.

6. Conclusion

We present a new method to compute the SFR of individual star-forming galaxies based on their location in the (NUV − r) versus (r − K) color–color diagram. We show that the NUVrK diagram provides an efficient way to separate quiescent and star-forming galaxies, an alternative to the UVJ diagram proposed by Williams et al. (2009). For the star-forming galaxies, the UV/optical luminosities in this diagram are highly sensitive to the shape of the dust attenuation laws. On the other hand, the infrared excess, IRX = L_IR/L_UV, as the net budget of the absorbed versus unabsorbed UV light, is weakly dependent on these effects. We combine the two dust diagnostics by analyzing the distribution of the mean infrared excess (i.e., ⟨IRX⟩) in the NUVrK diagram for a large sample of star-forming galaxies at redshift 0 ≤ z ≤ 1.3 selected from the COSMOS 24  μm and low-z GALEX-SDSS-SWIRE (Johnson et al. 2007a) samples. We observe the presence of stripes with constant ⟨IRX⟩, associated with a small dispersion around the mean, which allows us to describe ⟨IRX⟩ with a unique vector (i.e., NRK , a combination of (NUV − r) and (r − K) colors). We derive a simple relation between ⟨IRX⟩ and NRK, the norm of the vector NRK , and redshift, valid for star-forming galaxies with M ≥ 2 × 10⁹ M_⊙ and 0.05 − 0.1 ≤ z  ≤  1.3. This relation allows us to predict the IR luminosity of individual galaxies with an accuracy of ~0.2 dex (up to 0.27 dex for the local sample), which is better than the accuracy obtained with the SED fitting method based on 31 COSMOS pass-bands for the same galaxies.

We perform extensive comparisons of the L_IR and SFRs derived with the 24  μm, NRK and SED fitting methods. We find that the three methods provide consistent results for the vast majority of star-forming galaxies (~85–90%). The methods diverge for highly star-forming galaxies (SFR ≥ 100 M_⊙ yr^-1), which remain a negligible population at z ≤ 1.3, and for more evolved galaxies (sSFR ≲ 10^-10 yr^-1). For the latter, we describe a sSFR-, redshift-dependent limit below which the NRK method becomes unreliable and we also show that this population with inconsistent SFR estimates can be easily isolated and discarded in the NUVrK diagram.

By using the NRK method, we reconstruct the relationship between SFR and stellar mass for a K-selected sample of star-forming galaxies and find an excellent agreement with previous results over the entire redshift range.

Finally, we investigate the physical origin of the ⟨IRX⟩ stripes in the NUVrK diagram by appealing to a library of model SEDs based on the population synthesis code of Bruzual & Charlot (2003). We find that this library of models is able to qualitatively reproduce the location and shape of ⟨IRX⟩ stripes in the NUVrK diagram if we adopt a realistic prescription for dust attenuation (Chevallard et al. 2013). We show that to reproduce the observed stripes of ⟨IRX⟩ we must appeal to a two-component (i.e., birth clouds + diffuse ISM) dust model, which must account also for the effect on dust attenuation of galaxy inclination and geometry (i.e., the spatial distribution of dust and stars).

The method discussed in this work offers a simple alternative to assess the total SFR of star-forming galaxies in the absence of Far-IR observations or spectroscopic diagnostics. Because it directly predicts the infrared excess, no assumption on the dust attenuation curves is required to derive the SFR, in contrast to other methods such as the β-slope or SED fitting.

In a companion paper Le Floc’h et al. (in prep.), we extend our analysis toward lower stellar mass and higher redshifts based on the stacking technique of the Far-IR emission using the complete dataset available from Spitzer/MIPS at 24  μm to Herschel/SPIRE at 250, 350 and 500  μm.

---

Acknowledgments

We thank the anonymous referee for his useful comments. This work has been partially supported by the CNRS-INSU and the Progamme National de Cosmologie et Galaxies. J.C. acknowledges the support of the European Commission through the Marie Curie Initial Training Network ELIXIR under contract PITN-GA-2008-214227.

References

Online material

Appendix A: The SED fitting technique

The broad band SED fitting technique is a simple approach to infer the physical properties of a galaxy, such as stellar mass, SFR, amount of dust , age of stellar populations, by statistically comparing model and observed SED. The constraints on the physical parameters depend on the wavelength range spanned by the data and their quality. The COSMOS field, for which a wealth of multi-wavelength, high signal-to-noise ratio observations exist, is thus well suited for such modeling.

To derive the physical parameters, we adopt a library of SEDs based on the synthetic stellar population code from Bruzual & Charlot (2003, hererafter BC03). We describe the star formation history either with an exponentially declining function, with e-folding time 0.01 ≤ τ ≤ 15 Gyr, or with a constant. We adopt two metallicities, subsolar and solar (i.e., Z_⊙, 0.2 Z_⊙), and the IMF of Chabrier (2003), truncated at 0.1 and 100 M_⊙. Since the maximum redshift of the galaxies in our sample is z ~ 1.3, we force the age of galaxies to be larger than 100 Myr, computed from the onset of SF which, in our case, corresponds to the initial burst. We also constrain the age of the galaxies not to exceed the age of the universe at any redshift. We do not adopt rising star formation histories, since these were developed to improve the SED fitting at high redshift (z ≥ 2, see Maraston et al. 2010). The prescriptions for TP-AGB stars adopted in BC03 produces a lower near-IR luminosity for intermediate age stellar populations with respect to the prescriptions of Maraston (2005). This affects the galaxy mass-to-light ratio, and produces a difference in the stellar mass estimated with Maraston et al. (2010) of ~–0.15 dex. The dust attenuation curve encodes informations about the nature of the dust grains (sizes, chemical compositions) and the spatial distributions of dust and stars. Boquien et al. (2009) have shown the necessity of adopting a range of attenuation laws to reproduce the observed scatter in the IRX vs. β relation. In particular, they show that a gray (i.e., shallow) attenuation curve, such as the starburst curve of Calzetti et al. (2000), and a steeper curve, such the SMC extinction curve of Prevot et al. (1984), are required to span the observed distribution in the IRX versus β diagram of the starbursting and normal star-forming galaxies. Similar conclusions are reached by Ilbert et al. (2009), who find that a range of attenuation curves is required to improve the photometric redshift accuracy via the SED fitting method. For these reasons, we adopt three different attenuation curves: a starburst, SMC-like laws and a curve with a slope in between them (λ^-0.9). We then consider reddening excess in the range 0 ≤ E(B − V) ≤ 0.6, which allows us to explore the observed color distribution of our sample (see Fig. 3).

We use LePhare code (Arnouts et al. 1999; Ilbert et al. 2006) to compute the χ² for each observed galaxy and the entire model library, with all the photometric passbands from 0.15  μm to 4.5  μm. The physical parameters are derived by computing the median of the marginalized likelihood for each parameter and the errors corresponding to the 68% credible region.

In Fig. A.1, we compare the instantaneous SFR and sSFR derived with the SED fitting with the “total” SFR and sSFR derived from the observed IR and UV luminosities (i.e., Eq. (1)). The mean errorbars (based on 68% errors) for the SED parameters vary between 0.2 to 0.4, as shown by the gray region on the right side of the figures.

For the bulk of the 24  μm population the SFRs are in good agreement over ~3 order of magnitude, with a dispersion lower then a factor of two. The vast majority of sources with large discrepancies is located in the region occupied by passive galaxies (red triangles) or next to it, in the “intermediate” zone (green dots) as defined in Fig. 2. For those galaxies, the SED fitting predicts a low specific SFR (i.e., sSFR^SED ≤ 10^-10.5 yr^-1). The origin of this difference may have multiple causes: an inadequate description of dust attenuation may cause the SED fitting to reject highly attenuated models and prefers models with low or no dust content and a low specific SFR. Alternatively the adopted definition of the total SFR in Eq. (1) may over-estimate the SFR, since we neglect the contribution of old stars to the dust heating ((1 − η)L_IR). This contribution is often considered as 30% (η = 0.3) for star-forming galaxies (Inoue 2005) but it can be higher for the most evolved galaxies (e.g., Cortese & Hughes 2009). However, even an extreme value of η ~ 0.9 will not reconcile the estimated SFRs with the two methods. Another possibility is that the extrapolation of the 24  μm flux into the total IR luminosity could fail for galaxies with low specific SFR, if for exemple, a warmer dust temperature is associated to the same 24  μm flux with respect to galaxies in the star-forming main sequence, as recently reported by Skibba et al. (2011); Smith et al. (2012). It is beyond the scoop of this paper to address this issue, since we focus on the star-forming galaxies, but we find that ~7% of the entire 24  μm sample is affected by this mismatch in the SFR estimates.

In Fig. A.2, we show the relative contribution of the different attenuation laws, corresponding to the best-fit templates, as a function of the SFR and specific SFR. As mentioned above, the SMC-like extinction law is favored for galaxies with low SFR and/or low sSFR while the starburst law better fits the active/starbursting galaxies with high SFR and sSFR (sSFR ≥ − 9). Our results agree with Wuyts et al. (2011) with a transition for a steeper attenuation law at SFR ≤ 20 M_⊙ yr^-1. The most active galaxies are consistent with a mixed distribution of the dust and star resulting in the gray attenuation law (Calzetti et al. 2000), while “normal” star-forming galaxies prefer the SMC-like attenuation consistent with a simple dust screen model. See Sect. 5 and Chevallard et al. (2013) for a purely geometric origin of different attenuation laws.

In conclusion, our analysis shows that SFRs estimated via the SED fitting technique are in good agreement with those derived from the UV+IR contribution. This validates our method to derive the L_IR from the 24  μm flux and the use of Eq. (1) as a good measure of the SFR. we have noted that Eq. (1)may become inadequate to describe the SFR in more evolved galaxies, possibly because of the presence of a larger population of old stars or the inadequate conversion of L_24 μm to L_IR. While this problem can not be easily solved, in this paper we show that we can isolate in the NUVrK diagram the region occupied by galaxies for which we obtain inconsistent SFR estimates with the different methods.

Appendix B: Separation between passive and active

As the UVJ diagram proposed by Williams et al. (2009), in Sect. 2.3 and in Fig. 2, we define a criterion based on the NUVrK diagram to separate the passive and star-forming galaxies. To test the validity of the above color criterion, we show the distribution in the NUVrK diagram of galaxies with morphological information provided by the Zurich estimator of structural type (ZEST) catalog (Scarlata et al. 2007). The ZEST classification uses three morphological classes: Early type, Disk and Irregular, with subclasses describing the degree of “irregularity” in the early-type class (i.e., 0 for regular, 1 for irregular), and the contribution of the bulge for disk galaxies (i.e., from 0 for bulge dominated galaxies, to 3 for pure, bulge-less disks). We also consider the ellipticity class for the galaxies classified as disks. This traces the galaxy inclination with an ellipticity of 0 corresponding to a face-on galaxy and up to 3 for an edge-on galaxy. The ZEST catalog includes galaxies down to I_AB ≤ 24 and we restrict the sample to galaxies with good flags. The distributions of the different morphological classes in the NUVrK diagram are shown in Fig. B.1 for the whole ZEST sample (as density contours) and for the 24  μm selected subsample (small dots). We detect a clear evolutionary sequence in this diagram, with the Irregular (IRR) and Spiral-disk dominated (Sp-3, and Sp-2) galaxies showing only blue colors, typical of active star-forming galaxies. The Spirals with a growing contribution of bulge (Sp-1 and Sp-0) and the early-types (ELL) show an increasing fraction of their population to lie in the passive region (i.e., top-left side) of the diagram. The 24  μm sample tends to lie in the blue plums of the early type (ELL) and bulge dominated spiral (Sp-0) samples. The blue plum in the early type class could be due to some residual of star formation activities (the plum is present in the two subclasses based on the regularity criterion: ELL-0 and ELL-1).

The Spiral disk-dominated population (Sp-3, and Sp-2) extend to relatively red colors in (NUV − r) and (r − K) (top-right part of the diagram), where high values of IRX are also observed. As discussed by Patel et al. (2011), in the UVJ diagram, the disk inclination can be indeed responsible for this extreme reddening. This is indeed supported by Fig. B.2, where we show the mean values of the ellipticity parameter for the disk population, in different redshift bins (see also the discussion in Sect. 5).

Appendix C: The average spectral energy distributions along the vector NRK

The evolution of ⟨IRX⟩ with NRK vector should be reflected in the shape of the galaxy SEDs when moving from low to high value of NRK. In Fig. C.1, we have reconstructed the averaged rest-frame SEDs in different bins of NRK, redshift and stellar mass, by using the 31 broad and medium bands available in the COSMOS catalog. They are reconstructed in logarithmic wavelength bins (Δlog   (λ) = 0.15). Before computing the average and the dispersion in each wavelength bin, the low resolution spectrum of all the galaxies was first normalized at λ ~ 0.5 μm (dashed lines). In each panel, we show the number of galaxies used to reconstruct the averaged SEDs and as an indicative value of the IRX, we show the ratio log   (F_IR/F_NUV), where the F_IR is the redshift corrected 24 μm flux of the averaged SED (F_{24 μ/(1 + z)}) and F_NUV the rest-frame flux at 0.23 μm (F_0.23 μ). As a reference, we overplot the template of a young star-forming galaxy (SB2⁴ from Ilbert et al. 2009). For the first bin 0.5 ≤ NRK ≤ 1, corresponding to the bluest population in our sample, no attenuation is applied to the template. For the three other bins of NRK, we increase the amount of dust reddening to E(B − V) = 0.2,0.4,0.6, assuming a starburst attenuation law (solid black lines) and a SMC extinction law (gray dashed lines).

A strong evolution of the SED shape is observed with increasing NRK. We can quantitatively reproduce this evolution by increasing the reddening excess applied to our SB2 star-forming template. On the other hand, redshift and stellar mass appear to play a secondary role in shaping the SED properties. This provides an additional evidence that NRK is sensitive to the global budget between UV and Far-IR emission in galaxies, regardless of other properties. It also supports our choice of neglecting the dependence on the stellar mass in Eq. (3). However, in Le Floc’h et al. (in prep.) we extend our analysis to a complete mass selected sample using the stacking technique at 24, 220, 350 and 500  μm and we show that the ⟨IRX⟩ vs. NRK relation should include a corrective term based on stellar mass when considering galaxies with stellar masses M_⋆ ≤ 10^9.3 M_⊙.

Appendix D: Dependence with sSFR

As observed in Fig. 8, the difference between the IR luminosity derived from the 24  μm observations and the NRK method varies with the specific SFR. In particular, we have shown that the NRK method over-estimates L_IR for galaxies with low sSFR. In Fig. D.1 we compare, in different redshift bins, the predicted specific SFR (sSFR^NRK) with the reference sSFR (derived with the 24  μm, sSFR_tot), with both sSFR estimated with Eq. (1). The sSFR^NRK tends to saturate while the true sSFR keeps declining toward galaxies with low star formation activity, and this deviation varies with redshift. We propose an analytical fit which allows us to reproduce this deviation over the entire redshift range: (D.1)\vspace*{1mm}where, y₀ = log  (sSFR^NRK) + 9.5, a(z) = −1.4 + 0.8 z and b(z) = −0.1 (1 − z) and z is the redshift. By mean of this equation, we can define a sSFR threshold below which the NRK method over-estimates the SFR by a given factor. By chosing a factor 2, the redshift evolution of this threshold can be simply described by sSFR^NRK_lim(z) = −10.6 + 0.8 z. These thresholds are shown as shaded yellow regions in Fig. D.1.

All Figures