© ESO, 2012

1. Introduction

It has been nearly 20 years since the first surveys explicitly dedicated to detect galaxies at z ≃ 3 were completed (Steidel & Hamilton 1992, 1993; Giavalisco et al. 1994; Steidel et al. 1995). In the early 90 s, these seminal works have introduced the so-called “Lyman break” method to select rest-frame UV bright objects, i.e. star-forming, moderately absorbed galaxies, at z ≃ 3 and at higher redshifts. In subsequent years, this technique has been widely applied by surveys covering either progressively wider areas (e.g. Steidel et al. 2003; Ouchi et al. 2004; McLure et al. 2009) to enlarge the volume sampled, reaching fainter magnitudes (e.g. Bouwens et al. 2007; Reddy & Steidel 2009) to detect the faintest galaxies, or extending to longer wavelengths to explore the highest redshifts (e.g. Bouwens et al. 2004; Castellano et al. 2010; Bouwens et al. 2011a). In the meantime, spectroscopic follow-up campaigns of Lyman break galaxy (LBG) samples have shown that this technique is extremely efficient and reliable up to the highest redshifts (e.g. Steidel et al. 1996; Vanzella et al. 2009; Stark et al. 2010; Vanzella et al. 2011; Pentericci et al. 2011).

Some of the primary goals of these surveys have been the statistical analysis of high-redshift galaxies, in particular of their UV luminosity function (LF) from z ≃ 3 to z ≃ 7 and beyond (Steidel et al. 1999; Bouwens et al. 2007; McLure et al. 2009; Bunker et al. 2010; Grazian et al. 2011), and their contribution to the star formation rate density (SFRD).

Understanding the physical properties of UV-selected galaxies is also fundamental to derive the global properties of the Universe at these high redshifts. In particular, to estimate their contribution to the total SFRD one needs to apply the proper conversion between the observed UV rest frame luminosity and the ongoing star formation rate (SFR). This conversion factor depends on the physical properties of the stellar population (IMF, metallicities, and ages) and on the amount of dust extinction. The slope of the UV continuum in LBGs provides precious information for measuring these quantities. Apart from weak absorption lines and other minor features, the spectrum of a star-forming galaxy longward of the Lyα wavelength can be reasonably approximated by a power law until 3200 Å for most combinations of stellar properties and attenuation curves (Calzetti et al. 1994; Meurer et al. 1999; Calzetti et al. 2000; Gallerani et al. 2010). The actual slope is the result of several factors involved in the stellar mixture in the galaxy: the total amount and composition of dust grains, the metallicity Z, the initial mass function (IMF), and the star formation history (SFH). Lower metallicity stars, top heavy IMFs, and young ages all tend to produce bluer spectra than higher metallicities, standard IMFs and aging stellar populations. Nebular continuum emission can also contribute to redden the slope for galaxies. We note that nebular emission is expected only if the escape fraction of UV ionizing photons is small, as currently observed in most of the galaxies at z ≤ 3.5 (Vanzella et al. 2010; Boutsia et al. 2011). However, dust absorption is likely the strongest reddening factor, and the observed slope is often directly converted into dust extinction assuming standard stellar populations.

In all cases, the UV spectrum is approximated well by a simple power law f_λ = λ^β: the accurate determination of the exponent β is the central goal of the present paper. For reference, a dust-free stellar population with solar metallicity and constant star formation rate has β ≃ −2.2. Clearly, it is impossible to constrain all the physical quantities mentioned above from a single parameter. Nevertheless, a robust measure of β is the first mandatory step in constraining these free parameters and deriving the global properties of LBGs.

Despite the large samples collected in past years, only a few estimates of the statistical distribution of β are available today, and no firm conclusion has been reached on its evolution with luminosity and redshift. Previous determinations have been provided by Bouwens et al. (2009), which found a strong decrease in the average UV slope both with rest-frame magnitude and redshift from z ~ 2 to z ~ 6, and by Bouwens et al. (2010b) presenting evidence of very blue UV slopes in their z ~ 7 sample. At z ~ 6–8 Finkelstein et al. (2010) find blue UV slopes consistent with dust-free normal stellar populations. On the other hand, Dunlop et al. (2012) found that UV slopes are consistent with β ≃ − 2 at z ~ 5–7, and argue that very steep values found at high redshift are due to the effects of photometric scatter and uncertainties. Wilkins et al. (2011) reach a similar conclusion from their z ~ 7 sample, but at z ~ 5 find average slopes redder than both Dunlop et al. (2012) and Bouwens et al. (2009). At z ~ 4 Ouchi et al. (2004) find an average E(B − V) ≃ 0.15 and no significant dependence on UV magnitude for L > L_∗ LBGs, while Lee et al. (2011) found a decreasing extinction, from E(B − V) ≃ 0.22 at L ≈ 4−5L_∗ to E(B − V) ≃ 0.07 at L ≈ L_∗.

Discrepancies among different works are probably due to the known limits of LBG selection techniques (e.g. Stanway et al. 2008) and to uncertainties in the measurement of the UV slope, often based on a single colour. These limitations can be overcome by assembling deep near IR observations of large LBG samples, which are required to strengthen colour selections and to fully sample the UV continuum up to ≃ 3000 Å providing a proper estimate of β. For this reason, early results based on optically selected samples are intrinsically prone to stronger systematic effects. Even more important, however, they suffer from a conceptual limitation. For a given dereddened luminosity (or equivalently a given SFR), even a modest dust absorption produce a significant dimming of the galaxy, which may escape detection in optically selected surveys.

In this paper we take advantage of the new IR observations carried out with WFC3 on board HST. We show that these deep observations in the Y, J and H bands provide significant advantages over previous ones in three different aspects:

• They make a much cleaner selection of LBGs possible, compared to the traditional criteria that only use optical bands. This allows us to reject lower redshift interlopers with much higher reliability.

• They yield a robust determination of the UV slope, providing both a large leverage in wavelength and a sizeable number of independent determinations of the flux (of each object) at different wavelengths.

• They allow us to detect galaxies with a larger slope for a given dereddened luminosity. In this respect, the H-band that we use to detect galaxies constitutes the optimal choice since it corresponds to the reddest portion of the UV rest-frame domain (3200 Å rest-frame at z ≃ 4), and significantly redder than the 1500–1800 Å of optically selected samples. Our samples are thus intrinsically more complete in terms of content of reddened galaxies than optically selected ones, although this is currently limited by the current depth of optical observations used for the colour selection.

To fully exploit these potential advantages, we apply a tailored colour selection criterion to our H-band selected catalogue, and measure β with a straightforward linear fit to the observed magnitudes. In principle, similar results could be derived using an SED-fitting approach on the same H-selected sample (e.g. Fontana et al. 2003), where the slope depends on the free parameters mentioned above: IMF, metallicity, age, and E(B − V). We have decided to adopt a simpler method for two reasons. First, as mentioned above, the four free parameters are to a large extent degenerate, even using additional constraints from the Spitzer magnitudes. Our approach determines a single parameter (β) from four observed magnitudes. Further physical interpretation can be derived from the observed values of β alone. In addition, a complete analysis requires extensive Monte Carlo simulations, which are required to take the involved systematic effects into account. Selecting galaxies with a simple colour selection criterion and adopting a single free parameter makes these simulations affordable.

