© ESO, 2014

1. Introduction

The star formation rate (SFR) is a key parameter to understand galaxy formation and evolution. If the SFR plays such a central role, this is because of its fundamental relation to the stellar mass and the gas reservoir on which our understanding of galaxy formation and evolution is based. Yet, despite its considerable importance, measuring star formation accurately remains a challenging task (see Kennicutt & Evans 2012, for a recent review).

One of the key assumptions to measure the SFR of a galaxy is its star formation history (SFH). Most classical estimators are based on the assumption of a constant SFR over a period of 100 Myr (e.g. Kennicutt 1998). If this assumption seems reasonable for low redshift spiral galaxies evolving secularly, it is unlikely to hold true for interacting systems or at higher redshift where the SFR necessarily varies on timescales that can be similar or shorter than 100 Myr. With the increasing availability of observations spanning a broad range of wavelengths, one possibility for waiving the assumption on the SFH is to carry out multi-wavelength spectral energy distribution (SED) modelling. The SFH can then be left as a free parameter to compute the SFR. Unfortunately, the SFH is degenerate with other parameters such as the attenuation. For instance, a galaxy can exhibit a red ultraviolet (UV) slope because it is dusty and actively forming stars, or because it is not forming stars any more. If multi-wavelength data partly alleviate this problem when UV and infrared (IR) data are simultaneously available, the SFH still remains poorly constrained.

Even if we cannot uncover the SFH of galaxies in detail, it is nevertheless possible to investigate the impact of short- and long-term variations of the SFH on the measure of the SFR using numerical simulations of galaxies (e.g. Wuyts et al. 2009; Wilkins et al. 2012; Lanz et al. 2014; Simha et al. 2014; Domínguez et al. 2014). However, simultaneously taking both the attenuation and a variable SFH into account makes it difficult to isolate the exact impact of the latter on the accuracy of classical SFR estimators. Another potential issue that can be encountered when exploiting numerical simulations is the time resolution of reconstructed SFH. If it is not significantly smaller than the typical time sensitivity of an SFR estimator, then short-term variations are in effect smoothed out. The impact of these variations cannot be explored, and thus remains unknown.

The aim of this article is to explore the impact of long-term and short-term variability of the SFH on the measurement of the SFR using classical estimators, isolated from any other perturbing effect. In Sect. 2, we lay out our method, including how we combine state-of-the-art hydrodynamical simulations of galaxies with an SED modelling code to generate synthetic observations from which we estimate the SFR. We present our results in Sect. 3, and discuss the best course of action to measure the SFR while taking its variability into account in Sect. 4. We conclude in Sect. 5.

2. Method

To understand the impact of the variability of the SFH on the measurement of the SFR, our main guideline is to compare the estimated SFR from synthetic observations of simulated galaxies to the true SFR we know from the simulations.

In this section, we first present the set of simulations of high redshift galaxies that provides us with a sample of SFH with a time resolution of only 1 Myr. Then we show how we exploit these SFH to create simulated observations. Finally, we present the set of SFR estimators considered in this study and how they are applied to estimate the SFR from simulated observations.

2.1. The MIRAGE simulations sample

The use of SFH with plausible variations over all timescales, from rapid stochastic variations to slow variations over long periods, is especially important if we want to provide meaningful insight into the impact of the variability of the SFH on the measure of the SFR. This means that we need SFH generated by means of hydrodynamical simulations with a very detailed modelling of baryonic physics and, in particular, of feedback that can have a dramatic effect in governing the SFR in galaxies. Following these requirements, we rely on the MIRAGE simulations (Merging and Isolated high redshift Adaptive mesh refinement Galaxies, Perret et al. 2014). It comprises a set of 20 idealised merger simulations sampling four initial disk orientations and five total baryonic masses. The whole merger sample is built using three distinct idealised disk models G1, G2, and G3, which are also simulated in isolation with stellar masses of log M_⋆ = 9.8, 10.2, and 10.6 respectively. These stellar masses were selected to sample the stellar mass histogramme of a representative spectroscopic sample of 1 <z< 2 galaxies, namely the MASSIV sample (Contini et al. 2012). Each model is composed of a stellar disk, a gaseous disk, a stellar bulge, and a dark matter halo. The two disk components are built using an exponential density profile, while the stellar bulge and the dark matter halo follow a common Hernquist (1990) profile. The initial gas fraction of the three disk models is f_g = 0.65, which is in agreement with the most gas-rich 1 <z< 2 population observed (Daddi et al. 2010a). These simulations were carried out using the adaptive mesh refinement code RAMSES (Teyssier 2002), allowing us to reach a physical resolution of 7 pc in the most refined cells. The modelled gas physics in the simulations includes metal lines cooling, star formation using a local Schmidt law, metallicity advection as well as a recent physically-motivated implementation of stellar feedback (Renaud et al. 2013) that encompasses OB-type stars radiative pressure, photo-ionization heating, and supernovae. The SFH of each simulation has been obtained by computing a mass-weighted histogramme of the age of the stars particles using a bin of 1 Myr over about 800 Myr.

Because of the high gas fraction of the galaxy models, Jeans instabilities arise in the first dynamical times of the simulations and trigger very high SFR, more than 2σ away from the normal star formation regime of the simulations. This transitory relaxation phase that occurs within the first 250 Myr is excluded from our analysis to avoid spurious effects due to its non-representativeness.

