© ESO, 2012

1. Introduction

One of the important, recent results of multi-wavelength galaxy surveys is the finding of a well-defined relation between the star formation rate (SFR) and stellar mass of galaxies – now often called the main sequence of star forming galaxies – and the evolution of this relation with redshift, at least from the nearby Universe up to redshift z ~ 2 (Noeske et al. 2007; Daddi et al. 2007; Elbaz et al. 2007; Rodighiero et al. 2011). The observed redshift evolution corresponds to an increase of the typical specific star formation rate (sSFR) by more than an order of magnitude from z ~ 0 to 2, indicating a higher star formation activity of galaxies in the past.

Recent analysis of Lyman break galaxies (LBGs), combining often HST and Spitzer observations, have tried to examine whether a SFR-mass relation is already in place at higher redshift, and how the specific star formation rate of galaxies behaves at z ~ 4–7 (e.g. Stark et al. 2009; Schaerer & de Barros 2010; Labbé et al. 2010; McLure et al. 2011; González et al. 2011). The first studies have found a relation similar to the z ~ 2 SFR-mass relation, which would indicate no evolution of the sSFR beyond redshift  ≳ 2 – a plateau in sSFR – (e.g. Stark et al. 2009; González et al. 2011). This result appears difficult to reconcile with most theoretical models, which successfully explain the behavior of star formation properties from z ~  0 to 2, but predict a continuing rise of the specific star formation rate towards higher redshift (Bouché et al. 2010; Dutton et al. 2010; Weinmann et al. 2011; Davé et al. 2011), although others propose alternatives (e.g. Krumholz & Dekel 2012). More recent work, including different spectral energy distribution (SED) models, or allowing for a revision of dust attenuation have shown that the observed sSFR at z ~ 5–7 could be higher than previously thought (Schaerer & de Barros 2010; de Barros et al. 2011, 2012; Bouwens et al. 2012; Yabe et al. 2009).

If all galaxies at z ≳ 2 obey a mass–SFR relation with star formation rate increasing with galaxy mass and its normalization remains basically unchanged with redshift, it is evident that their star formation rate must increases with time. The apparent, small scatter in the SFR-mass relation, and the value of the exponent α close to unity in this relation SFR found by various studies (e.g. Daddi et al. 2007; Elbaz et al. 2007; González et al. 2011) has led several authors to conclude that high redshift galaxies must have gone through a phase of quasi-exponential growth (e.g. Renzini 2009; Maraston et al. 2010). Earlier studies have independently advocated rising star formation histories (SFHs) for high redshift galaxies from hydrodynamic simulations and semi-analytical models (Finlator et al. 2007, 2011). Arguments in favor of rising SFHs, at least on average, have been put forward by Finkelstein et al. (2010); Papovich et al. (2011) to explain the evolution of the UV luminosity function. If representative of the star formation history of individual galaxies, such star formation histories would indeed by quite different from exponentially declining or constant SFHs, most frequently assumed in the literature to model/analyze the observed SED of distant galaxies.

Assumptions on the star formation history affect the physical parameters of galaxies derived from SED fits. This is therefore generally treated as a free parameter. See e.g. Reddy et al. (2010); Maraston et al. (2010); Wuyts et al. (2011); Reddy et al. (2012b) for examples on z ≳ 2 star forming galaxies. Determinining star formation histories and their associated timescales is also important, as this may provide constraints on different modes of star formation and on feedback processes at high redshift, recognized as key features for our understanding of galaxy evolution. For example, star formation may proceed on different timescales if regulated by cold-accretion or by star formation feedback or triggered by mergers (cf. Khochfar & Silk 2011; Wyithe & Loeb 2011). These points already illustrate the interest and importance of determining the star formation histories and typical timescales of star formation at high z. We here wish to discuss and reexamine these issues on the basis of an analysis of a large sample of LBGs covering redshifts from  ~3 to 6, and to present observational tests which should be able to distinguish different star formation histories.

One of main arguments often invoked to argue for rising star formation histories is the small scatter of the SFR-mass relation. However, it should be recognized that generally the SFR is derived from observables – the UV or IR luminosity – which depend on relatively long (≳100 Myr) timescales, and which are assumed to be at an equilibrium value, which is only reached after this timescale and for constant SFR. With these assumptions entering e.g. the commonly used SFR(UV) or SFR(IR) calibrations or Kennicutt (1998), it is natural that the “observational” scatter is smaller than the true scatter in the current SFR, if typical ages are less than 100 Myr and/or the timescale shorter than this. Especially at high redshift, where timescales are shorter (e.g. the dynamical timescale decreases with (1 + z) ^− 3/2, cf. Wyithe & Loeb 2011) one should therefore carefully (re)examine the SFR-mass relation and its “tightness” using consistent diagnostics, and with all the available observational constraints including on age, star formation timescales and histories. In any case, the conclusion that rising SFHs are favored is a priori inconsistent with the assumptions on the SFH for the determination of the SFR-mass relation and its scatter, where most studies simply assume a constant star formation rate. Although this has been re-examined for z ~ 2 samples (Maraston et al. 2010; Wuyts et al. 2011; Reddy et al. 2012b) this obviously calls for a revision at higher redshifts.

Most SED studies of LBGs at z ≳ 3 have assumed constant star formation rates, or exponentially declining SFHs to determine the physical parameters of these galaxies (e.g. Egami et al. 2005; Schaerer & Pelló 2005; Eyles et al. 2005; Verma et al. 2007; Yabe et al. 2009; Stark et al. 2009; Lee et al. 2010; Schaerer & de Barros 2010; González et al. 2011), some of them imposing no dust attenuation, motivated by the blue UV slopes observed at high redshift. Notable examples are the work of Finlator et al. (2007) who analyzed 6 galaxies at z > 5.5 with different SFHs, including rising ones taken from they hydro-dynamic simulations. The physical parameters they derive are consistent with those using simple parametrised histories, albeit with reduced uncertainties. Their study also shows that SED fits with rising star formation histories yield a higher attenuation than inferred assuming constant SFR, a result not yet appreciated enough, which we also find in de Barros et al. (2012) and in this paper. Most recently, other groups have also analyzed high-z samples with rising SFHs (e.g. Curtis-Lake et al. 2012; Gonzalez et al. 2012). However, dust extinction is, for example, not treated consistently with the star formation history in the approach of Gonzalez et al. (2012).

Another drawback of most earlier studies is that the contribution of nebular emission (most emission lines) is not taken into account, an effect which can significantly alter the ages, masses, and other physical parameters derived from SED fits, as demonstrated by Schaerer & de Barros (2009, 2010). Indeed, there is now clear evidence for the presence of nebular lines affecting the broad-band photometry of high-z star-forming galaxies (Lyman alpha emitters and LBGs in particular), as discussed e.g. in Schaerer & de Barros (2012). The best demonstration comes from LBGs with spectroscopic redshifts between 3.8 and 5, among which a large fraction shows a clear excess in the 3.6 μm filter with respect to neighboring filters (K and 4.5 μm), as shown by Shim et al. (2011). At these redshifts Hα is located in the 3.6 μm filter, whereas very few lines are expected at 4.5 μm, showing that the observed 3.6 μm excess is naturally explained by strong Hα emission.

Given these limitations of published SED studies and the interesting results obtained from our work on a small sample of z ~ 6–8 LBGs (Schaerer & de Barros 2010), we have recently undertaken an extensive study of the physical parameters of a large sample of  ~1400 LBGs at z ~ 3–6, using our state-of-the art photometric redshift and SED fitting model including nebular emission, and considering a range of different star formation histories. Among the numerous detailed results obtained in this work (de Barros et al. 2012, hereafter dBSS12) we find in particular from the preferred models: 1) a large scatter in the SFR-mass relation for the preferred, variable star formation histories, 2) higher dust attenuation than obtained from models assuming SFR = const. and from standard methods using the UV slope, and 3) a higher sSFR than commonly obtained, and a rising sSFR with redshift. Our models therefore reconcile the observationally-inferred specific star formation rate with predictions from cosmological models predicting a continuous rise of the sSFR with redshift (Bouché et al. 2010; Dutton et al. 2010; Weinmann et al. 2011; Davé et al. 2011; Krumholz & Dekel 2012). The results from dBSS12 have also other important implications, e.g. on the typical star formation timescales at high redshift, and on the cosmic star formation rate density.

Given the importance of these findings, it is of interest to carry out additional independent tests of the models and to provide further constrains on the ages, extinction, and star formation histories of high redshift galaxies. In this paper we present predictions allowing such tests. Furthermore, we extend the study of dBSS12 by exploring other star formation histories, not considered there. Indeed, while dBSS12 adopted the fixed, rising star history of Finlator et al. (2011) obtained from the hydrodynamic simulations, we here explore exponentially rising SFHs and so-called “delayed” histories with variable timescales. These histories were previously applied to the SED fits of galaxies at lower redshift, e.g. at z ~ 2 by Maraston et al. (2010), Wuyts et al. (2011), Reddy et al. (2012b) and found to be preferred over other simple star formation histories.