In this paper we apply our method to z ≃ 4 LBGs. In future papers, we plan to extend the analysis to higher redshifts, also using the forthcoming data sets provided by the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011).

Fig. 1 The (B − V) vs. (V − H) colour selection window adopted in the present paper (right) along with the commonly adopted one based on the (V − Z) colour (left). Blue filled squares and lines are for a CB07 LBG model with Z/Z⊙ = 0.2 and Age = 100 Myr at increasing redshift and E(B − V) respectively. Models having solar metallicity and Age = 100 Myr and Age = 300 Myr are represented in magenta (crosses and long-dashed lines), and red (empty squares and short-dashed lines), respectively. Red and green points indicate lower redshift passively evolving and dusty star-forming galaxies, respectively. See main body for details.

The plan of the paper is the following. In Sect. 2 we present our data-set from which we select a sample of B-dropout galaxies, applying a newly defined IR-based colour selection, discussed in Sect. 3. In Sect. 4, we present the estimate of UV slopes and rest-frame magnitudes for our objects, along with the simulations used to evaluate systematic effects and uncertainties. These simulations are exploited in Sect. 5 to constrain the distribution of β at different UV magnitudes, from which a dust-corrected value of the SFRD is estimated in a straightforward way (Sect. 6). In Sect. 7, we discuss how the UV slope of z ~ 4 galaxies in our sample changes as a function of their dereddened UV emission. Finally, a summary is given in Sect. 8.

Throughout the whole paper, we adopt a Calzetti et al. (2000) extinction law, unless otherwise specified. Observed and rest-frame magnitudes are in the AB system, and we adopt the Λ-CDM concordance model (H₀ = 70 km s^-1 Mpc^-1, Ω_M = 0.3, and Ω_Λ = 0.7).

2. Data and photometry

In the present paper we exploit data obtained in the WFC3/IR ERS2 programme (Windhorst et al. 2011) and the ACS and WFC3/IR UDF programmes (Beckwith et al. 2006; Bouwens et al. 2011b), where the data were processed and combined as described in Koekemoer et al. (2011), along with the full public ACS data available in the GOODS-South field (Giavalisco et al. 2004). Full details on image reduction and catalogue building are given in Grazian et al. (2011) and Santini et al. (2012): we summarize the main properties of our dataset here.

The GOODS-ERS WFC3/IR dataset (HST Programme ID 11359) covers ~40 sq. arcmin. in the Y₀₉₈, J₁₂₅, H₁₆₀ filters. It is made of ten contiguous pointings in the GOODS-South field observed during 60 HST orbits, for a total of 4800–5400 s per pointing and filter. The astrometry of WFC3 ERS images was bound to the GOODS-S v2.0 z band (Giavalisco et al. 2004). The limiting magnitudes are Y = 27.3, J = 27.4, and H = 27.4 at 5σ in an area of ~0.11 arcsec² (corresponding to 2FWHM).

The HUDF WFC3/IR dataset (HST Programme ID 11563) is a total of 60 HST orbits in a single pointing of 4.7 sq. arcmin. (e.g. Oesch et al. 2010; Bouwens et al. 2010a) in three broad-band filters (16 orbits in Y₁₀₅, 16 in J₁₂₅, and 28 in H₁₆₀). The HUDF images were reduced in the same way as the ERS ones, and they were astrometrically matched using the ACS HUDF sources. For the Y₁₀₅ filter, we selected a subset of seven visits out of nine that were not badly affected by persistence from bright sources. The magnitude limits of HUDF WFC3 observations are Y = 29.0, J = 29.2, and H = 29.2 total magnitudes for point like sources at 5σ in an aperture of ~0.11 arcsec².

In both fields we used the available BViz ACS public data: the latest V2.0 version of ACS BViz images of GOODS-South released by the STSci (Giavalisco and the GOODS Team, in prep.) for the ERS, and the deep BViz ACS data presented in Beckwith et al. (2006) for HUDF. To ensure precise colours, the ACS bands were smoothed with an appropriate kernel to match the resolution of the WFC3 images.

We performed object detection on the H₁₆₀ band, obtaining 10 946 and 2418 sources for the ERS and HUDF, respectively. Their total magnitudes were computed using MAG_BEST of SExtractor (Bertin & Arnouts 1996) for galaxies more extended than 0.11 arcsec², while circular aperture photometry (diameter equal to two times the WFC3 FWHM) is used for smaller sources. The magnitudes that were computed in a circular aperture were corrected to total by applying an aperture correction of 0.4 mag. Colours in BVIZYJ were measured running SExtractor in dual image mode, using isophotal magnitudes for all the galaxies.

3. A new approach to selecting z  ~  4 LBGs

We define the sample of z ~ 4 galaxies by applying the Lyman-break technique. The main spectral feature that enables the identification of high-z galaxies is the drop shortwards of the Lyman-α, which, at z ≳ 3.5, is redshifted at λ ≳ 5500 Å. The Lyman break can thus be effectively sampled by the B and V bands (B-dropout galaxies). The usual colour-selection of z ~ 4 LBGs is carried out on a (B − V) versus (V − Z) diagram (e.g. Giavalisco et al. 2004; Bouwens et al. 2009; Vanzella et al. 2009) to separate LBGs from lower redshift Balmer break and dusty galaxies, which show redder colour at longer wavelengths. However, as already noted by Beckwith et al. (2006), the redshift selection range provided by these criteria has a non-negligible dependence on the amount of extinction. For example, Beckwith et al. (2006) (Tables 8 and 9) show that a 100 Myr star-forming model with E(B − V) = 0.15 can be included in the selection from lower redshifts (z_min = 3.3) than can a dust-free one (z_min = 3.5).

3.1. Colour selection