While the physical properties of the three disk models of the MIRAGE sample are defined to be comparable to the properties of galaxies observed in the redshift range 1 <z< 2, the star formation implementation is not specifically designed to perfectly match the normal star formation efficiency (SFE) of the observations. Indeed, MIRAGE simulations display a SFE slightly under the Daddi et al. (2010b) relation, within 1σ. Nevertheless, this feature is not really relevant to this analysis since we are interested in the relative SFR fluctuations more than the absolute value. Moreover, the merger simulations do not display any SFE enhancement that is interpreted as a saturation mechanism of the star formation in Perret et al. (2014). The remarkable homogeneity of the observed specific SFR in high redshift galaxies (Elbaz et al. 2007, 2011; Nordon et al. 2012) suggests that the starburst regime might be infrequent, with most galaxies being on the main sequence at any given time. The MIRAGE galaxies always lie on the main sequence, which is consistent and suggests that they do not represent rare objects but rather common galaxies.

One of the main findings of Perret et al. (2014) is the variation of the SFR on short timescales, with a dispersion σ_SFR ≃ 0.3 × ⟨SFR⟩. They attribute this to the presence of multiple long-lived giant star-forming clumps wandering in a turbulent ISM. If they are indeed widespread, these rapid variations of the SFR show the importance of fully understanding their impact on classical SFR estimators.

There are two caveats that should nevertheless be noted. First, because these simulations are computationally very demanding, galaxies were followed over less than 1 Gyr. This means that we do not have access to SFH on longer timescales, which prevents us from modelling the contribution of an old stellar population. Cosmological simulations naturally provide SFH over 13 Gyr if they are run until z = 0, but by definition such simulations (e.g. Pontzen & Governato 2012; Dubois et al. 2014; Vogelsberger et al. 2014) cannot reach the numerical resolution achieved in idealised simulations (Renaud et al. 2013). A numerical resolution that allows us to resolve the vertical structure of the disk combined with a gas treatment that allows for cooling below 10⁴ K is essential to understanding the nature of star formation. Hydrodynamical simulations of a cosmological volume are able to display a cosmic SFH coherent with observations (Hopkins & Beacom 2006), but at the scale of individual galaxies the SFH is smoothed because of the lack of both spatial and temporal resolution. Using idealised simulations in this analysis, we trade off a realistic long-term SFH for more realistic SFR fluctuations on shorter timescales. The second point is that these simulations do not aim to be representative of all main-sequence galaxies at the aforementioned redshifts. In reality, the global population of high-redshift galaxies is likely to exhibit a broader diversity of SFH than can be explored here. However, the sample we have adopted is nevertheless very useful to exemplify some of the issues with the measurement of the SFR, even if the exact amplitude of these issues may, in some circumstances, be dependent on the considered SFH.

Fig. 1 Left: SFH of the MIRAGE G2_G2_90_90_90 merger simulation (Perret et al. 2014). The time resolution is 1 Myr. The first 250 Myr are not shown as the initial relaxation is handled (Sect. 2.1). The coloured arrows at 400, 500, and 600 Myr indicate at which time we compute the spectra shown in the right-hand panel from the Lyman break to the U band. These spectra are created with the cigale SED modelling code (Boquien et al., in prep.; Burgarella et al., in prep.). In the present case we assume a Chabrier (2003) IMF and a metallicity Z = 0.008. Nebular emission is also included. It takes into account recombination lines and the nebular continuum from free-free, free-bound, and two-photon processes. The continuum contributes most in the red part of the shown spectra. We have computed the spectra for all 23 MIRAGE simulations (three isolated galaxies and 20 interacting galaxies) every 1 Myr.

2.2. Simulating observations

For each of the 23 MIRAGE simulations, we compute the attenuation-free spectra every 1 Myr from the SFH, foregoing all spatial information. To do so, we use the cigale SED modelling code (Boquien et al., in prep.; Burgarella et al., in prep.). We adopt the stellar population models of Bruzual & Charlot (2003) with a Chabrier (2003) initial mass function (IMF) between 0.1 M_⊙ and 100 M_⊙ and a metallicity Z = 0.008, consistent with the sub-solar metallicity of galaxies at the considered redshifts (e.g. Yuan et al. 2013). We will discuss the impact of the metallicity on our results in Sect. 4.2. The emission from H ii regions includes recombination lines and the nebular continuum via free-free, free-bound, and two-photon mechanisms.

The luminosity in individual bands is computed by integrating the spectrum in the corresponding passbands except for the number of Lyman continuum photons, which is an output parameter of the model, and for the total infrared (TIR) luminosity, which is computed by integrating over the entire spectrum, based on the assumption that the galaxy is entirely buried, following Kennicutt (1998).

An example of a MIRAGE merger SFH and the corresponding spectra from the Lyman break to the U band is shown in Fig. 1 at 400, 500, and 600 Myr. While the details of the SFH differ somewhat from one simulation to another, for instance in terms of normalisation, they are qualitatively similar. There is no long-term trend after the initial star formation episode and peak-to-peak variations are of the order of a few. Overall, the SFH shown in Fig. 1 is typical of the MIRAGE SFH. In terms of relative scatter, we find σ_SFR ≃ 0.36 × ⟨SFR⟩ versus 0.34 ± 0.06 for the entire sample starting from 250 Myr. All the SFH are shown in Fig. 7 in Perret et al. (2014).

2.3. Star formation rate estimators

2.3.1. Converting a luminosity to an SFR