First, we examine the effect of the different SFHs on the derived physical parameters of LBGs. We then present the implications these different model assumptions have on the SFR-mass relation. After that we present consistent predictions for the infrared luminosity of these galaxies, based on their SED fits and the assumption of energy conservation, abandoning the standard SFR–L_IR conversions. We show in particular that such a consistent prediction yields a smaller scatter in the observables (L_IR) than in the current SFR. The true star formation rate of LBGs at z ≳ 3 could thus very well show a large scatter around a “main sequence” even if the UV and/or IR luminosities show a small scatter. Our models also show that star formation histories can, to some extend, be distinguished from measurement of IR luminosities, which will become possible with ALMA, since different amounts of UV attenuation are expected for different SFHs. We also present the predicted strength of some selected emission lines, which can be used to test our models and the different star formation histories.

Our paper is structured as follows. The observational data and the method used for SED modelling are described in Sect. 2. The dependence of the physical parameters on the assumed star formation histories is discussed in Sect. 3. Our general predictions for the IR emission are presented in Sect. 4. Specific predictions for z ~ 4−6 LBGs and ALMA are given in Sect. 5. The predicted strengths of optical emission lines are shown in Sect. 6. We discuss our results in Sect. 7, and Sect. 8 summarises our main conclusions. We adopt a Λ-CDM cosmological model with H₀ = 70 km s^-1 Mpc^-1, Ω_m = 0.3 and Ω_Λ = 0.7.

2. Observational data and SED modelling

de Barros et al. (2012, hereafter dBSS12) have analysed a large sample of z ~  3–6 dropout-selected galaxies in depth using an up-to-date photometric redshift and SED-fitting tool, that treats the effects of nebular emission on the SEDs of galaxies. In their homogeneous analysis they determine the main physical properties, such as the star formation rate (SFR), stellar mass, age, and reddening. They assess carefully their uncertainties, and discuss the evolution of these properties with redshift. We here extend these simulations to include other star formation histories and we examine the predicted IR luminosities and equivalent widths of selected emission lines from these galaxies.

2.1. Photometric data and sample selection

We have used the GOODS-MUSIC catalogue of Santini et al. (2009), which provide photometry in the U, B₄₃₅, V₆₀₆, i₇₇₆, z_850LP, J, H, and K bands mostly from the VLT and HST, and the 3.6, 4.5, 5.8, and 8.0 μm bands from the IRAC camera onboard Spitzer. Using standard criteria as in Nonino et al. (2009) and Stark et al. (2009) we then selected U, B, V, and i-drop galaxies. To reduce the contamination rate (typically  ~ 10–20%) we only retained the objects whose median photometric redshifts agree with the targetted redshift range. This leaves us with a sample of 389, 705, 199, and 60 galaxies with median photometric redshifts of z_phot = 3.3, 3.9, 4.9, and 6.0. Typicallly, the completeness limit of our sample is M₁₅₀₀ ≈  −19 (corresponding to L_UV ≈ 8 × 10⁹ L_⊙) at z ~ 4 and up to M₁₅₀₀ ≈  −19.8 at z ~ 6. It is discussed in Stark et al. (2009). See dBSS12 for more details.

2.2. SED models

Our SED-fitting tool, described in Schaerer & de Barros (2009) and (2010), is based on a version of the Hyperz photometric redshift code of Bolzonella et al. (2000), modified to take nebular emission into account. In dBSS12 we considered a large set of spectral templates based on the GALAXEV synthesis models of Bruzual & Charlot (2003) and covering different metallicities and a wide range of star formation histories. A Salpeter IMF from 0.1 to 100 M_⊙ was adopted. Our results can easily be rescaled to other IMFs with the same power law slope at high masses. Nebular emission from continuum processes and numerous emission lines were added to the spectra predicted from the GALAXEV models as described in Schaerer & de Barros (2009), proportionally to the Lyman continuum photon production. Following the results of dBSS12 we set the Lyα line to zero for z ~ 3–5, and we assume the case B value for the i-drop sample. A more detailed treatment of Lyα, following e.g. Schaerer et al. (2011), is not necessary and clearly beyond the scope of this paper, as the main physical parameters show generally only a small dependence on the treatment of Lyα. The intergalactic medium (IGM) was treated with the prescription of Madau (1995).

The free parameters of our SED fits are: redshift z, metallicity Z (of stars and gas), the age t_⋆ defined since the onset of star-formation, the star formation history, and attenuation A_V described by the Calzetti law (Calzetti et al. 2000). Stellar mass and the current star formation rate are derived from the absolute scaling of the SEDs. For comparison, a subsample of the data was also modeled using the SMC extinction law (cf. Sect. 7).

Since one of the main objectives of this paper is to explore a variety of star formation histories (SFHs), we have adopted five different histories, which are summarized in Table 1, and illustrated in Fig. 1. Note that for exponentially rising SFHs the SFR is set to zero, after a growth by more than 20 decades. The bolometric luminosity output per unit SFR, respectively its inverse SFR/L_bol, and the mass-to-light ratio in the V-band corresponding to these simple star formation histories are shown in Figs. 2 and 3. The behavior of these quantities are fundamental to explain the main differences in star formation rate and stellar mass obtained from different models, discussed in Sect. 3.

The first three cases, constant star formation, exponentially declining SF, and the average rising SFH predicted by the hydrodynamical models of Finlator et al. (2011) have already been used by de Barros et al. (2012). Two additional SFHs, exponentially rising, and so-called “delayed” SF showing a initial increase and subsequent decrease of the star formation, are also considered here. The way these assumed SFHs modify the derived physical parameters of LBGs and how they could be distinguished is discussed below.

In practice we compute for each model set (i.e. star formation history) SED fits for all combinations of z ∈  [0, 10]  in steps of 0.1 Z =  (0.02 = Z_⊙, 0.004, 0.001), τ =  (0.01, 0.03, 0.05, 0.07, 0.1, 0.3, 0.5, 0.7, 1., 3., ∞) Gyr for model A¹, 51 age steps from 0 to the age of the Universe (see Bolzonella et al. 2000), A_V ∈  [0,2]  mag in steps of 0.05. Minimisation over the entire parameter space yields the best-fit parameters and SED, along with other properties such as the stellar mass and star formation rate (SFR) and UV magnitude. Other physical parameters, briefly described below, are determined from the best-fit SED. To determine the uncertainties of the physical parameters, we used Monte-Carlo simulations to generate 1000 realizations of each object perturbing the photometry in all filters according to the tabulated errors. Each of these 1000 realizations of each galaxy was fit. From the best-fit results for each realization we derived the probability distribution function for each parameter/quantity, either for each individual object or for (sub)samples.

Simplifying assumptions made for our SED fits and possible caveats are discussed in Sect. 7.4.

2.3. Predicted physical parameters and observables

In addition to the fit parameters of our models, redshift z, metallicity Z, age t, SF timescale τ, stellar mass M_⋆, and current star formation rate SFR, our models also allow us to determine other parameters and observables of interest. We briefly describe them here, for sake of clarity. Note that for all the physical parameters and observables we determine the detailed probability distribution function (pdf) from the full MC simulations. From the 1D pdf we derive in particular the median value and the 68% confidence interval, which we use to illustrate the uncertainties on these parameters. For most of the paper we use for simplicity the median value of the physical parameter of interest.

2.3.1. Emission lines

From the numerous emission lines included in our spectral templates we select some lines, primarily Lyα, Hα, Hβ, [O ii] λ3727, [O iii] λλ4959,5007, for which we save the emission line fluxes and equivalent widths predicted by the photometric SED fits. These lines can subsequently be used for comparison with spectroscopic observations (existing or future).

2.3.2. UV and IR emission

For straightforward comparison with observations we use the absolute UV magnitude M₁₅₀₀ defined at 1500 Å, and derived from the SED using a square filter of 300 Å width centered on this wavelength. To define an emergent (observed) UV luminosity L_UV we follow common practice by computing L_UV = λF_λ, where λ = 1800 Å, and F_λ is the average flux between 1400 and 2200 Å.

Assuming energy conservation, i.e. that the attenuated/extinct stellar light is reemitted by dust in the IR, we can predict the amount of IR radiation expected from each galaxy. The corresponding IR luminosity L_IR is simply computed from the difference between the intrinsic, unattenuated SED given by the spectral template, and the best-fit SED with the parameters. In practice L_IR is calculated from the attenuated energy integrated from 913 Å to 3 μm following the same procedure as SED models carrying out complete from the UV-visible-IR fits (cf. da Cunha et al. 2008; Noll et al. 2009). This quantity also includes the nebular emission (lines and continuum) if present in the SED, consistently with our assumption of identical attenuation for both stars and gas. We have checked that the integral over the Balmer continuum and up to 3 μm yield basically identical results.

The method used to compute the IR luminosity does not make any explicit assumption on the origin of the stellar sources contributing to it and on the timescale over which the IR emission is produced. All the energy emitted between the UV and 3 μm and attenuated by dust according to the adopted extinction/attenuation law contributes to the IR emission, irrespective of stellar type, age etc. In practice, however, the timescale of the IR emission is not very different than that of UV emission for the cases modeled here, since the bulk of the emitted energy originates from the UV, where the attenuation is also strongest. In models assuming a specific dust geometry, time-dependent attenuation, multiple populations or alike (e.g. the models of da Cunha et al. 2008), the timescale of UV and IR emission could be different. However, such models include additional free parameters, and have, to the best of our knowledge, not yet been used to fit distant galaxies. For simplicity we here fit the SEDs with single stellar populations and a single extinction/attenuation law.