In this paper, thanks to the very deep WFC3 near infrared images, we have adopted a modified version of the B-dropout selection that is based on the (V − H) colour instead of the (V − Z) one. We used models taken from Charlot and Bruzual 2007 (Bruzual 2007a,b, hereafter CB07) to empirically define colour criteria allowing for a clean selection of z ~ 4 LBGs: $\begin{array}{ccc}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{)}& \mathit{>}& \mathrm{1.1}\\ \mathrm{(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{)}& \mathit{>}& \mathrm{1.2}\mathrm{+}\mathrm{0.8}\mathrm{(}\mathit{V}\mathrm{-}\mathit{H}\mathrm{)}\mathit{.}\end{array}$In Fig. 1 we compare the commonly adopted selection window (left panel) to the one used here (right panel) for the CB07 models, showing the position of LBG models with Z/Z_⊙ = 0.2 − 1.0 and Age = 100–300 Myr from z = 3 to z = 5 at 0.5 intervals, at an increasing extinction from E(B − V) = 0.0 to E(B − V) = 1.0 at 0.1 intervals. We also show the position of galaxies at z = 0.0–3.5 with short star formation exponential timescales (0.1–1 Gyr) and ages > 1 Gyr, and constant star-forming models with 0.5 < E(B − V) < 1.5 in the same redshift range. It appears evident that the (V − Z) colour cannot efficiently distinguish between B-dropout and lower redshift contaminants, since z ≳ 3.5 LBGs with moderate E(B − V) and both lower redshift dusty star-forming galaxies and passive ones occupy the same region of the colour diagram. This effect potentially leads to a contamination in the sample of the more extincted B-dropout galaxies, biasing the estimate of the UV slope distribution at z ~ 4, which cannot be completely avoided by redefining the relevant colour criteria. The (V − H)-based criterion instead allows us to remove the reddest potential interlopers from the sample.

Moreover, the (V − H) colour provides a cleaner LBG selection window. A careful look at the reference LBG tracks in Fig. 1 shows that the (B − V) − (V − H) criteria are, in principle, unbiased since they are met by LBGs at z > 3.5 regardless of their E(B − V), age and metallicity. Specifically, there is no need to a priori adopt a cut in the (V − H) colour, unlike the (V − Z) one, and any intrinsic constraint on the (V − H) colour due to observational limits corresponds to the same cut in E(B − V) for galaxies of different redshifts. On the other hand, the usual (B − V) − (V − Z) cut is redshift-dependent, with galaxies at increasing E(B − V) meeting the selection criteria at decreasing redshifts, in agreement with the findings of Beckwith et al. (2006).

3.2. The ERS+HUDF B-dropout sample

Our LBG selection is naturally limited by the capability of accurately measuring large (V − H) colour terms, that is by the relative depth of the optical and IR data-set. We limited our sample to S/N(H) >  10σ to consider only objects with accurately measured colours. This choice implies that the V-band magnitude of an object with E(B − V) = 0.5 (corresponding to (V − H) ~ 2, i.e. β ~ 0) is measured at ~3σ. On the other hand, to exclude z ≳ 4.5 galaxies, we did not consider objects with S/N < 3 in the V band: in practice, given the conservative S/N cut in the H band, the latter criterion only removes one bright V-dropout galaxy from the ERS selection. Finally, to ensure that no AGN are present among the selected objects, we removed three point-like sources (Sextractor CLASS_STAR=0.99) from the ERS sample and check on the 4Ms CDFS images that none of our candidates present X-ray emission (Fiore et al. 2012). The final sample consists of 142 objects at 23.9 < H < 26.5 from the ERS field and 25 objects (24.8 < H < 27.4) from the HUDF.

Fig. 2 Upper panels: colour–colour diagrams for B-dropouts in the ERS (left) and HUDF (right) fields. Selected objects are marked in red, and the blue points indicate B-dropout candidates in the (B − V) vs. (V − Z) selection window excluded by the present analysis. Lower panel: combined photometric redshift probability distribution of our ERS B-dropout candidates (black line) and of galaxies selected from the (B − V) vs. (V − Z) diagram (blue dashed line): the (V − H) selected objects have a lower probability of having zphot ≲ 3.5 solutions.

We show in Fig. 2 the relevant colour-colour diagrams and the differences between the criteria adopted in the present paper and the usual B-dropout selection. Our selected candidates are indicated by red points, while the blue symbols mark the position in the diagram of objects that would be included by a standard (V − Z) colour-selection, but are excluded in our sample. Following the previous analysis of the performance of the two different selections on CB07 models, the position of the latter in the diagrams is consistent with their being either low redshift interlopers or U-dropouts at z ~ 3.0−3.5.

This can be verified through a simple test by exploiting our ERS photometric redshift catalogue that also exploits deep K band and IRAC observations (Santini et al. 2012). We note that the photometric redshifts in the ERS field have been obtained through the well-tested χ² minimization procedure outlined in Grazian et al. (2006) and Santini et al. (2009), which is based on a different set of synthetic spectral templates (PEGASE 2.0 by Fioc & Rocca-Volmerange 1997) than those (CB07) used in Fig. 1. We determined the combined photometric redshift distributions P(z) for the two selection criteria as the sum of the redshift probability distribution of the relevant B-dropout candidates. As shown in the lower panel of Fig. 2, when objects selected through the (V − Z) colour are considered, P(z) displays a lower redshift tail extending out to z ~ 2.5. On the other hand, the (V − H) selected sample is characterized by a sharper P(z) at high redshift and has a lower probability for secondary solutions at z < 1: by integrating P(z) we estimate that the contamination in a sample selected through the (V − Z) colour is more than twice what it is in a (V − H) selected one, with 32% objects at z < 3.5 (14% at z < 3.0), compared to 15% (5% at z < 3.0) in the latter case. Photometric redshifts solutions at z < 1 are due to the similarity between the optical and near-IR SEDs of Balmer break galaxies at these redshifts, and those of B-dropout galaxies. In this respect, an improved colour-based selection is not feasible, since it would require deeper Spitzer data to enable a thorough study of the LBG population.

We also exploited ERS photometric redshifts to explore different criteria for both the (B − V) − (V − H) diagram and the (B − V) − (V − Z) one. In particular, we considered 1) 172 candidates selected through a less conservative (V − H) criterion, i.e. having the redder selection boundary (B − V) > 1.1 + 0.7(V − H), and 2) 129 B-dropouts isolated through a “steeper” (V − Z) colour cut aimed at avoiding a redshift-dependent selection of galaxies with different E(B − V)¹. In the first case we find an increase in the contamination from z ~ 3−3.5 LBGs (P(z < 3.5) = 18%), indicating that this new boundary is indeed extending towards lower values the redshift selection window, consistently with the LBG tracks displayed in Fig. 1. In the second case, we still find higher contamination than in our original (B − V) − (V − H) sample, with a 26% (12%) probability of objects at z < 3.5 (z < 3.0). These findings agree with the tests on CB07 models (Sect. 3.1) and show that our proposed criteria provide a more reliable selection for B-dropout galaxies.

Objects having spectroscopically confirmed redshift further validates the efficiency of our method. Out of ten objects with 3.5 ≤ z_spec ≤ 4.5 and S/N(H) > 10, the (V − H) selected sample includes nine with 3.58 ≤ z_spec ≤ 4.4 and only fails to select an LBG at z_spec = 4.41. Instead, the (V − Z) colour criterion also includes a galaxy at z_spec = 3.11 along with those already present in our sample. This is also consistent with the result of intensive follow-up spectrocopic campaigns of B-dropout samples across all the GOODS-South field: Vanzella et al. (2009) finds that eight out of 46 B-dropout galaxies initially selected through the (B − V) − (V − Z) diagram lie in the range 3.3 ≤ z_spec ≤ 3.5. A similar result is evident from Fig. 7 of Stark et al. (2010) showing the redshift distribution of their spectroscopic sample over GOODS-North. We show in the next section and in Fig. 6 how including these potential interlopers in our sample would affect the determination of the UV slope distribution.

4. Estimating UV slopes and rest-frame magnitudes

We adopt the common power-law approximation for the UV spectral range F_λ ∝ λ^β, and estimate the slope β by fitting a linear relation through the observed I,Z,Y,J AB magnitudes of the objects: ${\mathit{M}}_{\mathit{i}}\mathrm{=}\mathrm{-}\mathrm{2.5}\mathrm{(}\mathit{\beta }\mathrm{+}\mathrm{2.0}\mathrm{)}\log \mathrm{\left(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\right)}\mathrm{+}\mathit{c}$(1)where M_i is the magnitude in the ith filter at effective wavelength λ_i. The filters we use span the rest-frame wavelength range λ ≃ 1500 − 2500 Å at z ~ 4, providing a β estimate that is consistent with the original definition by Calzetti et al. (1994) and that is independent of the exact redshift of the sources, which only affects the intercept c in Eq. (1). We underline that this estimate of β is unbiased with respect to the colour selection being based on a different set of filters. We present in Fig. 3 four examples taken from the ERS sample, along with a comparison between the UV slopes determined by this method and the relevant best-fit SEDs.