The computation of the SFR from the luminosity in a given band relies on several assumptions such as a proper correction for partial obscuration by dust, an IMF, a metallicity, and as we have mentioned in Sect. 1, an SFH. Because we want to completely isolate effects due to the variability of the SFH from other effects, the simulated observations and the SFR estimators must be generated with the exact same set of parameters and models. We therefore refrain from adopting SFR estimators from the literature, rather we devise a fully consistent set of estimators. As we will see in Sect. 2.3.2, different SFR estimators probe different timescales. A constant SFR over 100 Myr or 1 Gyr is often assumed (e.g. Hao et al. 2011). We choose the former as our baseline because it typically corresponds to the age at which a constant SFR produces 90% of the steady state luminosity (Kennicutt & Evans 2012; Boissier 2013) and it is widely used in the literature.

To summarise, we compute the SFR estimators building a synthetic model with the following assumptions: 1) constant SFR over 100 Myr; 2) a Chabrier (2003) IMF between 0.1 M_⊙ and 100 M_⊙; and 3) a metallicity Z = 0.008. The calibration coefficients are presented in Appendix A. Conceptually similar SFR estimators spanning a broad range in terms of age and metallicity have been presented in Otí-Floranes & Mas-Hesse (2010).

In the context of this study, we consider five of the most popular SFR estimators: the Galaxy Evolution Explorer (GALEX) far-ultraviolet (FUV), GALEX near-ultraviolet (NUV), and Bessel U bands, which trace the photospheric emission from massive stars, the Lyman continuum, which is responsible for ionising the hydrogen in star-forming regions and can be traced through recombination lines or free-free emission, and the total IR emission from dust in a completely buried galaxy.

2.3.2. A short note on timescales probed by SFR estimators

Different SFR estimators are generally not sensitive to the same stellar populations. For instance, recombination lines are due to massive stars that live of the order of 10 Myr. However, stars that dominate the UV domain live up to a few hundred Myr. Yet, these timescales should rather be seen as an indication of the time required to reach a luminosity steady state for a constant SFR. To quantify the typical age of stars dominating the luminosity in a band for a given SFH, we introduce luminosity weighted ages: $\frac{{\sum }_{\mathit{t}}\mathit{t}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\sum }_{\mathit{t}}\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}}$, with t the age and L(t) the luminosity of a simple stellar population of age t in a given band. In practice, we sum over t, with each t representing a bin of 1 Myr. In Fig. 2, we present the luminosity weighted ages versus the wavelength for a constant SFR over 100 Myr and over 1 Gyr.

Fig. 2 Luminosity weighted ages versus the wavelength for a constant SFR over 100 Myr (solid lines) and 1 Gyr (dashed lines). The magenta and red lines indicate the luminosity-weighted age for the Lyman continuum and the TIR for a completely buried galaxy, respectively. We see that the typical timescale probed is not only strongly dependent on the estimator and on the wavelength, but also on the SFH, with longer wavelengths and the TIR being the most sensitive to the duration of star formation episode in this simple scenario.

We find that there is a strong dependence of the luminosity weighted ages on the length of the star formation episodes. Unsurprisingly, an episode of 1 Gyr systematically leads to larger ages than in the case of 100 Myr. This also depends on wavelength: longer wavelength bands yield larger ages and on average the TIR is more sensitive to older stars than the UV or the Lyman continuum. This shows the importance of older stellar populations for dust heating when considering timescales longer than 100 Myr, in line with well-known results (e.g. Lonsdale Persson & Helou 1987; Sauvage & Thuan 1992; Bendo et al. 2010, 2012; Boquien et al. 2011; Crocker et al. 2013).

Even for a naïve SFH, we see the complex relation there is between 1) the length of the star formation episode; 2) which stars contribute most to a given tracer; and 3) the strong differences from one tracer to another in terms of age sensitivity. More realistic SFH certainly render these relations even more complex. This shows that the presence of long-lived stars also has an impact on SFR estimates. We will examine these issues in detail over the following sections.

3. Results

3.1. Comparison of estimated and true star formation rates

For each of the 23 simulations, we have estimated the SFR using the conversion factors devised in Sect. 2.3.1 and listed in Appendix A, which assume a metallicity Z = 0.008, a Chabrier (2003) IMF between 0.1 M_⊙ and 100 M_⊙, and a constant SFR over 100 Myr. We compare the estimated SFR with the true SFR for one of the simulations in Fig. 3.

Fig. 3 Left: true SFR (black), SFR(Ly) (magenta), SFR(FUV) (blue), SFR(NUV) (cyan), SFR(U) (green), and SFR(TIR) (red) for the MIRAGE G2_G2_90_90_90 simulation (Perret et al. 2014). Note that the true SFR and SFR(Ly) are nearly blended. Right: zoom between 350 Myr and 450 Myr. We see that SFR(Ly) follows the variations of the true SFR extremely well even for high frequency variations. SFR estimators at longer wavelengths are not able to capture these variations though, leading to important discrepancies in the estimates of the instantaneous SFR. Beyond the capture of short-term variations, we see that there is also a systematic offset that is visible in the FUV and NUV bands but which is especially prominent in the U band. For the present simulation, we obtain ⟨[SFR(FUV) − SFR(true)] /SFR(true)⟩ = 0.27 ± 0.30 between 250 Myr and 790 Myr with a peak of ~1.5. For the U band, we obtain ⟨[SFR(U) − SFR(true)] /SFR(true)⟩ = 0.70 ± 0.46 with a peak of ~2.8. This is representative of the entire MIRAGE sample as can be seen in Fig. 4.