From the UV and IR luminosity just defined, we also compute L_IR + UV, the emergent luminosity in the UV + IR domain, which is a good proxy of the total bolometric luminosity and to some extent also of the SFR (but cf. below). Similarly, we will use the ratio of the IR to UV luminosity, L_IR/L_UV, a quantity known to be a good measure of dust attenuation (cf. Meurer et al. 1999).

2.3.3. Predicted IR, sub-mm, and mm fluxes

Given the IR luminosity L_IR predicted consistently from the amount of attenuation by dust, it is straightforward to predict the flux at various IR wavelengths, provided a corresponding spectral template for dust emission is adopted. We simply assumed standard modified black-body spectra described by F_ν ∝ ν^βB_ν(T_d), with a fixed value of β = 2 and three different dust temperatures T_d = (25,35,45) K, which span the range of observed dust temperatures in nearby and distant galaxies. The modified black-body is scaled to L_IR computed from 8 to 1000 μm (rest frame). From this, we determine the predicted IR fluxes in bands corresponding to observations with Herschel (PACS and SPIRE bands), several sub-mm bands (350, 450, 850, 870 μm), 1.2 and 2 mm bands, and ALMA bands 1 − 10 (corresponding to  ~8.1 to 0.35 mm). Empirical or semi-empirical spectral templates covering the IR domain could also be used to convert the predicted IR luminosities to observed fluxes. However, since it is not known how well the empirical templates based on low redshift and/or high IR luminosity galaxies describe typical LBGs at z > 3, and to use a simple prescription, we have used modified black-body spectra. In any case, our conclusions rely on the total IR luminosity, which is template independent. Of course, the predicted L_IR depends indirectly on the stellar templates and attenuation law assumed for the SED fits.

3. Dependence of the physical parameters on the assumed star formation histories

To illustrate and discuss the dependence of the physical parameters on the assumed star formation histories we here choose the largest sample comprising 705 B-drop galaxies. The general behavior discussed here is also found for most of the other samples at z ~ 3, 5, and 6. Some results for other redshifts are discussed in Sect. 5. First we show briefly how the SFHs affect the derived physical parameters. Then we compare the fit quality of the models, and discuss implications on the SFR-mass relation.

3.1. Comparison of physical parameters

It is well known that the physical parameters determined from broad-band SED fits depend on assumptions such as the star formation history. In dBSS12 we have shown how quantities such as the stellar mass, age, dust reddening, and SFR depend on this assumption and how they are modified with the inclusion of nebular emission. Here we briefly illustrate the dependence of the main physical parameters on SFH, including in particular the new histories explored in this paper (exponentially rising and delayed SFH; cf. Table 1). The effects of nebular emission on the SED fits and the resulting physical parameters have been discussed in detail in dBSS12. In particular, we have shown there that nebular emission affects approximately two thirds of the galaxies at all the redshifts considered (z ~ 3−6), quite irrespectively of their UV magnitude and galaxy mass. For more details see dBSS12.

Figures 4 to 7 show how the M_⋆, age, A_V, and SFR derived for different assumptions on SFH compare to models assuming exponentially declining star formation histories (our model A) for our largest sample, the 705 B-drop galaxies. This model is taken as a reference, since it generally provides the best fits for the vast majority of galaxies analyzed here (cf. below). For all models except for SFR = const. (model B), we include the effects of nebular emission. Model B, called REF in dBSS12, serves as a “reference” model to ease comparisons with other studies in the literature assuming constant star formation, and neglecting nebular emission. When nebular emission is added to model B, this does not lead to very different parameters (cf. dBSS12).

Figure 4 shows that most models yield similar stellar masses, as already known from earlier studies (Yabe et al. 2009). The largest differences are found with respect to models assuming constant star formation and a mimum age of 50 Myr (model B), which yield generally higher stellar masses. On average we find masses lower by a factor  ~3 at z ~ 4 for declining SFHs (model A), where half of this effect is due to inclusion of nebular emission (cf. de Barros et al. 2012). The difference is mostly due to the average shift towards younger ages, which have a lower mass-to-light ratio, as shown in Fig. 3. It is immediately apparent that the models with delayed star formation histories (model E) yields very similar results to model A. Exponentially rising SFHs with variable timescales (model D) and the rapidly rising history from Finlator et al. (2011) (model C) give also comparable stellar masses, although with a larger scatter, when compared to model A. In short, all the models with variable SFHs and nebular emission considered here yield quite similar stellar masses, which are on average lower than those obtained from standard models without nebular emission.

Age is notoriously known to be sensitive to models assumptions (e.g. Papovich et al. 2001; Yabe et al. 2009). Indeed, absolute ages predicted by different models can differ by orders of magnitude, as shown in Fig. 5. Overall models with delayed SF (E) show again quite similar results to those with exponentially decreasing histories (A). The behavior of the rising SFHs (C and D) is quite comparable. By construction, model B yields only stellar ages  ≥ 50 Myr² Since recent populations dominate for exponentially rising star formation, ages are more difficult to constrain. This is reflected in increased uncertainties, e.g. by wider confidence intervals (not shown here), as we have verified from our models. Whereas for exponentially declining and delayed SFHs most objects are best fit with relatively short and comparable timescales τ, the opposite is true for fits with exponentially rising histories. In this case the probability distribution function for τ clearly peaks at the maximum value allowed here (τ = 3 Gyr), which means that a relatively slow growth is preferred for most galaxies. This seems quite natural given the finding that the rising SFH from Finlator et al. (2011) fits on average less well than declining models.

The attenuation predicted by the various models is shown in Fig. 6. As already discussed by dBSS12, constant SF models show generally lower dust attenuation than the exponentially declining models, whereas the rising SFH of Finlator et al. (2011) predicts a higher attenuation on average. Rising star formation histories with variable timescales globally scatter around those assuming declining SF. Again, delayed SFHs give results very similar to those with declining SFHs.

The predicted median SFR of all galaxies are compared in Fig. 7. As expected due to lower dust attenuation and older ages, the models with constant SFR show generally the lowest (current) star formation rates (cf. dBSS12). As before, models A and E are very similar. Finally, assuming rising star formation histories generally implies a higher SFR, due to, first, the lower L_UV output per unit SFR (cf. Fig. 2), and, second, to a higher average dust attenuation, needed to redden the UV-to-optical colors, which are always dominated by the youngest stars. Again, this result is discussed in dBSS12, and has already been pointed out earlier by Schaerer & Pelló (2005) and shown in the analysis of Finlator et al. (2007).

For all models with no age constraint, fits at young age will generally push the derived SFR upward. For this reason we discuss below (Sect. 3.3) how the SFR changes when an age prior is adopted. The age distributions of our samples found for different SFHs and at different redshifts are extensively illustrated in Paper II (dBSS12), which also discusses possible caveats about the age. Overall, our models show a trend of increasing median age with decreasing redshift (e.g. at fixed absolute UV magnitude). At z ~ 3 e.g., the typical median ages are 30 (90) Myr for models with declining (rising) SFH and including nebular emission, compared to  ~250 Myr for model B (constant SFR, no lines), comparable to others in the literature. Obtaining independent age constraints should also help testing the various models, and narrowing down the SFR differences between models.

3.2. Comparison of the fit quality

Figure 8 shows the fit quality of the models expressed by the reduced value³ for the z ~ 4 sample. Overall, declining star formation histories (model A) give the best fit to most galaxies, and delayed SFHs very similar fits. Exponentially rising histories provide better fits for some objects, but a slightly higher χ² on average. Models B and C with constant SFR and rapidly rising star formation yield fits of lower quality for most galaxies. Overall the significance of the fit quality is relative small, as can be seen by the χ² differences (Fig. 8). Since model B (constant SFR) is a subset of model A (exponentially declining), the difference between two models can e.g. be compared using the likelihood ratio test. The average (median) χ² difference between the two models corresponds to a probability of  ~25–30% for model B being a better fit than model A. Statistically speaking, based on the fitted part of the SED, it is therefore not possible to firmly distinguish different SFHs. Other constraints, as the IR emission or emission line measurements, will be needed, as discussed in the rest of the paper. See also dBSS12 for other arguments on different star formation histories.

As also done in other studies, we have examined the influence of another reddening law (cf. Sect. 7). Applying the SMC law of Prevot et al. (1984); Bouchet et al. (1985) we find that the majority of galaxies fit less well than with the Calzetti law, both the declining and rising star formation histories (model A, C, and D), although the differences are generally relatively small (cf. Sect. 7). Some studies have favored the SMC law for a fraction of LBGs e.g. for young objects (cf. Reddy et al. 2010; Shim et al. 2011), whereas usually the Calzetti law is thought to be a good description of the observations.

3.3. Implications for the SFR-mass relation and influence of the minimum age

Figure 9 shows the SFR-mass relation for the 705 B-drop galaxies as derived from the SED fits assuming the five different star formation histories. As already discussed in dBSS12, the models show overall a large scatter and a median specific SFR (SFR/M_⋆) above the mean relation derived by Daddi et al. (2007) at z ~ 2, also shown here comparison with other studies and redshifts. In part the scatter in this diagram is determined by the possible range of SFR/M_⋆ predicted for stellar populations with different SFHs and different ages. Indeed, for a given L_bol/L_V, probed observationally by L_UV/L_V, the ratio SFR/M_⋆ varies by  ~2 orders of magnitude between young ages and t ≲ 1 Gyr for constant SFR models, and more for declining or delayed star formation histories (Figs. 2 and 3). Of course, since reddening also modifies L_UV/L_V, this adds further scatter.