The large wavelength baseline we adopt gives β estimates that are significantly more stable than estimates based on the I − Z colour alone (e.g. Bouwens et al. 2009): as an example, we compare in Fig. 4 the β uncertainty distribution of the two methods for H ~ 27 objects from the HUDF mock catalogues presented in Sect. 5.1.

Fig. 3 A comparison of the fitted UV slope (red continuous line) and best-fit SED of four objects in our sample having best-fit photometric redshift in the range zphot = 3.5 − 4.1. The observed magnitudes (filled squares and errorbars) are, from left to right: B,V,I,Z,Y,J,H. The linear-fit is performed on I,Z,Y,J. The lower and upper horizontal axes show the observer-frame and the rest-frame wavelengths, respectively.

Fig. 4 Comparison of uncertainties on β derived from a linear fit to the observed I,Z,Y,J magnitudes (histogram) and from the (I − Z) colour only (red curve) for simulated objects with observed UV slope βfit. The real UV slope βinput of each model is measured on its unperturbed magnitudes.

Fig. 5 The UV slope β of the ERS (open triangles) and HUDF (filled squares) B-dropout objects as a function of their observed H band magnitude. Error bars are given by the uncertainty on the linear fit on their observed I,Z,Y,J AB magnitudes.

Figure 5 shows the resulting UV slope β for the objects in our sample as a function of their H band magnitude: most of the galaxies have blue UV slopes (β ~ −2), apart from a few red LBGs that are mainly populating the brightest end of the diagram. We explored how the distribution of β would change when adding objects to our sample that are included by a standard B-dropout selection diagram, but that are not selected when their (V-H) colour is considered. We performed this test on the ERS to also explore their photometric redshift distribution. As shown in Fig. 6, these candidates are on average redder than those already present in our sample, with 38 out of a total of 47 having a best-fit photometric redshift z_phot < 3.5. Moreover, the only spectroscopically confirmed object at z_spec = 3.11 that we exclude through the (V − H) criterion has a very red UV slope with β = −0.75. Notably, all the reddest objects (β ~ 0) of the (V − Z) selected sample can be considered interlopers on the basis of their photo-z. These results are consistent with the analysis of LBG models presented in Sect. 3.1, showing that the usual (B − V) − (V − Z) cut mainly selects objects with large E(B − V) at z ~ 3−3.5.

A more detailed analysis of the LBG population can be given by looking at how UV slopes change as a function of the rest-frame UV magnitude: to this aim, a linear fit of the UV SED gives us a further advantage, namely we can interpolate the slope of each object to estimate the rest-frame magnitude at 1600 Å. To provide an estimate for M₁₆₀₀, we consider all objects as located at the same redshift, which for simplicity we assume to be z = 4: the effects of this assumption are self-consistenly taken into account in the statistical analysis discussed in the following sections.

Fig. 6 Distribution of the UV slope values from the ERS sample (black hollow histogram). The magenta dashed histogram represents the objects that are not included in our sample but are selected through a standard (V − Z) colour cut. Most of them have zphot < 3.5 (blue filled histogram). The blue arrow marks the position of a zspec = 3.11 object included in this last sample.

5. A statistical analysis of the UV slope distribution

5.1. Uncertainties and systematics in β and M₁₆₀₀ estimates

We first start with a set of Monte Carlo simulations that reproduce the overall properties of our observations. We use the CB07 library to produce a large set of simulated galaxies at 3.0 < z < 5.0 with expected magnitudes in the same filter set as in the ERS and HUDF observations. We consider models with constant star formation histories and within the following range of physical parameters: Z/Z_⊙ = 0.02,0.2,1.0, 0.01 ≤ E(B − V) ≤ 1.0, 0.01 ≤ Age ≤ 1 Gyr. We assume a Salpeter (1955) IMF and a Calzetti et al. (2000) attenuation law, while the transmission of the IGM is treated according to Fan et al. (2006). We exploit this library to generate mock catalogues of 5 × 10⁵ objects for each field by perturbing magnitudes in all the bands to match the relevant depths in the observed ones. Our simulations are designed to reproduce both the average and the scatter of the S/N distribution as a function of magnitude to provide the closest match between observed and simulated data. The simulated catalogues are then analysed like the real ones to determine the observational effects on galaxies of different magnitudes and spectral slopes.

While being much less time-consuming than two-dimensional imaging simulations, our Monte Carlo approach allows us to estimate the effect of magnitude and colour uncertainties in the dropout selection and UV slope analysis. A further advantage of such a purely empirical approach is that it is free of assumptions on the typical size and size-luminosity relation of the LBGs, which, as shown in Grazian et al. (2011), are a critical factor in estimating the completeness correction for faint objects. We randomly select a subsample of 20 000 model LBGs from the same CB07 sample and perform a simple two-dimensional simulation to compare its results to those obtained through our empirical method. We adopt as morphology template the stacking of ten bright candidates from the ERS field, and perform various runs inserting 200 objects in the images each time after masking the regions where real objects were detected. The inserted objects are then recovered and analysed with the same procedure as used for the real catalogues. A comparison between simulated objects obtained in this way and the mock catalogues discussed above shows good consistency, therefore we base our analysis on the latter ones, which allow us to perform larger and less time-consuming simulations.

Fig. 7 From top to bottom: detection completeness and difference between output and input β and M1600 from our ERS (left) and HUDF (right) mock catalogues. All quantities are computed for objects in three observed magnitude bins centred at: H = 24.5 (ERS) and H = 25.5 (HUDF) (black continuous lines and filled squares); H = 25.25 (ERS) and H = 26.3 (HUDF) (dashed blue lines and open squares); H = 26.0 (ERS) and H = 27.0 (HUDF) (red dot-dashed lines and crosses).

In the uppermost panel of Fig. 7, we show for both fields the detection completeness at different observed fluxes of galaxies with a given spectral slope, as indicated by the β_input of the model measured on its real unperturbed magnitudes. For the ERS sample the completeness curves remain approximately flat for β_input ≲ −1.5 and rapidly decrease at β_input ≳ −1, especially for faint objects. The completeness in the HUDF is higher because of the deeper optical images: it remains close to one up to β_input ≲ −0.5 for bright objects, while fainter ones are still efficiently selected as long as β_input ≲ −1.0. The combination of the two fields thus allows us to sample a wide range of UV slopes and magnitudes.