We find that SFR(Ly) follows very closely the true SFR and is able to capture all but the highest frequency variations for which it suffers from a short time delay. Conversely, SFR(FUV), SFR(NUV), SFR(U), and SFR(TIR) do not properly capture these variations on timescales of no more than a few tens Myr. A related and important question is whether the estimated SFR are accurate on average on long timescales for which the short term variations could be averaged out. We find that in reality the estimated SFR systematically overestimate the true SFR. If the overestimate is very small for SFR(Ly) it is clearly visible for SFR(FUV) and SFR(NUV), and it is particularly strong for SFR(U): we find ⟨[SFR(FUV) − SFR(true)] /SFR(true)⟩ = 0.27 ± 0.30 and ⟨[SFR(U) − SFR(true)] /SFR(true)⟩ = 0.70 ± 0.46 between 250 Myr and 790 Myr (we forego the first 250 Myr to avoid being affected by the highly transitory nature of the onset of star formation). To examine this issue further, in Fig. 4 we compare the ratio of the difference between the estimated and true SFR to the true SFR versus time for the entire sample. We find that this excess can also be retrieved when considering the entire sample:

• ⟨[SFR(Ly) − SFR(true)] /SFR(true)⟩ = 0.02 ± 0.11,

• ⟨[SFR(FUV) − SFR(true)] /SFR(true)⟩ = 0.24 ± 0.26,

• ⟨[SFR(NUV) − SFR(true)] /SFR(true)⟩ = 0.33 ± 0.29,

• ⟨[SFR(U) − SFR(true)] /SFR(true)⟩ = 0.65 ± 0.42, and

• ⟨[SFR(TIR) − SFR(true)] /SFR(true)⟩ = 0.38 ± 0.29.

This means that except for the Lyman continuum, classical estimators calibrated assuming a constant SFR over 100 Myr all clearly overestimate the SFR. We will discuss the origin of this effect in Sect. 4.1 and its impact in Sect. 4.3.

3.2. Timescales

We have seen in Sect. 2.3.2 that the typical timescales probed by different estimators are strongly dependent on the duration of the star formation episode. However, we examined the case of a constant SFR over a long period. It is unclear how these timescales compare with those for more realistic SFH. In Fig. 5, we have computed the luminosity weighted ages versus the simulation age and versus the wavelength. There is generally a progressive increase of the luminosity weighted ages with time. However, after 400 Myr, FUV progressively reaches a steady state. When considering longer wavelengths, the weighted ages keep increasing because of the continuing accumulation of long-lived stars. Short-term variations of the SFR seem to have little overall impact. At 790 Myr, the Lyman continuum has a luminosity weighted age of ~3 Myr, the FUV band ~40–60 Myr, the NUV band ~60–90 Myr, and the TIR and the U band tracers more than 100 Myr. Finally, when compared to Fig. 2, we find that the dependence of the ages on the wavelength is remarkably similar to the case of a constant SFR over 1 Gyr.

4. Discussion

4.1. Why the SFR is overestimated

We have found that except for the Lyman continuum, other SFR estimators systematically overestimate the true SFR over long periods, from 24% on average for the FUV band up to 65% for the U band. We examine here the fundamental reasons for this discrepancy and under which conditions standard estimators could provide reasonably unbiased results.

Intuitively, SFR estimators sensitive to short timescales such as SFR(Ly) provide us with accurate estimates because they tightly follow the variations of the true SFR both when it increases and when it decreases, albeit with a slight delay. However, as other tracers are sensitive to longer timescales we are facing the combined impacts of short- and long-term variations of the SFR. To investigate these issues, we have decomposed the MIRAGE SFH into low and high frequency components using a Fast Fourier Transform. While the choice of the threshold in terms of variation period between these components is somewhat arbitrary, we have chosen a separation period of 40 Myr, as it is similar to the FUV age-weighted luminosities. That is to say that all variations faster (resp. slower) than a timescale of 40 Myr are assigned to the high (resp. low) frequency component. From each of these two SFH, we have recomputed the estimated SFR independently. We show the Fourier-filtered SFH along with the corresponding recomputed estimated SFR for one of the MIRAGE simulations in Fig. 6. Note that the estimated SFR are not Fourier-filtered, they are only computed from Fourier-filtered SFH.

Fig. 4 Ratio of difference between the estimated and true SFR to the true SFR versus time for the entire MIRAGE sample. The solid lines indicate the mean value across all simulations whereas the shaded areas of the same colour indicate the standard deviation. The dashed lines indicate the mean value between 250 Myr and 790 Myr. The colour coding is the same as for Fig. 3. We see that except for the Lyman continuum, the SFR is systematically overestimated, from 24% on average for the FUV band to 65% on average for the U band. Peak relative differences reach ~0.7 for the Lyman continuum, ~1.9, ~2.1, and ~3.2 for the FUV, NUV, and U bands, and ~2.0 for the TIR.

First, we only consider the high frequency component. We see that only SFR(Ly) can follow the rapid variations of the true SFR. But even then, there is always a slight delay. Other tracers cannot follow rapid variations accurately but no systematic bias is visible. If we look at the entire sample between 250 Myr and 790 Myr, we find:

• ⟨[SFR_HF(Ly) − SFR_HF(true)] /SFR(true)⟩ = 0.01 ± 0.11,

• ⟨[SFR_HF(FUV) − SFR_HF(true)] /SFR(true)⟩ = 0.03 ± 0.17,

• ⟨[SFR_HF(NUV) − SFR_HF(true)] /SFR(true)⟩ = 0.03 ± 0.18,