As expected, the delayed SFH yields basically identical results to the declining histories. Rising histories (model D) occupy a similar domain as the other variable SFHs (models A, C, E), except for the absence of galaxies below the mean relation from Daddi et al. (2007). This is simply due to the fact that star formation cannot decrease in this model and for model B (constant SFR). Hence situations where the current SFR is low, but the UV flux is still present from somewhat older stars (e.g. B-type stars with lifetimes up to  ~100 Myr) are not possible for these SFHs (cf. Fig. 2). Objects well below the SFR-mass relation require obviously decreasing star formation histories, or some shut-off of star formation.

The upper boundary of the SFR–M_⋆ plot, parallel to the indicated z ~ 2 relation, is determined by the minimum age, and by dust attenuation. Indeed, young ages imply lower masses and a higher SFR from simple considerations of SFR/UV and mass/light ratios (cf. Figs. 2 and 3, also McLure et al. 2011).

Since age is kept free in all our models⁴ except for model B where we assume a minimum age of 50 Myr, and since some the median ages obtained are relatively young (cf. dBSS12) it is interesting to explore how the main physical parameters and the SFR–M_⋆ relation change when introducing a prior on the minimum age (or equivalently on the formation redshift) for all models. The physical motivation for such a minimum age of 50 Myr is to avoid young ages, which may be unrealistic for the entire stellar populations of a galaxy, in particular in view of longer dynamical timescales (e.g. Reddy et al. 2010).

The effect of fixing the minimum age to t_min = 50 Myr in the SED fits with exponentially declining and rising star formation histories (models A and D) on the SFR-mass relation is shown in Fig. 10. In short, for the declining SFHs the upper boundary becomes now identical to the one obtained with SFR = const. (model B), as expected, since declining SFHs include long timescales τ for which the star formation history becomes undistinguishable from constant, and since t_min = 50 Myr sets the maximum SFR/UV output. In contrast, for rising SFHs the upper boundary is found to be higher by a factor  ~3 approximately. Again, this simply due to the less efficient UV output per unit SFR (or higher SFR/UV ratio) of the youngest stars, which always dominate the SFR for such star formation histories (Fig. 2). The lower boundary of the SFR-mass relation remains of course unchanged by changes of the minimum age. In conclusion, introducing a minimum age affects in different ways the SFR–M_⋆ relation, depending on the SFH.

Among the z ~ 4 sample, 50 − 70% have median ages less than 50 Myr for models with nebular emission and for the different SFHs. As already mentioned, the uncertainties on ages are large. The median ages are also found to increase with decreasing redshift and with increasing mass, as shown in dBSS12. On average, i.e. for the entire sample of 705 B-drop galaxies, the median age and mass increase by a factor  ~2, whereas SFR and A_V decrease by a factor  ~2, when the prior t_min = 50 Myr is introduced. For the rising histories (model D) the median age increases by a factor  ~5, SFR decreases by 1.5, stellar mass increases by a factor 2. On the other hand the dust attenuation, both of individual galaxies and on average, remain very similar, with interesting implications on the predicted IR luminosity (see Sect. 4.4).

To conclude this section we note that SED fits using models with variable star formation histories imply a fairly large scatter in the SFR-mass relation for LBGs at z ≥ 3, as shown here and in dBSS12. However, it must be pointed out that the amount of scatter cannot directly be translated into observable quantities such as UV and/or IR luminosity, which are commonly used as SFR indicators. Indeed, as we will show below (Sect. 4.3), a much smaller scatter is expected in terms of observed UV + IR luminosity.

4. Predicted IR-mm emission

As shown in the previous section, the physical parameters inferred from SED fits depend in particular on the assumed star formation histories and on possible age priors (assumptions of a minimum age). We now present the IR properties predicted from the SED fits and discuss their dependence on these parameter. The effects of different extinction laws are addressed in Sect. 7.

4.1. IR luminosity and constraints on UV attenuation by dust

Figure 11 shows the infrared luminosity predicted by the different models for the z ~ 4 sample, Fig. 12 the corresponding histogram of L_IR/L_UV⁵. Since some galaxies have a median A_V = 0 and hence no UV light is absorbed, their predicted median IR luminosity is L_IR = 0. Such galaxies are not shown on these plots; their percentage ranges from 24% for model B (with constant SFR) to 10% for model C, which has the highest median L_IR. Due to the adopted discretisation in A_V, the lowest, non-zero extinction A_V = 0.05 corresponds to log (L_IR/L_UV) ≈  − 0.85, as seen e.g. in Fig. 27.

As expected from the behavior of the UV attenuation discussed above and in dBSS12, the models with constant SFR predict generally significantly lower L_IR than models with different star formation histories. On average, the highest L_IR is predicted for steeply rising SF of Finlator et al. (2011), followed by the exponentially rising model D, and model A and E with very similar average IR luminosities. Clearly, if L_IR can be determined with sufficient accuracy from observations, it should be possible to distinguish different models for individual galaxies, as well as statistically for large enough samples.

Indeed, the IR/UV ratio allows one to measure the UV attenuation as well known (cf. Buat et al. 2005; Iglesias-Páramo et al. 2007, and references therein), and this independently of the assumed SF history and of the extinction law. The corresponding relation between the UV attenuation factor f_UV and the IR/UV ratio is shown in Fig. 13 for all galaxies of the z ~ 4 sample, for all star formation histories and for the Calzetti attenuation and SMC extinction curve. Here f_UV is defined by the ratio of the intrinsic (emitted) over the observed UV luminosity, i.e. (1)The UV attenuation factor (defined here at 1800 Å) can simply be described by (2)where x = log (L_IR/L_UV). For a specific attenuation law this can of course be translated into quantities such as A_V. For example, for the Calzetti law one has A_V = 2.5(R_V/k_λ)log f_UV = 1.08log f_UV. Such a simple relation is quite naturally obtained, since the bulk of the energy reemitted in the IR after absorption by dust is also produced by the same stars contributing the bulk of the UV emission. More complicated relations could be expected when e.g. older stellar populations such as AGB or RGB stars contribute to dust emission, as appears to be the case in low redshift galaxies (Buat et al. 2005).

In short, measurements of L_IR/L_UV of individual LBGs provide a direct measure of their UV attenuation, which should allow us to distinguish certain star formation histories, and in particular to test the systematic shift in physical properties predicted between models assuming long star formation timescales (constant SFR) and variable (rising or declining) star formation histories. However, additional observational constraints, such as emission line strengths discussed later, may be necessary to distinguish between among variable SFHs, since both rising and declining histories yield fairly similar predictions for many physical parameters (cf. dBSS12 and above).

4.2. IR and UV luminosity as tracers of the SFR in light of different star formation histories

Often the IR luminosity or the sum of IR+UV luminosity are used as a measure of the star formation rate (e.g. Kennicutt 1998; Burgarella et al. 2005; Reddy et al. 2010, 2012b). The basic assumptions made for such conversions are that of constant SFR and an age of  ≳50−100 Myr, necessary to reach an equilibrium output for the UV and bolometric luminosities (see also Fig. 2). The determination of the SFR from L_IR also assumes that the total (stellar) bolometric luminosity is reemitted in the infrared domain, which requires significant dust attenuation. Obviously the above assumptions cannot be valid in general, and in particular they are not compatible with the different star formation histories explored in this paper.

In Figs. 14 and 15 we show the predicted IR and IR + UV luminosity respectively as a function of the (current) SFR determined from fitting our models to the sample of z ~ 4 galaxies. The “standard” SFR(IR) calibration given by Kennicutt (1998) – based on the bolometric luminosity – is also indicated for comparison. Overall quite large deviations from the simple calibration are obtained for L_IR–SFR, and smaller deviations for (L_IR + L_UV)–SFR. As expected model fits assuming constant SF (model B) agree very well with the Kennicutt relation when the IR+UV luminosity, providing a very good measure of the bolometric luminosity, is used. However, for galaxies with SFR  ≲  20–30 M_⊙ yr^-1, the predicted IR luminosity is decreased with the respect to this relation. This is simply due to the fact that these galaxies are not found to be dusty enough to reprocess a significant fraction of their UV light into IR radiation. This is also clearly illustrated by the large uncertainties in L_IR, which can reach values down to 0 if the pdf of A_V does so. This effect, shown here for SFR = const., is present for all star formation histories.

More generally, for varying star formation histories, both the IR and the IR + UV luminosity predicted by model fits are lower than expected from the Kennicutt relation between SFR and luminosity (Figs. 14 and 15). Two effects explain these differences. First, dust extinction needs to be high enough for L_IR to be a good tracer of the total SFR (cf. above). Second, the ages are younger than the time needed to reach the bolometric output assumed by Kennicutt (1998). Typically, this is the case for ages t/τ ≲ 1 for declining SFHs, somewhat longer for delayed and rising star formation histories, and  ~100 Myr for constant SF, as illustrated in Fig. 2. The second effect explains the values of L_IR + L_UV below the canonical L_IR–SFR relation. Added to that, the first effect further diminishes L_IR for galaxies with lower SFR, as already explained above. In some cases, although for a small fraction of galaxies, we find that the predicted IR or IR+UV luminosity is higher than what is expected from the Kennicutt (1998) relation using the model SFR. This is found with declining or delayed star formation histories for galaxies where t/τ ≫ 1, i.e. where star formation starts to decrease but longer-lived stars still contribute to the bolometric luminosity. These objects are expected to show weak emission lines, i.e. they should include galaxies from the category of “quiescent” LBGs identified by dBSS12, and which can also be recognized by their 3.6 micron excess, if at z ~ 3.8–5.