In Fig. 7, we also compare input and measured values of M₁₆₀₀ and β for galaxies taken from our simulations. In particular, we explore the dependence of uncertainties and systematics on the spectral shape of the galaxy, as indicated by β_input, and on the observed magnitude. The latter test is of particular importance since we are analysing the colour distribution of a flux-limited sample. The uncertainties on both M₁₆₀₀ and β are relatively small, with σ_M1600 ≃ 0.1−0.2 and σ_β ≃ 0.1−0.3. However, they become larger for faint and red objects in the sample with respect to blue/luminous ones. The measured UV slope β_fit does not show any large systematic offset as long as β_input < −1, a small trend towards underestimated values appears for the faintest red objects in the ERS sample. This is due to the shallower optical data in that field. The estimated UV rest-frame magnitude M₁₆₀₀ shows a negligible offset of Δ ~ 0.05−0.1 mags with respect to the real value for β_input < −1 with redder objects being underestimated by Δ ~ 0.15−0.2 mags. The conservative selection S/N(H) > 10 helps keep the uncertainties small. A similar test for objects of S/N(H) ~ 5−10 shows that the uncertainty on the UV slope is as large as 0.6 for faint galaxies.

Fig. 8 Maximum likelihood contours (68%, 95%, and 99% c.l. in green, blue, and red, respectively) for the Gaussian and lognormal fit at different magnitudes. The histograms are the observed data. The red continuous line in each bin is the relevant best-fit distribution. The blue line show how it appears after being “convolved” by observational effects.

5.2. Reconstructing the UV slope distribution: the method

The tests discussed in Sect. 5.1 show that the completeness and the accuracy on both β and UV rest-frame magnitude estimates depend on the observed magnitude, as well as on the colour of the galaxies. This implies that any measured UV slope distributions at different luminosites arise from the underlying intrinsic one after being convolved by a combination of observational effects. To estimate the dependence of β on UV magnitude, it is therefore not sufficient to apply simple corrections to the observed distributions, but it is necessary to adopt a more comprehensive approach. Such an approach has to be based on assumptions on the properties of the intrinsic UV slope distribution, as for any “de-convolution” procedure. We apply here a maximum likelihood technique which is analogous to the one used to estimate the galaxy luminosity function (LF) where a Schechter LF is assumed, and its best-fit parameters are obtained by comparing simulated number counts to the observed ones. We start from the assumption that the probability distribution function of UV slopes in a given magnitude range, PDF(β), follows a given functional form, whose average and standard deviation are to be determined. For any given set of parameters and for each magnitude bin, we extract from our CB07 libraries 50 000 objects reproducing the given distribution and perturb them with noise as discussed in Sect. 5.1, to introduce the effect of observational uncertainties and biases. Finally, we compare the output distribution of β values to the observed one in each bin through a maximum likelihood estimator ℒ (e.g. Bouwens et al. 2008; Castellano et al. 2010) and minimize Δχ² = −2.0ln(ℒ). Unfortunately, while a Schechter form in the case of the luminosity function is suggested by theoretical arguments and has been verified through many samples at low and high redshift, the present work constitutes the first attempt of a parametrical reconstruction of the UV slope distribution PDF(β). We are thus forced to start from a reasonable hypothesis on the analytic shape of PDF(β), as suggested by our data and by previous non-parametric analysis. We adopt two different analytic PDF(β) as a working hypothesis. We first assume a Gaussian distribution for the β values. We note that such a distribution has already been suggested by other authors in the past (e.g. Bouwens et al. 2009). As an alternative, we assume a log-normal distribution PDF(β) with β = −2.5 as the lowest, asymptotic value. This is motivated by the asymmetric distribution observed in the brightest magnitude bins (Fig. 8) and corresponds to the hypothesis that the high-z galaxy population is made of a bulk of low-extincted star-forming LBGs along with a number of rare objects at increasing E(B − V) values.

We change the average and dispersion of both input distributions across a 25 × 25 grid in the range ⟨ β ⟩ = −2.5, −1.5 and rms (β) = 0.03,0.7.

Fig. 9 The β − M1600 diagram for our B-dropout sample. Open triangles and filled squares represent individual objects in the ERS and HUDF fields, respectively. Error bars on the observed points are given by the uncertainty on the linear fit on their observed magnitudes. Blue and red circles and lines indicate the result of the statistical analysis described in Sect. 5, which takes systematic effects into account. Blue open circles and error bars indicate the average ⟨ β ⟩  and intrinsic r.m.s of the best-fit Gaussian distributions in each of the four bins: uncertainties on ⟨ β ⟩  are much lower and are not visible in this scale. Red filled circles show the best-fit log-normal distributions. The Gaussians best-fit values are plotted at the average magnitude of the observed points in the relevant bins, log-normal ones with a small magnitude offset for clarity. The green continuous line is the β − UV relation by Bouwens et al. (2009). We show for comparison a linear fit on the Gaussian and log-normal best-fit points (blue dashed and red continous lines), but underline that our analysis indicates a significant difference only between the first and second magnitude bins.

5.3. Results

To analyse the relation between magnitude and UV slope, we divide the sample into four M₁₆₀₀ bins (Fig. 8). This choice is somewhat arbitrary and is motivated by the competing needs of having magnitude bins that are as narrow as possible, but with enough objects to allow the statistical analysis described in the following. In practice, given the different magnitude ranges probed by the two fields, the ERS objects will be used to constrain the first and second bins, the third one is analysed by combining data from both fields, while HUDF objects is used for the faintest magnitude bin.

The resulting best-fit distributions (Table 1) indicate that most of the B-dropouts have blue UV slopes: in particular, at L < L^∗ galaxies tend to be very blue regardless of their M₁₆₀₀, while the brightest objects only (first bin) are, on average, slightly redder than faint ones. Indeed, the best-fit Gaussians are remarkably similar in the three magnitude bins at M₁₆₀₀ > − 21.0 with an average of ⟨ β ⟩ ≃ −2.1 and an intrinsic rms ≃ 0.3, while brighter galaxies have a best-fit Gaussian with ⟨ β ⟩ ≃ − 1.9, but similar dispersion. Given the low uncertainties, this small difference between bright and faint galaxies is significant at ~2σ. The log-normal fit yields similar results, with bright galaxies (first bin) redder than M₁₆₀₀ > −21.0 objects, the only difference being that in this case, given the asymmetric shape of the distribution, the bulk of the objects would have β ~ −2.0 in the first bin and β ~ −2.2 at M₁₆₀₀ ≳ − 21.0. We verified that these results are robust against small changes in the colour selection criteria: objects selected with the redder boundary in the (B − V) vs. (V − H) diagram discussed in Sect. 3.2 yields best-fit values within the 68% c.l. region given by our “fiducial” sample.

The binned data and best-fit distributions along with the maximum likelihood contours for the parameters are shown in Fig. 8. The need for an approach like the one carried out here can be appreciated by looking at how the input best-fit distributions in the various magnitude bins (red curves in Fig. 8) are modified when they are self-consistently “convolved” with observational effects through simulations. Finally, Fig. 9 shows the distribution of observed points in the β − M₁₆₀₀ diagram along with the average and intrinsic dispersion of the best-fit Gaussian and lognormal distribution.