• ⟨[SFR_HF(U) − SFR_HF(true)] /SFR(true)⟩ = 0.03 ± 0.17, and

• ⟨[SFR_HF(TIR) − SFR_HF(true)] /SFR(true)⟩ = 0.03 ± 0.16.

On average, the bias due to rapid variations of the SFR is limited to 3% and even just 1% for SFR(Ly). However, rapid variations of the SFR increase the scatter on the instantaneous SFR as even the SFR(Ly) cannot follow the variations closely enough. If we compare these values with those obtained in Sect. 3.1, it appears that rapid variations are responsible for nearly all of the scatter for SFR(Ly). For other tracers, however, this is not true as the scatter due to rapid variations of the SFR alone cannot explain all of the total scatter. This suggests that low frequency variations play an important role.

Fig. 5 Left: luminosity weighted ages versus the simulation time. The solid lines indicate the mean value of for the 23 simulations and the shaded areas of the corresponding colours indicate the standard deviation. The colour scheme is the same as that of Fig. 3. Right: luminosity weighted ages versus wavelength (black) for a simulation age of 790 Myr. The values for the Lyman continuum and the TIR at 790 Myr are also indicated. For comparison, the corresponding values for a constant SFR over 1 Gyr shown in Fig. 2 is indicated in dashed lines.

Fig. 6 Left: Fourier-filtered SFH of the MIRAGE G2_G2_90_90_90 simulation (black lines) along with the estimated SFR recomputed from each of the two filtered SFH. The top (resp. bottom) lines represent the low (resp. high) frequency component. The colours follow the same scheme as that of Fig. 3. Note that the true SFR and SFR(Ly) are nearly blended. Right: zoom between 350 Myr and 450 Myr. We find that if there are important discrepancies in the SFR, on average at high frequency the true SFR and the estimated SFR give similar results, both close to 0 Msun/yr. The behaviour at low frequency is markedly different. We see that the estimated SFR show both a systematic overestimate and a delay in time from the extrema of the true SFR. The smallest delay occurs for SFR(Ly). Conversely, SFR(FUV), SFR(NUV), SFR(U), and SFR(TIR) all show approximately similar delays of a few Myr.

If we now consider only the low frequency component, we can see two main features. First, in the right panel of Fig. 6 we see a clear offset between the extrema of the true SFR and that of the estimated SFR. This delay is as expected minimal for SFR(Ly), around 1 Myr. For other estimators it is typically delayed by a few Myr. This delay appears remarkably independent of the wavelength from SFR(FUV) to SFR(TIR). The second obvious feature is that except for SFR(Ly) there is an excess in the estimated SFR. It appears that tracers having longer timescales have a stronger excess. If we consider the entire sample, we find:

• ⟨[SFR_LF(Ly) − SFR_LF(true)] /SFR(true)⟩ = 0.00 ± 0.03,

• ⟨[SFR_LF(FUV) − SFR_LF(true)] /SFR(true)⟩ = 0.21 ± 0.17,

• ⟨[SFR_LF(NUV) − SFR_LF(true)] /SFR(true)⟩ = 0.30 ± 0.20,

• ⟨[SFR_LF(U) − SFR_LF(true)] /SFR(true)⟩ = 0.62 ± 0.33, and

• ⟨[SFR_LF(TIR) − SFR_LF(true)] /SFR(true)⟩ = 0.35 ± 0.20.

This suggests there is an impact due to the long-term variations of the SFR and/or a contamination from stars living longer than the calibration timescale, here 100 Myr. To investigate the latter possibility, we have computed the fraction of the luminosity contributed by such long-lived stars for the entire sample in Fig 7. We find that except for the Lyman continuum, long-lived stars can account for a sizeable fraction of the luminosity in star formation tracing bands. The typical contribution ranges from slightly more than ~10% for the FUV band to ~35% for the U band. While such a contribution is not sufficient to explain the full extent of the bias shown in Sect. 3.1, it already explains at least half of the observed bias. We show the estimated SFR after correction for the presence of long-lived stars for one MIRAGE simulation in Fig. 8. We find that the SFR are in better agreement once they are corrected for the contamination by stars living longer than 100 Myr. The change is the most spectacular for the U band. To examine this issue further, in Fig. 9 we compare the ratio of the difference of contamination-corrected estimated SFR with the true SFR to the true SFR for the entire sample. We find that the bias is strongly reduced:

• ⟨[SFR(Ly) − SFR(true)] /SFR(true)⟩ = 0.02 ± 0.11,

• ⟨[SFR(FUV) − SFR(true)] /SFR(true)⟩ = 0.07 ± 0.20,

• ⟨[SFR(NUV) − SFR(true)] /SFR(true)⟩ = 0.08 ± 0.21,

• ⟨[SFR(U) − SFR(true)] /SFR(true)⟩ = 0.09 ± 0.22, and

• ⟨[SFR(TIR) − SFR(true)] /SFR(true)⟩ = 0.05 ± 0.18.

Of the remaining bias, a third to a half of it is due to rapid variations of the SFR. The rest of it is probably due to slower deviations from the assumption of a constant SFR.

Fig. 7 Fraction of the luminosity in star formation tracing bands contributed by stars living over 100 Myr versus simulation time. The colour scheme is the same as that of Fig. 3. The contribution of long-lived stars is negligible for the Lyman continuum but reaches upward 10% in the FUV, a small but non-negligible fraction. The contamination is stronger for the NUV and TIR tracers, and it is maximum for the U band with a typical value of ~35%.