To illustrate more quantitatively the possible deviations from the “standard” SFR calibration, we show in Fig. 16 the full probability distribution function (pdf) of the difference in IR+UV luminosity, ΔL_IR + UV, between the prediction from our SED fits for the entire B-drop sample and the luminosity for the same SFR following the Kennicutt relation. Negative values of ΔL_IR + UV correspond to cases where the model SFR is higher than what would be expected from the standard calibration, positive values the opposite. Of course, the pdf shows the same offsets already seen from the median values plotted for each galaxy in Fig. 15. The median offsets from the Kennicutt relation are  −0.21 dex with a distribution systematically shifted to lower luminosity (higher SFR) for the rising SFH with variable timescales (model D). For the declining histories (model A) the median offset is only  −0.07 dex, since the distribution extends to both higher and lower L_IR + UV. For model B (SFR = const.) a small shift (−0.14 dex) is also found, since L_UV does not completely measure the total bolometric luminosity for galaxies with low SFR and low dust content. More important, although the median offsets are relatively small, the models with variable SFHs predict deviations which can be quite high. For example, deviations between  − 0.83 (−0.66) and  −0.04 (0.46) are found within 68% confidence for model C (A).

As already discussed, part of the deviation is due to relatively young ages found. Indeed, imposing an age prior (t > 50 Myr, e.g.), reduces the median offset log ΔL_IR + UV to –0.15 and diminishes the scatter (68% confidence range now from –0.38 to 0) for rising SFHs, and to a median of –0.03 (68% confidence from –0.12 to 0.80) for exponentially declining histories. If only the UV or the IR luminosity, not both, are available, the deviations and scatter are obviously even larger than shown in Fig. 16, as already seen in Fig. 14. For comparison, Wilkins et al. (2012) have recently discussed the SFR(UV) relation expected at high redshift from their semi-analytical models of galaxy formation and evolution. They find essentially no deviation from the Kennicutt calibration (when the same IMF is used) and a scatter of  ~20%. These differences compared to our results are most likely due to the fact that semi-analytical models may not resolve relatively short time-scales and young ages for star-formation.

Of course it should be noted that the SFR shown in the above plots is not a direct observable, since it corresponds to the current SFR determined from the broad-band SED fits, and since the SFR derived in this manner depends itself on model assumptions, such as the assumed star formation history (cf. Sect. 3 and Fig. 7). Observationally the best measure of the current SFR obtained from the SED fits should be the one measured by hydrogen recombination lines (e.g. SFR(Hα)), once corrected for dust attenuation, since these emission lines trace best short-lived massive stars. From the Hα flux predicted by our SED models, the consistency with SFR(Hα) can indeed by verified.

The main conclusions to be drawn from Figs. 14 and 15 is that IR or (IR+UV) SFR calibrations assuming constant star formation over long timescales should be used with caution for the interpretation of high-z star forming galaxies (e.g. LBGs), since they are susceptible to variations of the star formation history on short timescales, for which there are good indications (see dBSS12 and this paper). In any case, the effects discussed here need to be taken into account to properly predict the IR emission of LBGs and for comparisons with observations. Possible problems related to the use of inconsistent SFR indicators, and implications on the SFR-mass relation are discussed below (Sects. 4.3 and 7.2).

4.3. Implications on the SFR-mass relation

Figure 17 shows the SFR-mass relation with SFR transformed to direct observables, namely the sum of the UV and IR luminosity. Compared to the SFR – M_⋆ relation plotted in Fig. 9, which shows a very large scatter (with an upper bound up to 2 dex above the mean relation at z ~ 2) for models allowing for young ages and for variable star formation histories, the scatter in this relation is now significantly reduced, approximately by one order of magnitude, as judged by the upper boundary. The reason is again due to the fact that the models with the highest current SFR/M_⋆ values correspond to relatively young ages, where the bolometric output per SFR is lower (Fig. 2). Hence the smaller dynamic range in bolometric luminosity, which is well approximated by L_IR + L_UV. In other words, although the intrinsic, current SFR may show a wide dispersion, the SFR one would infer from the UV + IR luminosity assuming “standard” SFR conversion factors would show a much smaller dispersion! This demonstrates that the finding of a small dispersion in the bolometric luminosity – M_⋆ plane does not exclude a high dispersion between current SFR and stellar mass, and objects with very high values of the specific star formation rate, SFR/M_⋆.

4.4. Influence of the minimum age on the dust attenuation and IR luminosity

How does the choice of a minimum age affect the derived attenuation and hence the predicted IR luminosity? Figure 18 shows the result of this exercise for the exponentially declining and rising SFHs (model A and D). We find that introducing a lower limit t_min (age prior) leaves the predicted IR luminosity of LBGs unchanged on average for the rising star formation histories, whereas L_IR diminishes for declining SFHs. The basic reason for this behavior is the following: since for rising SFHs the intrinsic UV emission does not decrease when the age is increased, the attenuation and hence their IR luminosity remains comparable for most galaxies. In contrast, for decreasing star formation histories an older age allows one to reproduce the SED with a lower attenuation (cf. the well known age-extinction degeneracy). Thus a lower L_IR, comparable to the values predicted for the constant SF model with the same minimum age, is expected for this case. In other words, although the predicted L_IR depends in general on assumptions on the minimum age of stellar populations, this is not the case for rising star formation histories, which consistently predict a higher UV attenuation and hence IR luminosities than constant SF or declining histories. A similar distinction was already seen in the SFR-mass relation (Fig. 10). Measuring L_IR can help test/constrain these different scenarios.

5. Predictions for ALMA

5.1. Continuum flux predictions

In Fig. 19 we show predicted fluxes for the B-drop samples in one of the ALMA bands (band 7 covering  ~0.8–1.1 mm) available during the early ALMA cycles. As explained in Sect. 2 these fluxes are derived from the predicted IR emission, L_IR discussed above, with the assumption of modified black bodies with a typical dust temperature of T_d = 35 K and β = 2. Changes of the dust temperature by  ± 10 K correspond to a change of the flux by a factor  ~2.5–2.8 approximately. At this redshift (z ~ 4) the emission increases (decreases) with decreasing (increasing) dust temperature in this band, since the peak of the spectrum lies at shorter wavelengths than band 7 considered here.

Figure 19 clearly shows the significant flux differences predicted by models assuming different star formation histories. If the dust properties (temperature and exponent β) do not change, the IR-mm fluxes simply reflect the differences in IR luminosity predicted by the models, as already discussed in Sect. 4. While the majority of the z ~ 4 LBGs show fluxes below the sensitivity limit reached with short (60 s) exposures with the full ALMA array, approximately half of the sample should be detectable with 30 min exposures. In any case, already with the current, incomplete ALMA array the brighter part of the distribution should be detectable. Such measurements will provide badly needed, direct measurements of the UV attenuation of z ~ 4 LBGs and very useful tests/constraints on the SED models and on the star formation histories of these distant galaxies.

When selecting galaxies by their IR luminosity, which part of the SFR-mass diagram such as Fig. 9 do we probe? Naively one would expect this to translate to a simple limit in SFR. However, since SFR/L_IR is a priori not constant, i.e. age dependent, and since L_IR also depends on the amount of dust attenuation, the limit is more complicated. This domain is e.g. illustrated for L_IR > 2 × 10¹¹ L_⊙ by the dashed line in Fig. 9. Whereas this criterium selects all massive galaxies with SFR  ≳ 35 M_⊙ yr^-1, such galaxies with M ≲ 6. × 10⁸ M_⊙ have relatively young populations and less dust, which implies a lower IR luminosity for the same SFR, below the IR luminosity limit. Although the detailed shape of curves of L_IR = const. in Fig. 9 depend somewhat on the value of L_IR, their qualitative behaviour is similar.

5.2. Evolution of the IR luminosity with redshift

The typical predicted IR luminosities at z ~ 4 are shown in Fig. 20 as a function of the absolute UV magnitude (here M₁₅₀₀ for comparison with dBSS12) for model A, i.e. declining star formation histories. The figure shows the 2D probability distribution function (pdf), the median values for each galaxy (white circles), and the median of the pdf in bins of UV magnitude (red stars). As mentioned earlier, the pdf also includes a tail to very low IR luminosities, set somewhat arbitrarily to L_IR values outside the plot, since some of the SED fits also allow for zero extinction. This explains why the median IR luminosity (red stars) is shifted towards lower values than the peak of the pdf. Incidentally, at z ~ 4, the median L_IR for model A is found to be quite close to the IR luminosity predicted from the UV magnitude using the standard SFR(UV) and SFR(IR) relations of Kennicutt (1998), and assuming SFR(UV)  =  SFR(IR), as shown by the dotted line. On average we find that the UV brightest galaxies are also expected to be brightest in the IR, in agreement with Reddy et al. (2012b); González et al. (2012) and other studies. Typical median luminosities are of the order of L_IR ~ 10¹⁰⁻¹¹ L_⊙ for model A (cf. Fig. 19).