While our tests consistently indicate that most of the B-dropouts have blue UV slopes, it does not allow us to determine which is the functional form that best reproduces the data. To this aim, we compared the observed UV slope distributions to the best-fit Gaussian and log-normal PDF(β) in each bin through a Kolmogorov-Smirnov test. We find a significant difference only in the first bin where the Gaussian fit yields a low KS probability ( < 50%) because of the presence of a tail with few red objects at β ≃ − 1.0. In the remaining bins the KS test does not allow us to distinguish between the two distributions. It is interesting to note that a constant distribution of β values within each bin is completely excluded by our maximum likelihood test at all magnitudes. However, we cannot exclude that the distribution is bimodal, with a “red” peak at values of β larger than accessible to Lyman Break selections, or that it is a more complicated one, as it would be the case of Gaussians with non-standard kurtosis and/or skewness. While the first case can only be probed by combining LBG-selected samples with far-IR and submillimetre selected ones (e.g. Wardlow et al. 2011), the latter would require fitting more parameters than is meaningful to do with the presently available data.

The physical interpretation of these results is not straightforward, since several factors determine the UV slope: dust reddening, varying ages, metallicity, IMF, and star formation histories (e.g. Bouwens et al. 2009; Wilkins et al. 2011). However, the effect of dust is by far the most relevant, and with the presently available data it is not possible to constrain how much other properties contribute to the variation in the UV slope distribution with magnitude, or if they only affect the intrinsic dispersion of the PDF(β) distributions. For this reason, and to provide comparison with other works, we interpret our results in the following by directly converting the UV slope into dust extinction.

In this context, our results imply that the average extinction at 1600 Å (following Meurer et al. 1999) is A₁₆₀₀ ≃ 0.65 in the first bin and A₁₆₀₀ ≃ 0.25 in the faintest bins. Deviations from the Meurer et al. (1999) relation have been found in z ≳ 4 QSOs by Gallerani et al. (2010). Their MEC attenuation curve would imply a higher attenuation, A₁₆₀₀ ≃ 1.4 for bright galaxies and A₁₆₀₀ ≃ 0.5 for faint ones.

5.4. Comparison with previous results

Previous determinations of the UV extinction of z ~ 4 LBGs have been provided by Ouchi et al. (2004), Beckwith et al. (2006), Overzier et al. (2008)Bouwens et al. (2009), and Lee et al. (2011). Ouchi et al. (2004) found an average E(B − V) = 0.15 ± 0.03 for their M < M^∗ B-dropouts and no dependence on the apparent (i.e. uncorrected) luminosity within the magnitude range they probe. Our estimate at similar magnitudes (M₁₆₀₀ < −21.0) is lower: when converting the relevant best-fit Gaussian PDF(β) into a corresponding E(B − V) distribution following the equations in Calzetti et al. (2000), we obtain an average $\mathrm{⟨}\mathit{E}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{)}\mathrm{⟩}\mathrm{=}\mathrm{0.0}{\mathrm{5}}_{-0.02}^{\mathrm{+}\mathrm{0.03}}$. On the other hand, Lee et al. (2011) find a strong dependence of β on UV magnitude in a large sample of M < M^∗z ~ 3.7 galaxies from the NOAO Deep Wide Field Survey. Although they probe much brighter magnitudes than those available in our sample, a comparison can be made between their faintest bin centred at M₁₇₀₀ = −21.43 and the brightest galaxies from our ERS sample. They find β = −1.78 ± 0.28 at these magnitudes, which falls within the 68% c.l. region given by our maximum likelihood test.

Beckwith et al. (2006) find blue slopes for both bright and faint z ~ 4 galaxies in the HUDF, suggesting that they are affected by little dust extinction. A quantitative comparison is not possible because Beckwith et al. (2006) do not explicitly compute the value of the UV slope in the different samples. However, their result agrees with our findings of blue UV slopes in all L ≲ L_∗ magnitude bins, which is the range probed by the HUDF.

Overzier et al. (2008) provides an estimate of the UV slope of LBGs and Lyman α emitters (LAE) in a z ~ 4.1 protocluster finding a global average ⟨ β ⟩ = −1.95 which agrees with our results. On the other hand, our results are markedly different from those obtained by Bouwens et al. (2009). They find redder β averages and a steep correlation between β and UV magnitude up to very faint luminosities: their relation is shown in Fig. 9, where for comparison we also show a simple linear regression of our best-fit Gaussian (log-normal) average slope ⟨ β ⟩  as a function of M₁₆₀₀.

Fig. 10 The relation between β and the dust-corrected UV magnitude for ERS (open triangles) and HUDF objects (filled squares) in our sample. The upper and the right axes show the conversion of unextincted magnitude and UV slope in SFR (Madau et al. 1998) and A1600 (Meurer et al. 1999), respectively. The red continuous lines and the blue dashed ones indicate the 90% and 30% selection completeness in the two fields.

These discrepancies can be explained in part by the differences between various approaches, especially by the different colour selections as discussed in Sect. 3. Specifically, we explored whether the greatest differences between our analysis and the Bouwens et al. (2009) one can be responsible of the discrepancies: we mimicked their approach by adopting the B − V, V − Z criterion (i.e. including in our sample the objects shown in Fig. 6) and we adopted the same definitions for β and M_UV based on the I and Z observed magnitudes. With the inclusion of the additional red candidates, at bright magnitudes we obtain average UV slopes redder (β ~ −1.7) than in the case when potential interlopers are excluded through their (V − H) colour. However, even with this inclusion, the values still are not as red as those found by Bouwens et al. (2009)¹.

A completely different approach in the estimate of dust extinction in high-redshift galaxies is provided by stacking analysis at radio wavelengths to determine their average uncorrected SFR. In this respect, the analysis of z ~ 4 LBGs by Ho et al. (2010), based on VLA 1.4GHz observations, favours lower values for dust extinction than those inferred from past UV analysis, in agreement with our findings. However, a comparison between the two approaches is complicated by several factors and assumptions (see discussion in Carilli et al. 2008).

6. The dust corrected star formation rate density

The parametric estimate of the UV slope distribution at different magnitudes gives us the remarkable opportunity of computing the SFRD in a straightforward way through the equation: $\mathrm{SFRD}\mathrm{=}\frac{\mathrm{1.0}}{\mathrm{8}\mathrm{\times }{\mathrm{10}}^{\mathrm{27}}}\mathrm{\int }\mathrm{d}\mathit{L}\mathrm{\int }\mathrm{d}\mathit{A}\mathrm{·}\mathit{PDF}\mathrm{\left(}\mathit{A,L}\mathrm{\right)}\mathrm{\times }{\mathrm{10}}^{\mathrm{0.4}\mathrm{·}\mathit{A}}\mathit{L}\mathrm{·}\mathrm{\Phi }\mathrm{\left(}\mathit{L}\mathrm{\right)}\mathit{,}$(2)where the constant factor is from Madau et al. (1998), A = A₁₆₀₀ such that PDF(A,L) is univocally related to the PDF(β,L) through the Meurer et al. (1999) equation, and Φ(L) is the UV luminosity function at 1600 Å. For the last we adopt the best-fit Schechter parameters by Bouwens et al. (2007) for consistency with previous results. In practice, we compute the integral in each luminosity bin by considering the relevant best-fit distributions PDF(β) and extrapolating the last one up to M₁₆₀₀ = −17.5 (corresponding to L = 0.04  L ^∗ ). When the Gaussian PDF(β) are considered we obtain log(SFRD) = $\mathrm{-}\mathrm{1.09}{\mathrm{3}}_{-0.18}^{\mathrm{+}\mathrm{0.18}} {\mathit{M}}_{\mathrm{\odot }}$/yr/Mpc³, taking both the 68% c.l. uncertainty on each PDF(β) and the uncertainty on the LF parameters into account. Using the best-fit lognormal PDF(β) in Eq. (2) we obtain a slightly higher value, log(SFRD) = $\mathrm{-}\mathrm{1.08}{\mathrm{6}}_{-0.18}^{\mathrm{+}\mathrm{0.20}} {\mathit{M}}_{\mathrm{\odot }}$/yr/Mpc³, because of the larger effective correction for dust extincion. Our SFRD estimates are ~55% lower than the one by Bouwens et al. (2009): this is not a surprise since we find bluer β at bright and intermediate magnitudes, so we apply a lower correction for dust. Interestingly, our SFRD values are consistent with the estimate obtained by Thompson et al. (2006) on a photometric-redshift selected z ~ 4 sample and through the use of SED-fitting to estimate extinction correction. This agrees with the discussion in Sects. 3 and 4 on the differences between our H-based selection and the usual optical-based one.