In summary, except for the Lyman continuum, standard SFR estimators calibrated over 100 Myr provide us with biased results: from ~25% for SFR(FUV) to ~65% for SFR(U). This is mainly because of the accumulation of stars living over 100 Myr. After correction for long-lived stars, the bias is smaller than 10% on average. The scatter is due to the combination of short-term and long-term variations of the SFR. These results are in line with the recent findings of Johnson et al. (2013). Reconstructing SFH from colour-magnitude diagrammes of a sample of nearby dwarf galaxies, they found that long-lived stars contribute from less than 5% to 100% for the FUV luminosity.

4.2. Effect of the metallicity

We have carried out the present study for Z = 0.008; however, the metallicity has a direct effect on stellar evolutionary tracks and on stellar atmospheres. To understand the exact impact of the metallicity on our results, we have recomputed the estimated SFR at different metallicities. All estimators have been systematically recalibrated for each metallicity to ensure the study remains self-consistent. We find that on average, as the metallicity decreases the bias increases. For the FUV (U) band, the relative bias ranges from 0.19 ± 0.23 (0.57 ± 0.38) for Z = 0.02 to 0.37 ± 0.32 (0.87 ± 0.52) for Z = 0.0001. Because low metallicity stellar populations tend to be bluer than more metal-rich populations of the same age, at low metallicity there is probably an increased accumulation of long-lived stars contributing to star-formation tracing bands. However, for such extreme metallicities, the time from the onset of star formation is necessarily much smaller than at 1 <z< 2 so the time available to accumulate long-lived stars is shorter and the metallicity rapidly increases, limiting the impact of the contribution of such stars.

4.3. Impact on practical measurements of the SFR

In the ideal case we have considered, we find that tracers based on the Lyman continuum, such as free-free emission or recombination lines estimate the SFR very accurately, within 2%. The reason is that tracers based on the Lyman continuum provide the instantaneous SFR, or the rate at which massive stars are being formed within a time lapse of a few Myr, while other tracers are giving the SFR averaged in a complex way over the typical timescale of the stars producing them. If all estimators are affected in reality at various degrees from dust obscuration or assumptions on the IMF, the Lyman continuum faces specific problems. The first issue is that it is difficult to measure observationally in the context of large surveys of high-redshift galaxies, be it from free-free emission for lack of sensitivity or from recombination lines. Indeed, measuring recombination lines requires spectra or specific narrow-band filters that severely limit the redshift range. In addition, at high redshift the Lyα line is the only widely usable line but it is very sensitive to resonant scattering by neutral hydrogen (e.g. Hayes et al. 2013), inducing large variations in the line strength even at fixed SFR. Another issue that has surfaced in recent years is the sensitivity of SFR calibration to models. In particular, the inclusion of stellar rotation can have a strong impact on the production of Lyman continuum photons (Levesque et al. 2012). If unconstrained such issues may outweigh the strong advantage of the Lyman continuum outlined in this paper, in which we have avoided the issue of model uncertainties on purpose.

Fig. 8 Same as Fig. 3 with the estimated SFR corrected for the contamination by long-lived stars. For reference, uncorrected luminosities are indicated with dashed lines. We see that once corrected, the estimated SFR are in much better agreement with the true SFR and that a large part of the bias has been eliminated. This is especially visible for the U band.

Even though it is much more affected by the variation of the SFR than the Lyman continuum, the FUV emission provides us with the second most accurate SFR estimator of the set we consider. This being said, the FUV overestimates the SFR by ~25% on average. The origin of this excess provides us with two important lessons to keep in mind for measuring the SFR in large surveys. First, unless the galaxy is currently undergoing a strong episode of star formation that completely outshines older stellar populations in star formation tracing bands, these long-lived stars need to be taken into account. In the case of the MIRAGE simulations, correcting for the presence of stars living more than 100 Myr lowers the bias by a factor ~3 in the FUV band and by a factor ~7 in the U band. We therefore recommend adopting SFR estimators calibrated over timescales longer than the standard 100 Myr. As we can see in Fig. 10, taking SFR estimators computed assuming a constant SFR over 1 Gyr would lower the bias from ~25% to ~10% in the FUV band as the difference in the calibration coefficients is ~10%. In the U band the difference is even more important, the calibration coefficient differing by ~50%. The second lesson, which is tied to the first lesson, is that the assumption of a constant SFR proves to be too simple. In the MIRAGE simulations there is a strong initial episode of star formation. This rapidly accumulates long-lived stars that contaminate star formation tracing bands for hundreds of Myr afterwards. Even if we have discarded this initial burst, it shows that strong episodes of star formation have lasting effects on the measure of the SFR. This exemplifies that variations of the SFR break the assumption of a constant SFR over long periods and that this has clear consequences. This calls for more realistic SFH to be considered or systematic biases may be left uncorrected. In other words, this recommendation to use SFR estimators calibrated on long timescales is valid for main-sequence galaxies with long star formation episodes and not for starbursting galaxies for which the relative contamination will be much lower and for which calibration on short timescales are adequate (Otí-Floranes & Mas-Hesse 2010). That is, we strongly recommend that SFR estimators are calibrated assuming a SFH consistent with the one expected for the sample under study, otherwise biases and errors will be introduced.