In Fig. 21 we show the IR luminosity for the i-drop sample (z ~ 6) predicted from model A. It should be noted that a higher fraction of the 60 objects of this sample are found with a median attenuation of zero (below our minimum non-zero value of A_V = 0.05), implying thus a median L_IR outside this plot. In other words the complete pdf has a secondary peak at low values of L_IR, not shown here, which explains – as for Fig. 20 – the shift of the median of the sample (red stars) with respect to the apparent peak of the pdf. Clearly, the predicted IR luminosity of z ~ 6 LBGs is lower by a factor  ~3–4 at most UV magnitudes, due to the lower dust attenuation at high redshift (see dBSS12 for more details). However, despite the lower reddening, the SED fits show that galaxies at z ≥ 6 with dust exist quite likely, even at these high redshifts (cf. dBSS12, also Schaerer & de Barros 2010). Indeed, these studies show a clear trend of increasing average dust attenuation with galaxy mass, at all redshifts.

6. Predicted (rest-frame) optical emission lines

Since one of the main features of our SED analysis is the treatment of nebular lines, whose effects turns out to be not negligible for approximately two thirds of the LBGs (cf. dBSS12), it is obviously of interest to examine the predicted strength of these emission lines, and to test these predictions observationally. From the 63 lines from the (rest-frame) UV, optical and near-IR domain included in our models, we illustrate here Hα, among the strongest optical lines, and [O ii] λ3727, the bluest optical line. Their predicted equivalent width distribution for the B-drop sample and for the different star formation histories are shown in Figs. 22 and 23.

The predicted Hα equivalent widths cover a wide range, with a maximum reaching up to  ~3000 Å for some rare cases of very young galaxies. Most SFHs, except for the model with constant SFR, predict a fairly wide distribution of equivalent widths. The rising star formation histories (models C and D) predict higher median equivalent widths, since the current star formation always produces Lyman continuum photons, hence also emission lines. As expected, models with constant SFR predict a relatively narrow range of equivalent widths, limited by a minimum W_Hα  ≳ 100 Å, which is reached close to the maximum age of  ≲ 1 Gyr. Only models allowing for some decline of star formation (e.g. model A, E), can predict objects with very weak or no emission lines, as clearly shown in Figs. 22 and 23. From these figures is evident that measurements of the equivalent width distribution in LBG samples can, in principle, distinguish different star formation histories. In practice, however, such line measurements are currently difficult or not yet possible, due to limited sensitivity/access at the observed near-IR wavelengths. For this task, recombination lines such as Hα should ideally be used, as they trace directly current star formation, whereas metal lines depend more strongly on nebular conditions (metallicity etc.).

The Hα fluxes predicted for z ~ 4 sample are between 10^-18 and 10^-16 erg s^-1 cm^-2, with the lowest median predicted for the SFR = const. model. The predicted [O ii] λ3727 fluxes are in a similar range, but less differences are expected between different SFHs, as already mentioned above.

Not surprisingly, we find a relatively strong (anti-)correlation of the Hα equivalent width with galaxy mass, shown in Fig. 24. Obviously different models occupy somewhat different regions of this plot, but in all cases do we find such a behavior. Indeed, this trend is expected, since the average age decreases and the specific star formation rate (sSFR = SFR/M_⋆) increases with decreasing galaxy mass for all SFHs (dBSS12). Remember, however, that selection effects and biases can affect such relations (cf. Stringer et al. 2011; Reddy et al. 2012b). W_Hα decreases also with galaxy age, but the spread in age can be very large, given various star formation timescales. The best correlation is indeed found between the Hα equivalent width and the specific SFR, shown in Fig. 25. Again, this is a natural trend, since the Hα equivalent width measures the ratio between the Hα flux – tracking the current SFR – and the underlying continuum flux, which depends on the accumulated stellar mass. We note that all models predict a very similar relation, which also holds for all samples from z ~ 3 to 6, although they occupy different regions of the relation. Measurements of W_Hα (from spectroscopy or estimates from broad-band photometry) may therefore be used to determine the sSFR of individual high-z galaxies. As Fig. 25 also shows, metal lines such as [O ii] λ3727, which are not recombination lines, trace the sSFR less well.

Our predicted equivalent widths are comparable to those found recently in emission line galaxies at z ~ 0.6–2.4 for [O ii] λ3727, [O iii] λλ4959,5007, and Hα (Straughn et al. 2009; Atek et al. 2011; van der Wel et al. 2011) An anti-correlation of W_Hα with galaxy mass is found at all redshifts, where data is currently available, from z = 0 to 2–2.6, as recently discussed by Fumagalli et al. (2012). Furthermore, the predicted Hα equivalent widths of the B-drops are in good agreement with the observations of z ~  2–2.6 LBGs by Erb et al. (2006); Mancini et al. (2011) in the same mass range (typically at M_⋆ ≳ 10¹⁰ M_⊙). Our models thus predict that this trend of increasing Hα equivalent width continues down to LBGs of lower mass.

The 47 strong Hα emitters identified by Shim et al. (2011) at z ~ 4 are included in our LBG sample, and amply discussed in dBSS12. Their selection, based on an excess in the 3.6 μm filter with respect to 4.5 μm, is sensitive to EW(Hα)  ≳  200 Å. Clearly, our method is able to detect/predict weaker lines, although of course with increasing uncertainty (not shown here). The highest inferred equivalent widths from Shim et al. (2011), well in excess of 1000 Å, are comparable to our values, as are the predicted Hα fluxes. Although very rare at low redshift, galaxies with EW(Hα)  >  500 Å also exist in nearby galaxies, as pointed out by Shim & Chary (2012). They are all found at masses well below log M_⋆ ≲ 9.5 M_⊙ and at high sSFR, quite compatible with our model predictions. Finally our predictions are also in agreement with the recent determination of Hα equivalent widths for z ~ 3.8–5 LBGs with spectroscopic redshifts from Stark et al. (2012), who find a mean value of log  EW(Hα) = 2.57 ± 0.25 Å.

We conclude that the emission line strengths predicted by our models from fits to the observed broad-band photometry and the correlations found are compatible with the currently existing data and trends. Future observations of emission lines in z ~ 3 − 6 LBGs with existing and new ground-based facilities, and with the James Webb Space Telescope will hopefully provide further tests on the importance of emission lines and constraints on the star formation histories and specific star formation rate of these galaxies.

7. Discussion

7.1. Influence of different extinction laws

All our SED fits have been carried out using the Calzetti attenuation law, commonly thought to be appropriate for starburst galaxies. To examine the effect of other laws we have rerun all our fits for the B-drop sample with the steeper attenuation/extinction law for the SMC Prevot et al. (1984); Bouchet et al. (1985). Other laws, intermediate between the relatively grey Calzetti law and the steep SMC curve, should yield intermediate results.

Overall, we find a lower fit quality (higher reduced ) for the vast majority of objects using the SMC law. However, for most of them the difference remains relatively small (). We verified that no systematic effect is found as a function of galaxy age, as could e.g. be expected if young galaxies (e.g. t < 100 Myr, as judged by the age derived assuming the Calzetti law) would be better fit with the SMC law, as suggested e.g. by Reddy et al. (2010, 2012b).

Since the SMC law is steeper than that of Calzetti, SED fits for LBGs are generally yield a lower UV attenuation f_UV (and even more so a lower A_V), older ages, and hence on average somewhat larger stellar masses (e.g. Papovich et al. 2001; Verma et al. 2007; Yabe et al. 2009; Reddy et al. 2010). As expected, our fits including nebular lines also yield the same average trends. The young ages favored with the Calzetti law for some objects by the strong Hα excess at 3.6 μm are in part compensated by increasing the star formation timescale τ for declining SFHs. The resulting Hα equivalent width distribution (cf. Fig. 22) remains quite similar, with the largest differences due to the assumed star formation histories.

Important differences of the SED fits with the SMC law are in particular a lower SFR and less IR emission. The former is due both to the older age (hence lower SFR/UV, cf. Fig. 2) and to the lower UV attenuation. The lower attenuation implies immediately a lower IR luminosity. For model A (declining SFHs) the median SFR of our B-drop sample is lowered by a factor  ~2, from 11.6 to 5.4 M_⊙ yr^-1, the median IR luminosity from 4. × 10¹⁰L_⊙ by a factor  ~5, compared to SED fits with the Calzetti law. The effect of adopting the SMC law on the SFR-mass on the (L_IR+L_UV)–mass relations for model A are shown in Fig. 26. The top panel clearly shows the reduced IR luminosity, implying a lack of B-drop galaxies with L_IR ≳ 2. × 10¹¹ L_⊙ (in this sample) if the SMC law was applicable. A large scatter is still found in the SFR-mass relation with some objects reaching large specific SFRs, as shown in the bottom panel. However, at masses  ≳ 10⁹ M_⊙ the SED fits with the SMC law predict basically no LBGs with high specific SFR, i.e. above the z ~ 2 relation of Daddi et al. (2007). This shift to lower SFR values is again mostly due to age shifts, leading to ages  ≳ 100 Myr for all galaxies with M ≳ 10⁹M_⊙. It must, however, be remembered that the model SFR shown here represents the current SFR, which shows a wider dynamic range (spread) than the observable UV and IR luminosity (cf. Sect. 4.3).