The largest uncertainties in the SFRD estimate are given by our lack of knowledge of the attenuation curve at these redshifts and by the possible presence of highly-obscured star-forming objects, like the submillimetre galaxies (SMGs) studied at lower redshifts (e.g. Chapman et al. 2005). As an example, the MEC attenuation curve by Gallerani et al. (2010), which is flatter than the Calzetti et al. one in the UV, yields an SFRD that is nearly three times larger: log(SFRD) ≃ − 0.65  M_⊙/yr/Mpc³. On the other hand, a significant contribution from z ~ 4 SMGs (e.g. Daddi et al. 2009) would imply higher values of the SFRD than those measured from LBGs alone.

7. Dust extinction as a function of intrinsic UV luminosity

Galaxy evolutionary models suggest that UV bright galaxies tend to be dustier than faint ones, as a natural outcome of the relation between SFR, mass, and metallicity. Our findings of a redder UV slope distribution at bright magnitudes apparently confirm this prediction. However, such analysis shows how galaxies of different UV slopes (which is a proxy of their dust content) are distributed as a function of their extincted UV magnitude. The effect of extinction is thus present on both quantities, so a clear characterization of how the dust content varies across the galaxy population is not possible. The usual β-UV magnitude relation, as presented in Sect. 5, is mainly useful for estimating a correction to the observed SFRD, as discussed above, but do not offer insight into the physical properties of Lyman break galaxies.

To provide a clearer view on this issue we show the relation between β and the dust-corrected UV magnitude M_1600corr = M₁₆₀₀ − A₁₆₀₀ in Fig. 10. We also show the conversion of both in more grounded physical quantities: A₁₆₀₀ (following Meurer et al. 1999) and SFR (Madau et al. 1998). Naturally, degeneracies between extinction and other physical parameters can bias this direct conversion, and quantities are not completely independent, since the extinction correction is itself related to the β of the single objects. Nonetheless, the simplicity of our analysis provides a more straightforward view of such observational degeneracies than in the case of common SED-fitting estimates. We also plot in Fig. 10 the 30% and 90% detection completeness levels, as estimated from our simulations, for the galaxies of the two fields.

From this plot it appears evident that at high intrinsic luminosities (M_1600corr ≲ −23.0, corresponding to SFR ≳ 80  M_⊙/yr), only dusty star-forming galaxies are found. Once the incompleteness of our selection is taken into account we obtain a surface density of 0.12 arcmin^-2 for objects with SFR = > 80  M_⊙/yr when 2 < A₁₆₀₀ < 3 and an upper limit of 0.02 arcmin^-2 for less extincted galaxies.

Among galaxies with lower SFR, we detect many with a much lower amount of reddening. As a result a trend toward decreasing extinction as a function of corrected UV luminosity seems to be in place. However, this effect may be entirely a systematic effect of the limits of our observations (as shown by the lines in Fig. 10) that do not allow us to ascertain whether this trend is real or it is due to selection effects even at moderately high luminosities.

8. Summary and conclusions

In this paper we have presented a new approach to the selection and to the analysis of the UV slope distribution of z ~ 4 LBGs. Thanks to the deep WFC3 near infrared observations of the ERS and HUDF fields, combined with the existing HST-ACS optical data, we performed a B-dropout selection that is based on a (B − V) vs. (V − H) diagram on H-detected samples of both fields. Following the analysis of galaxy spectral libraries and of our own photometric redshift catalogue on the ERS, we showed that this selection allows us to efficiently exclude low-redshift interlopers and dusty z ≳ 3 galaxies that might be present in B-dropout samples selected through optical colours. The sample we selected on the basis of the newly defined colour criteria comprises 142 and 25 objects at S/N(H) >  10 in the ERS and HUDF fields, respectively. We estimated the UV slope of our B-dropout candidates through a linear fit on their I,Z,Y,J magnitudes, and we interpolated the same fit to provide an estimate of their UV rest-frame magnitude after assuming that all sources are located at z = 4.

Despite the conservative selection and the greater accuracy in estimating β provided by a linear fit, we showed that uncertainties and systematics arise when measuring the colour distribution of such flux-limited samples. However, the effect of observational uncertainties can be effectively taken into account by adopting a parametric approach combined with detailed simulations, and the estimate of the typical UV slope at different magnitudes can be performed following hypothesis on the shape of the undelying intrinsic distribution. We carried out a parametric maximum likelihood analysis by assuming that the distribution of β is either Gaussian or lognormal. We allowed the relevant parameters of the two distributions to vary across a wide range of values and found the best-fit values by comparing the simulated and observed counts through a maximum-likelihood estimator.

We found that z ~ 4 LBGs are characterized by blue UV slopes, suggesting a low dust extinction. The best-fit Gaussians are remarkably similar in all the magnitude bins at M₁₆₀₀ > −21.0 with an average of ⟨ β ⟩ ≃ −2.1 and an intrinsic rms ≃ 0.3, while brighter galaxies have a best-fit Gaussian with slightly redder average (⟨β⟩ ≃ −1.9) but similar dispersion. A similar result was obtained when log-normal distributions are considered, the main difference being that in this case, given the asymmetric shape of the distribution, the bulk of the objects would have β ~ −2.0 in the first bin and β ~ −2.2 at M₁₆₀₀ ≳ −21.0. When interepreted in terms of varying dust extinction, our results imply that z ~ 4 LBGs have a low dust content, with a typical A₁₆₀₀ ≃ 0.65 at L < L^∗ and A₁₆₀₀ ≃ 0.25 for faintest galaxies, following Meurer et al. (1999).

Our results generally agrees with the previous findings of Beckwith et al. (2006) and Overzier et al. (2008), and are consistent within the uncertainties with the result on L ~ L^∗ galaxies by Lee et al. (2011). On the other hand, our results differ from those obtained by Ouchi et al. (2004) for bright galaxies and by Bouwens et al. (2009), that find redder β averages and a steeper correlation between β and UV magnitude. We find that these discrepancies are only in part due to the differences in the selection (IR-based vs. optically-based) and in the analysis of the B-dropout samples, since such red averages for the UV slope are not found in our sample even when potential interlopers selected through their optical colours are included.