If random errors can be a problem for the measurement of the SFR in individual galaxies, they will average out over large samples. However, such large samples will remain affected by systematic biases. There has been a long-standing, redshift-dependent, tension between the integral of the SFR density of the Universe corrected for the return fraction and the stellar mass density. In their review, Madau & Dickinson (2014) found an excess of 0.2 dex for the integrated SFR density. Numerous explanations have been put forward involving biases on the measurements of either quantity or in the extrapolation of luminosity functions (see the aforementioned review). The excess we find in the present study would explain a part of the discrepancy though not necessarily the entirety. They adopt SFR estimators calibrated over 300 Myr that reduce the bias we find, but do not eliminate it. However, the relatively smooth variations of the SFR over long periods in the MIRAGE simulations likely provide us with a lower bound of possible biases, with more violent variations yielding larger biases that cannot be corrected for by simply changing the assumption of the duration of the star formation episode. The use of numerical simulations may prove to be useful for determining new, statistically accurate SFR estimators across various redshifts and galaxy parameters. Such an approach was for instance adopted by Wilkins et al. (2012) using a semi-analytic model. Interestingly, contrary to what we find here they have determined that standard estimators from the FUV, NUV, and u bands provide unbiased estimates at z = 0, although with an uncertainty of 20% at z = 0, and 30% at z = 6. However, they find an evolution of the UV calibration coefficient with redshift, which may be due to the same processes that we have shown in the current study. Direct comparison between the two studies is difficult though. The fundamental difference is that they obtained the SFH of a large number of galaxies using the galform semi-analytic model (Cole et al. 2000; Baugh et al. 2005) and computed average calibration factors at a handful of redshifts whereas we followed a small sample of galaxies with their SFH determined over nearly ~700 Myr through the means of state-of-the-art hydrodynamical simulations with highly refined baryonic physics.

Fig. 9 Same as Fig. 4 with estimated SFR corrected for the contamination by long-lived stars. We see that the bias and the scatter are strongly reduced in all bands. This is especially visible in the case of the U band.

Fig. 10 Same as Fig. 4 with SFR estimators calibrated over 1 Gyr rather than 100 Myr. We see that this strongly improves the agreement between the estimated and the true SFR, and it provides estimates nearly as good as when correcting for the presence of stars older than 100 Myr. The remainder of the differences is likely due to variations of the SFR.

5. Conclusion

In this paper we have investigated the impact of the SFH properties on the measurement of the SFR in galaxies through classical estimators. To do so, we have used the state-of-the-art MIRAGE sample, which consists of a set of 23 detailed hydrodynamical simulations of gas-rich mergers and isolated galaxies, including highly refined star formation prescriptions. It is designed to reproduce the characteristics of MASSIV (Contini et al. 2012), a spectroscopic sample of 1 <z< 2 galaxies. Passed the initial relaxation, the SFH show little to no long-term evolution but have rapid variations of the SFR with an amplitude of σ_SFR ≃ 0.34 ± 0.06 × ⟨SFR⟩. Peak-to-peak variations are of the order of a few. We have then combined the simulated SFH with the cigale SED modelling code (Boquien et al., in prep.; Burgarella et al., in prep.) to compute the global spectrum for each simulation every Myr in an ideal case in order to isolate the impact of the variation of the SFR from any other potential effect: we took into account stellar populations and nebular emission free of dust obscuration while assuming a perfectly known IMF and no intrinsic uncertainties from the models.

Considering five different SFR estimators (Lyman continuum, FUV, NUV and U bands, and the TIR emission of buried galaxies), we have found that:

• 1.

  On average all estimators except for the Lyman continuumoverestimate the SFR from typically ~25% for SFR(FUV) to ~65% for SFR(U), when such estimators have been calibrated assuming a constant SFR over 100 Myr. This is chiefly because of the contribution of stars living longer than 100 Myr.

• 2.

  Rapid variations of the SFR contribute to an increase of the uncertainty on the instantaneous SFR but have little long-term effect on the measurement of the SFR.

• 3.

  Slow variations of the SFR on timescales of a few tens of Myr and longer lead both to an increase of the noise on the measurement of the instantaneous SFR, and to systematic overestimates of the SFR for all tracers except for the Lyman continuum.

• 4.

  The amplitudes of the systematic overestimates depend on the metallicity. For the same SFH, more metal-poor galaxies will be the most affected.

Finally, if we find that the Lyman continuum is the best tracer in the ideal case we have considered, it nevertheless suffers from specific uncertainties and observational difficulties. The widely used rest-frame FUV suffers from biases if applied naïvely. These biases can be at least partially compensated for by using specially tailored SFR estimators for given redshifts and galaxy properties. At the very least, for main-sequence galaxies we suggest adopting SFR estimators calibrated over 1 Gyr rather than the commonly used 100 Myr to take contamination by long-lived stars into account. Shorter timescales may still be considered for starburst galaxies or star-forming regions within galaxies for which the relative contamination from old stars should be lower.

Acknowledgments

We are particularly grateful to the referee for constructive and insightful comments that have helped clarify the article. We would like to thank Paola di Matteo for interesting discussions at the origin of this paper. MB would like to thank Robert Kennicutt and Michele Trenti for useful discussions that have helped improve this article.

References

Appendix A: SFR estimators

All Tables

All Figures