At low redshift, and ocasionally at z ~  2–3, observations of the IR/UV ratio and measurements of the UV slope β have been used to distinguish the attenuation law (e.g. Baker et al. 2001; Siana et al. 2008, 2009; Reddy et al. 2010). The so-called IRX–β plot, based on the predicted IR luminosity plus the observed L_UV and β, for our z ~ 4 LBG sample and for the two attenuation laws is shown in Fig. 27. As already mentioned above, a lower IR/UV ratio is found from SED fits using the SMC law. Besides an accurate measurement of the IR/UV ratio (soon feasible with ALMA), this method requires also a good knowledge of the intrinsic UV slope, and accurate measurements of β to constrain the attenuation law. Indeed, as clear from this plot, UV slope measurements of individual galaxies show a large scatter and their errors are relatively large since the β slope depends on one UV color, here β ∝ 5.3 × (i − z) color, following Bouwens et al. (2009). More accurate measurements are feasible, e.g. using more photometric bands or based on SED fits (Finkelstein et al. 2012; Castellano et al. 2012). In any case, the use of the IRX–β plot to distinguish attenuation laws also requires the knowledge of the intrinsic UV slope, which may not be unique, since the intrinsic UV slope depends on the star formation history, age, and metallicity of the stellar population.

Distinguishing different attenuation laws may thus not be straightforward. However, some combinations of star formation history and attenuation law are clearly distinguishable once IR luminosities are measured with ALMA. As already mentioned, our SED fits to the present data yield in general better fits with the Calzetti law for the vast majority of objects, although the difference in is not very large.

7.2. Comparison with earlier studies

We shall now briefly discuss other papers, which have presented tests for or constraints on different star formation histories for Lyman break galaxies using also the IR as a constraint.

To the best of our knowledge, our study is the first discussing SED fits with various star formation histories including both declining and rising parametrisations for a large number of z ~ 3−6 LBGs. At somewhat lower redshift Reddy et al. (2012b) have recently examined a sample of  ~300 LBGs with MIPS 24 μm detections and spectroscopic redshifts 1.5 ≤ z ≤ 2.6. The 24 μm flux is used to infer the IR luminosity, which in turn provides an independent constraint on the SFR and dust content (as discussed above and in earlier papers, cf. Reddy et al. 2008). From their work Reddy et al. (2012b) conclude that rising star formation histories provide a better agreement between different SFR indicators, and hence that declining star formation histories may not be accurate for typical galaxies at z ≳ 2.

This conclusion appears in contrast to ours, since we find that declining SFHs provide the best fits for the majority of z ~ 3−6 LBGs, and that delayed and rising star formation histories fare comparably. First, we have to note that based on the SED fit quality (χ²) alone Reddy et al. (2012b) also find quite similar results between their different star formation histories. They then use additional constraints (mid-IR observations) to examine SFR indicators for their consistency. The principle of their method is exactly what we’re promoting to use. However, such observations are not yet available for LBGs at higher redshift. We therefore think that it is premature to draw general conclusions on the star formation histories of LBGs at z ≥ 3 from the dataset of bright z ~ 2 galaxies studied by Reddy et al. (2012b).

We note that in our analysis of z ~ 3−6 LBGs with declining star formation histories we find values of t/τ < 1 for the majority of objects, in agreement with the observational arguments mentioned by Reddy et al. (2012b), and in contrast to their finding of t/τ > 1 for the majority of z ~ 2 sample. Furthermore, all the currently available constraints for z ~ 3−6 LBGs can be well reproduced with variable (i.e. non-constant) SFHs, whether exponentially declining or rising (cf. also dBSS12). In any case, stellar populations/SFHs of LBGs at higher redshift may be different, and from the result of Reddy et al. (2012b) at z ~ 1.6−2.4 one cannot generalize that exponentially declining SF histories are not appropriate for LBGs at z ≥ 3, since additional observational constraints – e.g. from the IR – are lacking at these redshifts. It is therefore important to obtain direct measurements of the IR luminosity at higher redshifts to try distinguish different star formation histories of high-z galaxies.

Another issue concerns the small subsample of “young” LBGs discussed by Reddy et al. (2012b). Indeed, for galaxies classified as “young” (t < 100 Myr) according to their SED fits with the Calzetti law, Reddy et al. (2012b) find that the SFR inferred from the fits, SFR(SED), is systematically larger than the SFR(IR+UV) using the standard calibrations of Kennicutt (1998). From this they conclude that there is an apparent conflict with young ages, which they propose to resolve by invoking a different, steeper attenuation law for young objects. We are not entirely convinced by the reality of this conflict, and hence by the proposed solution. Indeed, although they note that a young age is not compatible with assumptions made for the standard UV–SFR conversion, they neglect the fact that both IR and UV conversions depend on the SFH and age in a similar way (see their Figs. 25 and 2), modifying the (UV + IR) – SFR conversion upward by typically twice the amount they state. The effect they find is clearly shown to be predicted for our z ~ 4 sample, as illustrated in Fig. 15, where we note that for the majority of objects fitted with unconstrained ages and found to be younger than 100 Myr from our models, the SFR(SED) is larger than the SFR(IR+UV) obtained from an inconsistent application of the standard calibrations. Offsets by up to a factor  ~10 from the standard Kennicutt relation can be found (cf. Fig. 16). Whether this can quantitatively solve the (sometimes large) discrepancies found by Reddy et al. (2012b) remains, however, to be seen.

Other comparisons with Reddy et al. (2012b) regarding constant and rising star formation histories have already been discussed in dBSS12, and shall not be repeated here.

7.3. Derived ages and implications for the galaxy populations

At first sight the ages derived from our models and the non-evolution of the median mass at a given UV magnitude (cf. dBSS12) imply that a fraction/most of the LBGs observed in a given sample (e.g. at z ~ 4) would not have been seen at earlier epochs (e.g. in the z ~ 5 and 6 samples). This seems e.g. opposite to the finding of Papovich et al. (2011), who argue that LBGs are visible at multiple epochs. Furthermore our young ages and age spread are also related to the scatter we find in the SFR-mass relation and to the relatively high sSFR we find at high redshift (dBSS12), with important implications, if true. This raises two immediate questions: first, are our derived ages correct and physically plausible? And second, can our results be reconciled with the picture of LBG populations with star formation histories rising on average?

It must be recognized that stellar population ages and uncertainties are strongly model dependent, and that both the age uncertainties for individual galaxies and the age spread around the median for a galaxy sample are fairly large. Ages obtained with models including nebular emission are generally younger than without (cf. dBSS12). For the majority of galaxies the derived age is found to be higher than the dynamical time, and our models yield galaxy ages which increase both with galaxy mass and with decreasing redshift, as expected from the progressive build-up of galaxies. For a fraction of up to 20% of the LBGs our models give ages shorter than the typical dynamical time, estimated between 40 Myr at z ~ 3 and  ~20 Myr at z ~ 6 (cf. dBSS12), which may be unphysical. Imposing longer star formation timescales (larger τ) allows one to reduce this fraction. However, a fraction of galaxies in all our samples is best fit with young ages. Incidentally, these galaxies turn out to be among the least massive, which at z ~ 4 are characterized by a very blue (3.6 − 4.5) micron color, indicative of strong Hα emission. Besides this potential problem with young ages for a small ( ≳ 20%) fraction of LBGs, we should be aware that the SED fits may be mostly sensitive to the most recent SF event, due to the well known “outshining” effect (cf. Wuyts et al. 2009; Maraston et al. 2010). For this reason age and also the galaxy mass may be underestimated, although it remains difficult to establish by how much.

Our finding of relatively short SF timescales and young ages, and the non-evolution of the mass–UV magnitude relation from z ~ 5 to 3 (dBSS12), favours episodic star formation, in agreement with the findings of Stark et al. (2009). If successive SF events take e.g. place with increasing strength and emission from the most recent event dominates the SEDs at most wavelengths, this scenario could be reconciled with the picture of slowly rising star formation histories of LBGs inferred from the evolution of the UV luminosity function by Papovich et al. (2011). Whether this works out quantitatively remains, however, to be seen. In any case the observational tests proposed here, using dust emission in the IR and optical emission lines, will also help to answer these questions.

7.4. Possible caveats

Before closing this section, we wish to remind the reader that our SED fitting tool suffers from the same possible limitations as other models, including e.g. the unknown star formation histories of galaxies (modeled here with simple parametrisation), the question of multiple populations and the “outshining” problem due to young populations, uncertainties related to the extinction (which law? a unique attenuation for all stars? etc.), uncertainties in the stellar initial mass function, and of course also uncertainties in stellar evolution and atmosphere models. Some of these issues are e.g. addressed by Maraston et al. (2010); Wuyts et al. (2009); Charlot & Longhetti (2001); Walcher et al. (2011); Levesque et al. (2012); Schaerer (2012). The reader should also be aware that the timescale of UV and IR emission, which is fairly similar in our models and in the “standard” star formation rate calibrations, could be different in more complex situations, e.g. in cases of time-dependent attenuation, multiple populations or alike (e.g. the models of da Cunha et al. 2008). However, such models include additional free parameters, and have, to the best of our knowledge, not yet been used to fit distant galaxies. For simplicity we have therefore fitted the SEDs with single stellar populations and a single extinction/attenuation law.

In this respect our results should of course be considered as “differential”, i.e. with respect to other models making the same/similar assumptions, but neglecting the effects of nebular emission and of variable star formation histories, which are explored here systematically for the first time for a large sample of Lyman break galaxies. Concerning nebular emission our models make some simple assumptions, which may also influence the results. Our main assumptions are no loss of ionizing photons (i.e. no Lyman continuum escape, no dust inside the Hii regions etc.), and the same attenuation for the stellar continuum and nebular lines. Both assumptions maximize to some extent the possible contribution of nebular emission, and could therefore lead to an overestimate of the effects of nebular emission on the derived physical parameters. The results would then be intermediate between the present case with nebular emission and the models without, discussed in depth in dBSS12. Exploring variations of the above assumptions may be interesting and worth pursuing in future studies.