The low dust extinction we find points to an SFRD in the LBG population that is significantly revised downward with respect to previous estimates. If we compute the corrected SFRD from the UV LF by exploiting our parametric estimate (Eq. (2)), we obtain log(SFRD) ~ −1.09, which is ~ 55% lower than the Bouwens et al. (2009) estimate. However, this result is based on the assumption of a Calzetti et al. attenuation law: when the MEC attenuation curve of Gallerani et al. (2010) is assumed, we obtain a nearly three times higher SFRD at z ~ 4.

Finally, we discuss how the UV slope of z ~ 4 galaxies changes as a function of the dust-corrected UV magnitude, i.e. as a function of the SFR. Because of the difficulties in selecting highly extincted objects, we show that constraints can only be given on the number density of UV bright galaxies (SFR ≳ 80  M_⊙/yr). We find that highly extincted galaxies with 2 < A₁₆₀₀ < 3 are much more common (surface density of 0.12 arcmin^-2) than dust-free ones (<0.02 arcmin^-2). This implies that most galaxies with a high SFR (≳80  M_⊙/yr) are highly extincted objects. Among galaxies with lower SFR, we detect many with a much lower amount of reddening, although current observational limits prevent us from detecting those with high exctinction, if they exist.

Although the now available WFC3 observations allow for a refined selection of z ~ 4 samples, improvements are needed to further constrain the dust properties at high redshifts. Deeper spectroscopic follow-up observations are needed to validate both IR and pure optical colour selections. Such samples will also give the opportunity of complementing the study of UV slope distributions with the analysis of the Lyα EW distribution, which might also be related to the dust content of galaxies (e.g. Pentericci et al. 2007). Larger near IR data samples, such as those that will be available in the near future thanks to the CANDELS survey, will allow us to deepen this kind of analysis. While CANDELS/Wide fields will enable a comparison of the UV slope distribution of IR-selected LBGs to the past optical-based analysis of large L ≫ L ^∗  samples, the CANDELS/Deep observations will allow us to strengthen statistical constraints at faint/intermediate magnitudes.

Instead, a thorough characterization of the high-redshift population and a comprehensive analysis of the relation between SFR and dust content will require future instruments both in the optical and in the IR capable of observing highly extincted objects with SFR of few solar masses per year.

Acknowledgments

We acknowledge the support of ASI-INAF in the framework of programme I/009/10/0 supporting Data Analysis in the field of Cosmology.

References

All Tables

All Figures

Fig. 1 The (B − V) vs. (V − H) colour selection window adopted in the present paper (right) along with the commonly adopted one based on the (V − Z) colour (left). Blue filled squares and lines are for a CB07 LBG model with Z/Z⊙ = 0.2 and Age = 100 Myr at increasing redshift and E(B − V) respectively. Models having solar metallicity and Age = 100 Myr and Age = 300 Myr are represented in magenta (crosses and long-dashed lines), and red (empty squares and short-dashed lines), respectively. Red and green points indicate lower redshift passively evolving and dusty star-forming galaxies, respectively. See main body for details.

In the text

Fig. 2 Upper panels: colour–colour diagrams for B-dropouts in the ERS (left) and HUDF (right) fields. Selected objects are marked in red, and the blue points indicate B-dropout candidates in the (B − V) vs. (V − Z) selection window excluded by the present analysis. Lower panel: combined photometric redshift probability distribution of our ERS B-dropout candidates (black line) and of galaxies selected from the (B − V) vs. (V − Z) diagram (blue dashed line): the (V − H) selected objects have a lower probability of having zphot ≲ 3.5 solutions.

In the text

Fig. 3 A comparison of the fitted UV slope (red continuous line) and best-fit SED of four objects in our sample having best-fit photometric redshift in the range zphot = 3.5 − 4.1. The observed magnitudes (filled squares and errorbars) are, from left to right: B,V,I,Z,Y,J,H. The linear-fit is performed on I,Z,Y,J. The lower and upper horizontal axes show the observer-frame and the rest-frame wavelengths, respectively.

In the text

	Fig. 4 Comparison of uncertainties on β derived from a linear fit to the observed I,Z,Y,J magnitudes (histogram) and from the (I − Z) colour only (red curve) for simulated objects with observed UV slope βfit. The real UV slope βinput of each model is measured on its unperturbed magnitudes.
In the text

	Fig. 5 The UV slope β of the ERS (open triangles) and HUDF (filled squares) B-dropout objects as a function of their observed H band magnitude. Error bars are given by the uncertainty on the linear fit on their observed I,Z,Y,J AB magnitudes.
In the text

	Fig. 6 Distribution of the UV slope values from the ERS sample (black hollow histogram). The magenta dashed histogram represents the objects that are not included in our sample but are selected through a standard (V − Z) colour cut. Most of them have zphot < 3.5 (blue filled histogram). The blue arrow marks the position of a zspec = 3.11 object included in this last sample.
In the text

Fig. 7 From top to bottom: detection completeness and difference between output and input β and M1600 from our ERS (left) and HUDF (right) mock catalogues. All quantities are computed for objects in three observed magnitude bins centred at: H = 24.5 (ERS) and H = 25.5 (HUDF) (black continuous lines and filled squares); H = 25.25 (ERS) and H = 26.3 (HUDF) (dashed blue lines and open squares); H = 26.0 (ERS) and H = 27.0 (HUDF) (red dot-dashed lines and crosses).

In the text

	Fig. 8 Maximum likelihood contours (68%, 95%, and 99% c.l. in green, blue, and red, respectively) for the Gaussian and lognormal fit at different magnitudes. The histograms are the observed data. The red continuous line in each bin is the relevant best-fit distribution. The blue line show how it appears after being “convolved” by observational effects.
In the text

Fig. 9 The β − M1600 diagram for our B-dropout sample. Open triangles and filled squares represent individual objects in the ERS and HUDF fields, respectively. Error bars on the observed points are given by the uncertainty on the linear fit on their observed magnitudes. Blue and red circles and lines indicate the result of the statistical analysis described in Sect. 5, which takes systematic effects into account. Blue open circles and error bars indicate the average ⟨ β ⟩  and intrinsic r.m.s of the best-fit Gaussian distributions in each of the four bins: uncertainties on ⟨ β ⟩  are much lower and are not visible in this scale. Red filled circles show the best-fit log-normal distributions. The Gaussians best-fit values are plotted at the average magnitude of the observed points in the relevant bins, log-normal ones with a small magnitude offset for clarity. The green continuous line is the β − UV relation by Bouwens et al. (2009). We show for comparison a linear fit on the Gaussian and log-normal best-fit points (blue dashed and red continous lines), but underline that our analysis indicates a significant difference only between the first and second magnitude bins.

In the text

Fig. 10 The relation between β and the dust-corrected UV magnitude for ERS (open triangles) and HUDF objects (filled squares) in our sample. The upper and the right axes show the conversion of unextincted magnitude and UV slope in SFR (Madau et al. 1998) and A1600 (Meurer et al. 1999), respectively. The red continuous lines and the blue dashed ones indicate the 90% and 30% selection completeness in the two fields.

In the text