Fig. 1 Left: SFH of the MIRAGE G2_G2_90_90_90 merger simulation (Perret et al. 2014). The time resolution is 1 Myr. The first 250 Myr are not shown as the initial relaxation is handled (Sect. 2.1). The coloured arrows at 400, 500, and 600 Myr indicate at which time we compute the spectra shown in the right-hand panel from the Lyman break to the U band. These spectra are created with the cigale SED modelling code (Boquien et al., in prep.; Burgarella et al., in prep.). In the present case we assume a Chabrier (2003) IMF and a metallicity Z = 0.008. Nebular emission is also included. It takes into account recombination lines and the nebular continuum from free-free, free-bound, and two-photon processes. The continuum contributes most in the red part of the shown spectra. We have computed the spectra for all 23 MIRAGE simulations (three isolated galaxies and 20 interacting galaxies) every 1 Myr.

In the text

Fig. 2 Luminosity weighted ages versus the wavelength for a constant SFR over 100 Myr (solid lines) and 1 Gyr (dashed lines). The magenta and red lines indicate the luminosity-weighted age for the Lyman continuum and the TIR for a completely buried galaxy, respectively. We see that the typical timescale probed is not only strongly dependent on the estimator and on the wavelength, but also on the SFH, with longer wavelengths and the TIR being the most sensitive to the duration of star formation episode in this simple scenario.

In the text

Fig. 3 Left: true SFR (black), SFR(Ly) (magenta), SFR(FUV) (blue), SFR(NUV) (cyan), SFR(U) (green), and SFR(TIR) (red) for the MIRAGE G2_G2_90_90_90 simulation (Perret et al. 2014). Note that the true SFR and SFR(Ly) are nearly blended. Right: zoom between 350 Myr and 450 Myr. We see that SFR(Ly) follows the variations of the true SFR extremely well even for high frequency variations. SFR estimators at longer wavelengths are not able to capture these variations though, leading to important discrepancies in the estimates of the instantaneous SFR. Beyond the capture of short-term variations, we see that there is also a systematic offset that is visible in the FUV and NUV bands but which is especially prominent in the U band. For the present simulation, we obtain ⟨[SFR(FUV) − SFR(true)] /SFR(true)⟩ = 0.27 ± 0.30 between 250 Myr and 790 Myr with a peak of ~1.5. For the U band, we obtain ⟨[SFR(U) − SFR(true)] /SFR(true)⟩ = 0.70 ± 0.46 with a peak of ~2.8. This is representative of the entire MIRAGE sample as can be seen in Fig. 4.

In the text

Fig. 4 Ratio of difference between the estimated and true SFR to the true SFR versus time for the entire MIRAGE sample. The solid lines indicate the mean value across all simulations whereas the shaded areas of the same colour indicate the standard deviation. The dashed lines indicate the mean value between 250 Myr and 790 Myr. The colour coding is the same as for Fig. 3. We see that except for the Lyman continuum, the SFR is systematically overestimated, from 24% on average for the FUV band to 65% on average for the U band. Peak relative differences reach ~0.7 for the Lyman continuum, ~1.9, ~2.1, and ~3.2 for the FUV, NUV, and U bands, and ~2.0 for the TIR.

In the text

Fig. 5 Left: luminosity weighted ages versus the simulation time. The solid lines indicate the mean value of for the 23 simulations and the shaded areas of the corresponding colours indicate the standard deviation. The colour scheme is the same as that of Fig. 3. Right: luminosity weighted ages versus wavelength (black) for a simulation age of 790 Myr. The values for the Lyman continuum and the TIR at 790 Myr are also indicated. For comparison, the corresponding values for a constant SFR over 1 Gyr shown in Fig. 2 is indicated in dashed lines.

In the text

Fig. 6 Left: Fourier-filtered SFH of the MIRAGE G2_G2_90_90_90 simulation (black lines) along with the estimated SFR recomputed from each of the two filtered SFH. The top (resp. bottom) lines represent the low (resp. high) frequency component. The colours follow the same scheme as that of Fig. 3. Note that the true SFR and SFR(Ly) are nearly blended. Right: zoom between 350 Myr and 450 Myr. We find that if there are important discrepancies in the SFR, on average at high frequency the true SFR and the estimated SFR give similar results, both close to 0 Msun/yr. The behaviour at low frequency is markedly different. We see that the estimated SFR show both a systematic overestimate and a delay in time from the extrema of the true SFR. The smallest delay occurs for SFR(Ly). Conversely, SFR(FUV), SFR(NUV), SFR(U), and SFR(TIR) all show approximately similar delays of a few Myr.

In the text

Fig. 7 Fraction of the luminosity in star formation tracing bands contributed by stars living over 100 Myr versus simulation time. The colour scheme is the same as that of Fig. 3. The contribution of long-lived stars is negligible for the Lyman continuum but reaches upward 10% in the FUV, a small but non-negligible fraction. The contamination is stronger for the NUV and TIR tracers, and it is maximum for the U band with a typical value of ~35%.

In the text

	Fig. 8 Same as Fig. 3 with the estimated SFR corrected for the contamination by long-lived stars. For reference, uncorrected luminosities are indicated with dashed lines. We see that once corrected, the estimated SFR are in much better agreement with the true SFR and that a large part of the bias has been eliminated. This is especially visible for the U band.
In the text

	Fig. 9 Same as Fig. 4 with estimated SFR corrected for the contamination by long-lived stars. We see that the bias and the scatter are strongly reduced in all bands. This is especially visible in the case of the U band.
In the text

	Fig. 10 Same as Fig. 4 with SFR estimators calibrated over 1 Gyr rather than 100 Myr. We see that this strongly improves the agreement between the estimated and the true SFR, and it provides estimates nearly as good as when correcting for the presence of stars older than 100 Myr. The remainder of the differences is likely due to variations of the SFR.
In the text