In any case, it is clear that for a majority (approximately two thirds) of LBGs the inclusion of nebular effects improves the quality of the SED fits, and that models without nebular lines cannot reproduce the observed excess in certain photometric bands (most prominently the 3.6 μm excess for galaxies with redshift between 3.8 and 5). For the bulk of the LBGs our models should therefore provide a better description than previous models neglecting nebular emission, despite the simplifying assumptions made. To progress further on these questions additional observational constraints, such as direct measurements of IR luminosities, medium-or narrow-band photometry, or measurements of emission lines should be of considerable help. The predictions presented here should also provide a useful base for future comparisons and tests of the SED models. In parallel, our models will also be tested against observations at lower redshift where IR and emission line measurements are available for some galaxies.

7.5. IR-to-mm observations to determine the UV attenuation in LBGs and constrain their SFHs

One of main objectives of this paper has been to show how the assumption of different star formation histories for SED fits of Lyman break galaxies at z ≳ 3 affects their derived physical parameters, and in particular their inferred dust extinction. Independent, observational data is now needed to determine dust extinction. Three “classical” measurements are normally used to do this: the UV slope, a H recombination line ratio (e.g. the so-called Balmer-decrement for optical lines), or the IR/UV ratio. Since the UV slope has already been used for our SED models, and H emission lines beyond Lyα are currently difficult or impossible to obtain for z > 2.3, measurements of the IR/UV ratio are needed.

The deepest IR/sub-mm images available to date, taken with the Herschel satellite in the GOODS-S field, have detected few individual, very bright (L_IR > 10¹² L_⊙) LBGs out to z ~ 2.5−3 (e.g. Elbaz et al. 2011; Burgarella et al. 2011; Reddy et al. 2012b). Stacking techniques have so far been necessary to detect fainter, more typical galaxies, and to try extend these studies to higher redshift (Magdis et al. 2010; Rigopoulou et al. 2010; Reddy et al. 2012a). At z ~ 4, for example, the expected IR luminosity of a typical LBG with is L_IR < 10¹¹ L_⊙ (cf. Fig. 20), clearly beyond the reach of Herschel. However, detections of individual LBGs at these redshift are now becoming possible with ALMA. Waiting for such new observations we have recently analyzed stacks of the deepest IR/sub-mm data for the GOODS-S field, to examine whether some of the LBGs studied here are detected or whether they are detectable in stacks, and whether all of the star formation histories are compatible with the present data. A publication on this work is in preparation.

Other more indirect tests of the different IR luminosities predicted by the models discussed here could be carried out using IR galaxy counts and IR background measurement. Such tests will be carried out in the future, but are, however, clearly beyond the scope of the present work.

8. Summary and conclusions

Following up on our earlier detailed study (de Barros et al. 2012) of a large sample of LBGs from redshift z ~  3 to 6 located in the GOODS-South field, using for the first time an SED fitting tool including the effects of nebular emission on the synthetic photometry, we have examined the impact of different star formation histories (SFHs) on the derived physical parameters of these galaxies, on the SFR-mass relation, on different SFR indicators (UV and IR), on the expected dust extinction and the corresponding IR luminosity, and on the expected strengths of emission lines such as Hα and [O ii] λ3727.

To do so, we have carried out SED fits for five different SFHs including exponentially rising and so-called delayed SFHs, plus the three histories already considered in dBSS12 (see Table 1 and Fig. 1). Metallicity is also treated as a free parameter, and we have examined the effect of two different extinction/attenuation laws (Calzetti and SMC). The usual physical parameters, stellar mass, SFR, age, and attenuation are derived from the SED fits to the broad band photometry reaching from the U band to 8 μm, using Monte Carlo simulations to derive their median values (and the detailed probability distribution function, generally not discussed here). We have also computed consistently the predicted IR luminosities, L_IR, for all galaxies, assuming energy-conservation, i.e. that all the radiation absorbed by dust is reemitted in the IR. Finally, the predicted IR luminosities have been translated to flux predictions in various IR bands, assuming modified black body spectra. The L_IR predictions allow us in particular to examine in a consistent way the effects of variable SFHs and ages on this observable quantity, showing thus significant departures from results assuming inconsistent SFR(IR) or SFR(UV) calibrations.

Our main results concerning the impact of star formation histories on the physical parameters of LBGs, exemplified to a sample of 705 LBGs at z ~ 4 (B-drop galaxies), can be summarized as follows (see Sect. 3):

• Compared to commonly adopted SED fits assuming constantSFR, no nebular emission and an age prior of t > 50 Myr, models with exponentially declining SFHs, nebular lines and no age constraint yield younger ages, lower stellar masses, higher current SFR, higher specific star formation rates (sSFR = SFR/M_⋆), and higher dust extinction (A_V), as already shown in de Barros et al. (2012). Exponentially declining SFHs yield overall the best fits (in terms of χ² ) for the majority of LBGs. Based on the available SED constraints it is, however, difficult to distinguish different SFHs, although various arguments (in particular the distribution of emission line strengths inferred from broad-band photometry for z ~ 3.8–5 galaxies) favor clearly variable, i.e. non constant, histories (cf. dBSS12).

• Assuming delayed star formation histories one obtains basically identical physical parameters (and fit qualities) as for exponentially declining SFHs.

• Rising star formation histories with variable timescales imply generally a similar stellar masses, and comparable or somewhat higher dust extinction than models assuming declining SFHs. The latter leads to the highest star formation rates and to similar or higher IR luminosities as for declining histories.

• Overall “standard” models assuming constant and neglecting lines predict systematically higher stellar masses, lower extinction, lower SFR, lower IR luminosities, and more narrow range of equivalent widths for optical emission lines than all the other star formation histories considered here.

Combining these different physical parameters we obtain the following (Sect. 4):

• We find significant deviations between the derived SFR and IR luminosity from the commonly used SFR(IR) or SFR(IR+UV) calibration of Kennicutt (1998). Such differences naturally arise, due to differences in the derived ages and in the adopted star formation histories. In most cases (i.e. for most galaxies and SFHs) we find that the Kennicutt relation will underestimate the true, current SFR derived from the SED fits (Sect. 4.2). Consistent SED studies including also the IR are therefore necessary, if the SFHs and ages may differ from those assumed in standard SFR calibrations.

• A large scatter is found in the SFR-mass relation for models with declining and delayed SFH and no age prior (cf. de Barros et al. 2012). The same also hold for models with rising star formation histories. The scatter is reduced when a minimum age (e.g. t > 50 Myr) is adopted. Even in this case, models with rising SFHs can show a large scatter, since high SFRs are found to the high(er) extinction.

• A large scatter in the SFR-mass relation does not necessarily imply the same scatter in the L_IR (or (L_IR + L_UV)–mass relation and vice versa, when the IR luminosity is computed consistently from the chosen SED model (i.e. accounting for age and SFH effects). The same also applies to SFR(UV)–mass diagrams where the UV luminosity and the corresponding standard SFR conversion is used. We suggest that the true scatter in the SFR-mass relation obtained in this way may indeed be underestimated, if the true star formation histories are variable on relatively short timescales. Indeed such SFH variations can reproduce more successfully features related to emission lines, such as the observed 3.6 μm excess in z ~  4–5 LBGs. They may also be more relevant to at higher redshift, where the dynamical timescales decrease with (1 + z)^−3/2.

Our consistent predictions of IR luminosities (and fluxes) show that different SFHs lead to significantly different amounts of reddening and hence to different IR/UV luminosity ratios. Measurements of IR luminosities of individual LBGs or statistical samples of such galaxies can be used to distinguish different SFHs, and hence also different specific SFRs predicted by such models. ALMA observations will thus be able to provide independent constraints on behavior of the sSFR at high redshift, and on the scatter in the SFR-mass relation.

We show predictions for the IR luminosities of B-drop and i-drop galaxies for different star formation histories and as a function of UV magnitude. The typical/median L_IR is predicted to be L_IR ~ 10^{10 ...11} L_⊙ for LBGs with absolute UV magnitudes of M_UV ~  −22 to  −19.

Finally we also show the predicted strengths (equivalent widths) of the Hα and [O ii] λ3727 emission lines. Again, different star formation histories naturally lead to different EW distributions, which can in principle be used to constrain the SFHs. Our models predict on average higher equivalent widths in low mass galaxies, in agreement with currently available observations at z < 3, and a clear anti-correlation of W_Hα with the specific SFR.

Our predictions should in particular provide new tests using IR observations with ALMA and/or measurements of (rest-frame) optical emission lines to obtain a better insight on the star formation histories of high redshift LBGs, on the behaviour of the SFR-mass relations and on the evolution of the specific SFR with redshift.

---

Acknowledgments

We thank Mirka Dessauges-Zavadsky, Thomas Greve, and Michel Zamojski, for discussions and comments on the manuscript, and numerous other colleagues for stimulating and critical discussions during the last two years. We also thank the referee for useful comments which helped to improve the paper. This work is supported by the Swiss National Science Foundation.

References

All Tables

All Figures