© ESO, 2011

1. Introduction

The detailed understanding of the physical properties of the first building blocks of present-day galaxies is fundamental to set up the zero point for the theoretical models dealing with the formation and evolution of galactic structures in the Universe. In particular, the study of the faint end of the LBG LF at z = 7 is very important for understanding their contribution to the reionization processes of the high-z Universe.

At z ~ 6 the Universe is almost re-ionized (Fan et al. 2006; Becker et al. 2007) and there is both observational (Meiksin 2009; Songaila & Cowie 2010; Komatsu et al. 2011) and theoretical (Gnedin & Fan 2006; Cen 2010) mounting evidence that the crucial transformation from a neutral to an ionized state should happen around z ~ 9. The exact timescale of the reionization process, however, is not clear, the alternatives going from an extended period (z ~ 7−11, Dunkley et al. 2009) to a more sudden transition (Cen 2003).

The new Subaru (Ouchi et al. 2009), Hawk-I (Castellano et al. 2010a, 2010b, hereafter C10a and C10b), and WFC3 (Bouwens et al. 2011) surveys have started investigating the properties of galaxies during the reionization epoch, enlarging the small sample of galaxy candidates at z ~ 7 sketched by NICMOS (Bouwens & Illingworth 2006; Mannucci et al. 2007; Oesch et al. 2009). Recently, Oesch et al. (2010a) derived the number density of faint LBGs at z = 7 in the HUDF, concluding that these low luminosity galaxies are enough to reionize the Universe, assuming that the IMF, the clumpiness of the IGM, the escape fraction of ionising photons and their metallicity have not significantly evolved with respect to their local properties. Using the same data-set, Bunker et al. (2010) and McLure et al. (2010) found that it would be difficult for the observed population of high-redshift star-forming galaxies to achieve reionization by z ~ 6−7 without a significant contribution from galaxies well below the detection limits, together with significant variations in the escape fraction of ionising photons. Quantifying the number density of z ≥ 7 galaxies is therefore critical in order to check whether there are sources responsible for the re-ionization, such as Pop III dominated primordial galaxies, mini-black holes or others (Venkatesan et al. 2003; Madau et al. 2004).

The different results obtained for the LF of LBGs at z = 7 (C10a, C10b, Ouchi et al. 2009; Capak et al. 2011; Oesch et al. 2010a; McLure et al. 2010; Bunker et al. 2010; Wilkins et al. 2010; Hickey et al. 2010; Yan et al. 2010; Finkelstein et al. 2010) should be ascribed both to the different data-sets used and to the different techniques adopted to analyse the data. We want to study in detail the possible systematics/uncertainties in the LF estimate, in particular due to the choice of the color cut adopted, the morphology of the z = 7 galaxies assumed during simulations, and the statistical tools used in the LF derivation. Moreover, extending the LF at z = 7 down to faint magnitude limits is fundamental to break current degeneracies between M ^∗  and Φ ^∗ , and to put strong constraints on the number density of faint LBGs at high-z. Here we combine a re-analysis of the extremely deep WFC3 observations (Oesch et al. 2010a) with the selection of bona-fide z ~ 7 LBGs in wide areas of the sky (C10a, C10b, Ouchi et al. 2009) in order to derive stringent constraint to the reionization process at z ~ 7.

The paper is organised as follows. Sections 2 and 3 describe the imaging data used, the photometric analysis, and the selection of the z ~ 7 LBG candidates. The simulations carried out to derive incompleteness and contamination from lower-z galaxies are discussed in Sect. 4. In Sect. 5 we discuss the systematics affecting the LF (color criteria and morphology of the galaxies adopted during the simulations), with our results on the LF derivation. Sections 6 and 7 are devoted to our discussions and conclusions.

Throughout the whole paper, 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

We have analysed two different data-sets observed with the new NIR camera of HST, the Wide Field Camera 3 (WFC3¹): the Early Release Science on the GOODS-S field (GOODS-ERS) and the Hubble Ultra Deep Field (HUDF) programs.

2.1. ERS

The GOODS-ERS WFC3/IR dataset is a total of 60 HST orbits consisting of 10 contiguous pointings in the GOODS-South field (HST Program ID 11359), using 3 filters per visit (Y₀₉₈, J₁₂₅, H₁₆₀), and 2 orbits per filter (for a total of 4800−5400 s per pointing and filter). The total area covered by the GOODS-ERS is  ~40   sq.   arcmin down to Y = 27.3, J = 27.4, and H = 27.4 mag at 5σ in an area of  ~0.11 arcsec² (corresponding to 2 times the FWHM of the images). The ERS data were acquired and used to demonstrate the WFC3 capabilities on the study of low-z and high-z galaxies (Windhorst et al. 2011), and the raw images were publicly released by the instrument team. Thanks to its fine pixel sampling and high sensitivity in the NIR, WFC3 is an ideal instrument to study the z ≥ 7 galaxy populations. The raw images of the ERS dataset have been reduced and used to search for zdrop candidates by Wilkins et al. (2010; 2011) and Bouwens et al. (2011).

The WFC3 ERS images were reduced using the Multidrizzle software (Koekemoer et al. 2002), adopting on-orbit SPARS100 darks and also correcting for the gain differences between different quadrants. Large-scale residuals on each exposure which are apparently due to scattered light have been removed, and the satellite trails have been masked on affected exposures. The astrometry of the images is bound to the GOODS-S v2.0 z-band (Giavalisco et al. 2004). The matching is robust even for a single exposure, with 400−600 sources matching, so that the overall alignment seems to be robust to a level of  ~10   mas. The pixel size of the reduced images is 0.06 arcsec (a factor of 2 rebinning was applied with respect to the original pixel size of 0.135 × 0.121 arcsec²), with a PSF of 0.18 arcsec.

2.2. HUDF

The HUDF WFC3/IR dataset (HST Program ID 11563) is a total of 60 HST orbits in a single pointing (Oesch et al. 2010a; Bouwens et al. 2010a) in three broad-band filters (16 orbits in Y₁₀₅, 16 in J₁₂₅, and 28 in H₁₆₀). It is the deepest NIR image ever taken, reaching Y = 29.3, J = 29.2, and H = 29.2 total magnitudes for point like sources at 5 sigma (this S/N ratio is computed in an aperture of  ~0.11 arcsec², corresponding to one FWHM). The area covered by the WFC3-HUDF imaging is 4.7 sq.   arcmin, and the IR data have been drizzled to the ACS-HUDF data (Beckwith et al. 2006), with a resulting pixel scale of 0.03 arcsec, with a FWHM of 0.18 arcsec.

The HUDF images were reduced in the same way as the ERS dataset, i.e. background effects are properly corrected (on-orbit darks, quadrant gain corrections, multiplicative sky-flats and additive low-level residuals on large scales all removed), and the images have been astrometrically matched using the ACS HUDF sources. For the Y₁₀₅ filter, we have selected a subset of 7 visits out of 9 that are not badly affected by persistence from bright targets observed before these observations; the images in the other two filters (J and H) are not affected by this persistence problem. The magnitude limit of the cleaned Y-band image is 29.0, at 5 sigma in  ~0.11 arcsec². The GOODS-ERS and HUDF datasets will be described in detail in Koekemoer et al. (in prep.).

2.3. Photometry

The photometric catalog of the ERS and HUDF fields have been derived in a consistent way. Galaxies have been detected in the J₁₂₅ band, and their total magnitudes have been computed using the MAG_BEST of SExtractor (Bertin & Arnouts 1996), for galaxies more extended than 0.11 arcsec², while circular aperture photometry (diameter equal to 2 times the WFC3 FWHM) is used for smaller sources. The magnitudes which are computed in a circular aperture have been corrected to total magnitude by applying a correction of 0.4 mag. Colors in BVIZYH have been measured running SExtractor in dual image mode, using isophotal magnitudes for all the galaxies, independently from their area. Since the FWHM of the ACS images (0.12″) is better than the WFC3 ones, we smoothed the ACS bands with an appropriate kernel to reproduce the resolution of the NIR WFC3 images. This ensures both precise colors for extended objects and non-biased photometry for faint sources.

As a comparison, we produced, only for the HUDF field, an equivalent BVIZYJH multiwavelength catalog, using the SExtractor aperture magnitudes instead of the isophotal magnitudes to derive galaxy colors. In the resulting Z − Y vs. Y − J plot, indeed, galaxies are much more scattered than in Fig. 1, and the bulk of low-z galaxies are scattered to Z − Y ~ 1 when the aperture colors have been used, while it is limited to Z − Y ~ 0.5 if isophotal colors are adopted. This is due to the fact that the resulting isophotal areas are smaller than the circular apertures adopted, and thus are less noisy.

3. Selection of z = 7 candidates: color criteria

The selection of galaxies at z ~ 7 uses the well known “drop-out” or “Lyman-break” technique. At 6.5 < z < 7.5, this feature is sampled by the large Z − Y color, as shown in Fig. 1 for the HUDF field.

Fig. 1 The Z − Y vs. Y − J color − color diagram for the HUDF dataset limited at J < 29.2. Big points are observed galaxies, small black points are CB07 models at 0.0 < E(B − V) < 0.3, with small red points those at z > 6.5. The small blue points are the result of our simulations, with synthetic objects inserted in the FITS images and recovered using the same procedure used for the observed catalog. The solid line is the color cut of the A criterion, while the dashed line shows the B criterion (described in Sect. 3). Candidate galaxies at z ~ 7 selected with criterion C are shown in magenta. Vertical arrows are upper limits in the Z band, while the blue circle at Z − Y > 3.6 is the suspected SNa outlined by Oesch et al. (2010a).

Recently, a number of slightly different color criteria have been used to select bona fide galaxies at z ~ 7, rejecting lower-z objects and cold stars entering the z-dropout criteria due to photometric scatter. Despite the differences in these criteria, however, the results obtained i.e. by C10a, C10b, Wilkins et al. (2010), Bunker et al. (2010), Oesch et al. (2010a) are comparable, at least at the bright luminosities. We aim here at applying the different color criteria to check whether the LF, which is derived from the data after detailed simulations, is robust against all the random and systematic effects that can affect each selection criterion.

With WFC3 data (ERS and HUDF), we explore three different color criteria, one of Oesch et al. (2010a), one of Bunker et al. (2010), and one proposed by us:

• A) The Oesch et al.(2010a) criterion for the HUDF is:For the ERS, not analysed in Oesch et al. (2010a), an equivalent criterion has been adopted here, that takes into account the difference between the Y₀₉₈ filter instead of the Y₁₀₅ one: For the HUDF field, the S/N(J) > 5 criterion corresponds to J < 29.0, while for the ERS field it allows the selection of galaxies brighter than J = 27.4.

• B) The selection criterion of Bunker et al. (2010) is: The samples of Bunker et al. (2010) and Wilkins et al. (2010; 2011) are limited to Y = 28.5 for the HUDF and Y = 27.2 for the ERS field, respectively, aiming at a more robust selection of z ~ 7 candidates. We adopt the same magnitude limits in this case.

• C) We propose a modification of the Oesch et al. (2010a) criterion, based on simple considerations. We adopt the same color criteria: for the HUDF and for the ERS. To avoid contamination of faint spurious sources, we add the criterion J₁₂₅ − H₁₆₀ ≥ −1.0. This additional criterion removes only two extremely faint galaxies with 29.0 < J₁₂₅ < 29.2 mag, which are detected only in J but not in the H band. Indeed, the observed J − H color for our z-dropout candidates shows a bimodal distribution, with a main peak at J₁₂₅ − H₁₆₀ ~ −0.2 and two isolated sources at J₁₂₅ − H₁₆₀ ~ −1.2. Since LBGs are expected to be almost flat at λ ~ 2000 Å rest frame, we interpret these source with J₁₂₅ − H₁₆₀ ≤ −1.0 as spurious contaminants and remove them from the selected sample. As we have verified in the following (see Sect. 5.1), this further selection does not change the resulting luminosity function at z ~ 7. For the non-detection in bands bluer than Z, we adopt the same criteria used in Castellano et al. (2010a; 2010b) (S/N < 2 in all BVI bands and S/N < 1 in at least two of them). We apply however only a single cut in the J magnitude (J < 29.2 for the HUDF and J < 27.4 for the ERS), since this band corresponds to 1500 Å rest frame. We avoid the double selection S/N(J) > 5 and S/N(Y) > 5, as adopted by Oesch et al. (2010a), since it could be less sensitive to red galaxies in the Y − J color at faint J magnitudes.

Figure 1 shows the Z − Y vs. Y − J color plot for the HUDF field. The solid segments show the A and C criteria, while the criterion B is indicated by the dashed line.

We select 13, 6, and 20 candidates following the A, B, and C criteria in the HUDF field, while in the ERS we find 14, 6, and 22 z-dropout candidates, respectively. The small numbers for the candidates provided by the B criterion reflects the more conservative cut adopted by Bunker et al. (2010) to get rid of the lower-z contaminants with z − Y ~ 1. We have excluded from these candidates the bright transient (probably a SNa), outlined by Oesch et al. (2010a) in the HUDF field (blue circle in Fig. 1 at Z − Y > 3.6). This relatively bright transient was also identified by McLure et al. (2010) and Bunker et al. (2010), while Yan et al. (2011) showed its possible light curve.

Tables 1 and 2 summarize the properties of the z-dropout candidates in the HUDF and ERS field, respectively.

The number counts of the selected z > 7 candidates are different for the three different criteria. In the next paragraphs we will check if the simulations, by including the different selection effects, are able to correct for these differences, and if the resulting Luminosity Function depends on the color criteria adopted.

4. Simulations

While the selection criteria described above are formally designed to select a pure sample of high-z candidates, they are in practice applied to very faint objects, typically close to the limiting depth of the images. At these limits, systematics may significantly affect their detection and the accurate estimate of their large color terms. To take into account all the systematic effects (completeness, photometric scatter) involved in the LF estimate, we carried out a set of detailed simulations, with the aim of deriving the LF.

4.1. The Stepwise method to compute the LF

This method is similar to that used in Bouwens et al. (2008), improved to take into account the photometric scatter, as described in detail in C10a and C10b. We first produce a set of (M₁₅₀₀,z) according to a flat distribution both in redshift and in absolute magnitudes, spanning the redshift range 5.5 < z < 8.5 and the absolute magnitude interval −25 ≤ M₁₅₀₀ ≤ −16. We then convert the absolute magnitudes into observed magnitudes convolving the synthetic SEDs of the LBGs through the relevant filter response curves. For this purpose we have used the models of Bruzual (2007) (hereafter CB07) with the following range of free parameters: metallicity: 0.02, 0.2, and 1 Z_⊙; age from 0.01 Gyr to the maximal age of the Universe at a given z; E(B − V) from 0 to 0.3 following a Calzetti et al. (2000) extinction law. Since Lyman-α emission has an important influence on the selection function of Lyman break galaxies, as shown in lower redshift samples (e.g. Stanway et al. 2008; Dow-Hygelund et al. 2007), we explicitly take into account its effect by considering a Gaussian distribution of Ly-α rest-frame equivalent width with FWHM of 30 Å and centered at EW = 0 Å. We have also added the intergalactic absorption using the average evolution as in Madau (1995).

These synthetic galaxies are finally placed at random positions on the real images. We then recover the measured magnitudes in the synthetic images using SExtractor with the same parameters adopted for the real images. We simulate all the available bands, i.e. from the B to the H bands for the ERS and HUDF fields. To avoid an excessive and unphysical crowding in the simulated images, we have included only 200 objects of the same flux and morphology each time, after masking the regions of the images where real objects have been detected. We repeated the simulation until a total of at least 2  ×  10⁴ objects were tested for each WFC3 field and for each morphological template adopted. These simulations provide the transfer function between the input absolute magnitudes and the observed colors. In addition, the simulations can be used to evaluate the uncertainties in the estimate of the color criteria that we shall exploit to detect z > 6.5 candidates, especially in the Z − Y and Y − J combination, which is essential for a clear identification of z ~ 7 galaxies. These simulations indeed can be used to estimate the impact of different systematics in the study of the LF. For example, since the Y band filter used in the ERS field (Y₀₉₈) has a bluer effective wavelength than the one in the HUDF field (Y₁₀₅), it turns out that the selection function in redshift in these two fields is different. These systematics have been taken fully into account with these simulations.

The simulations described above have been used to derive an estimate of the z ~ 7 LF with the Stepwise method. This method assumes that the rest-frame luminosity function of galaxies can be approximated by a binned distribution, where the number density φ_i in each bin is a free parameter: this non-parametric approach allows us to constrain the number density of galaxies at different magnitudes without assuming an a priori functional form (i.e. a Schechter-like shape).

We have assumed that the LF is made of four bins in the interval −22.5 < M₁₅₀₀ < −17.0, corresponding to the range sampled by our observations. We also assume that galaxies are uniformly distributed within the bins, with number densities φ_i to be determined. With these simulated galaxies, we compute the distribution of observed magnitudes originated from each bin, scaled to the observed area in the WFC3 fields (subtracting the area covered by detected objects). We then find the combination of φ_i that best reproduces the magnitude distribution of our observed objects with a simple χ² minimization. Since the photometric scatter can move galaxies from a magnitude bin to the near ones, the χ² minimization ensures that this effect is taken into account. The resulting uncertainties can be higher than those derived assuming a simple linear relation between the apparent and absolute magnitude, as in Bouwens et al. (2008). These differences are discussed in detail in Castellano et al. (2010b).

Finally, the simulations are used to estimate the systematic effects acting when we use colors at shorter wavelengths, i.e. in the BVI bands, to reject possible interlopers. Because of the large IGM and internal HI absorption, the expected flux in these bands for z > 6.5 galaxies is far below the detection threshold, or even null. For this reason a stringent limit on the measured flux in these bands is adopted to remove lower redshift interlopers. However, the S/N estimated by SExtractor may be a poor representation of the actual photometric scatter at low fluxes, due to a combination of factors, such as uncertainties in the estimate of the local background, underestimates of the true rms, or chance superposition of faint blue galaxies along the line of sight. To account for these effects, we have measured the resulting signal-to-noise (S/N) ratio in the BVI images for each simulated object inserted in the J one, which should be zero on average. It turns out that the actual distribution of the S/N ratios is wider than the one obtained with SExtractor, which is computed scaling the input weight image. We thus estimate the “effective” rms σ_S/N, i.e. the rms of the signal-to-noise distribution in each of the BVI images, which is about 1.5 for the B and V bands and 1.8 for the I band. Even taking into account this wider distribution, we find that the tails of the S/N distribution (S/N > 2σ_S/N) contain more objects than in the case of a pure Gaussian distribution. As mentioned above, we will use the estimated σ_S/N in all BVI bands, requiring that high-z candidates have flux  <2σ_S/N in all BVI bands and  <1σ_S/N in at least two of them. With our simulation, we estimate that the fraction of true high-z galaxies lost because of this strict criterion is about  ~12% for the HUDF field. For the ERS, we find an effective S/N ratio of 1.8 for all the BVI bands and a rejection of true high-z galaxies of  ~24%. Such an aggressive requirement is needed in order to clean our sample by low-z interlopers that are contaminating our z-dropout sample. This technique has been extensively discussed in C10a and C10b. In particular, if we neglect the requirement that S/N must be  <1σ_S/N in at least two bands out of the BVI ones, we find 32 and 97 galaxy candidates for the HUDF and ERS fields, respectively, to be compared with our 20 and 22 master candidates. We thus prefer to reduce significantly the numbers of possible contamination by interlopers, at the expenses of being less complete by 24%, especially for the ERS field. This fraction is significantly higher than that found in the HUDF and it is probably due to the fact that the depth of the ACS GOODS images is not matched to the deep limiting magnitude one can reach with WFC3 in the ERS. In this sense, these simulations are essential to estimate the LF at z ~ 7.

4.2. The morphological profiles used during simulations

We expect that the output of the simulations depends critically on the assumed morphology. As realistic morphologies for the simulations related to the WFC3 fields, we check different profiles. In particular, we explored the following morphologies spanning from the most compact to the most extended objects available in the WFC3 fields, exploring also plausible extrapolations from lower-z studies:

• 1) We select four point-like objects in theHUDF field based on their FWHM, avoidingbright saturated objects. The spectroscopic classification is notavailable for these objects, and their colors are consistent withthose of stars. Since these are compact, we use them as an inputprofile for very compact galaxies at z = 7, as it is expected that faint galaxies at high-z are very small (Windhorst et al. 2002).

• 2) The transient object at Z − Y > 3.6 in Fig. 1 (probably a SNa) is used as an empirical profile for our simulations. Since this object was not detected in previous NICMOS images of HUDF (Oesch et al. 2010a), it is plausible that the host galaxy in the background is not able to change the SNa profile, and the result is a very compact object. We used this template to check whether it gives different results from the four stars described above.

• 3) We select the three brightest and most compact z ~ 7 candidates in the HUDF (J < 27.4) as input for our simulations, assuming that the physical dimension of high-z galaxies does not vary dramatically with their magnitudes. We avoid in this group the object HUDF-1095 (G2_1408 of C10a) since it is the most extended z ~ 7 candidate in our sample and we consider it in a separate group.

• 4) We study the effect of very extended objects or of merging at high-z using as input in our simulations the object HUDF-1095 (G2_1408 of C10a). It is very bright with two knots, possibly the result of a merging of two galaxies. Alternatively, its clumpy morphology in the rest frame UV could be due to extended star forming knots in a nascent galaxy. First results based on spectroscopic confirmation with FORS2 seem to indicate that it is a z = 6.972 LBG (Fontana et al. 2010), confirming the photometric redshifts from ACS and WFC3 based SED. Though not representative of the whole z = 7 population, this galaxy is an LBG at z ~ 7 and it is interesting to study the resulting LF assuming this morphology as input for our simulations.

• 5) There is the possibility that fainter galaxies are more compact than brighter ones, as expected in the hierarchical formation models. To check this hypothesis, we adopt as input profile the one obtained by stacking the 15 fainter galaxy candidates at z = 7 (27.5 < J < 29.2) in the HUDF field.

• 6) We used as input for our simulations the profiles of spectroscopically confirmed galaxies at z ~ 4 in the ACS I-band of the HUDF, scaled as (1 + z)^-1 in their physical sizes, as adopted in Bouwens et al. (2008, 2010a). Since there is a regular size evolution with redshift, as found by Ferguson et al. (2004) and confirmed by Bouwens et al. (2004) and Oesch et al. (2010b), it is useful to explore this case, in order to compare our results with earlier LFs that have appeared in the literature. These templates are also useful to avoid using self-selected z = 7 candidates as input for our simulations.

Figure 2 shows all the model galaxies used for the HUDF and ERS simulations.

Fig. 2 The different morphologies used for the simulations of WFC3 fields. From top left, the first four cutouts represent non-saturated stars in the HUDF survey (morphological template 1). The SNa candidate outlined by Oesch et al. (2010a) (template 2) is then shown, followed by the three brightest z ~ 7 candidates of the HUDF survey (template 3). The most extended galaxy is the brightest candidate in the HUDF (HUDF-1095 or G2_1408 of Castellano et al. 2010a, template 4). The first thumbnail on the fourth raw shows the stack of the faint candidates selected in the HUDF (template 5). All these thumbnails have been extracted from the J125 band of HUDF. The remaining objects are spectroscopically confirmed galaxies z ~ 4, extracted from the ACS-I band of the HUDF and shifted at z = 7 as described in the main text (template 6), and two artificial exponential disk galaxies with half light radius =0.3 arcsec (templates 7 and 8) and with half light radius =0.1 arcsec (templates 9 and 10).

In addition of adopting single objects for the simulations, we explore also different distributions of half light radii (rh) at z ~ 7:

• 7) Following Oesch et al. (2007), we explore auniform distribution in the half light radius of the syntheticgalaxies, assuming an exponential declining profile. The twoextremes of the uniform distribution are derived from the valuesused by Oesch et al. (2007) for their simulations torecover the z = 5 LF from HUDF data. They adopt a uniform distribution with half light radius from 0.05 to 0.5 arcsec. Converting these values to physical units at z = 5, and applying a (1 + z)^-1 evolution until z = 7, and reconverting in angular dimensions, we derive the limits to be 0.045 to 0.45 arcsec. From Fig. 3, it is clear that this distribution can be ruled out by the present observations. However, we decided to study this uniform distribution as a “thought experiment” on how hypothetical size distributions might change the derived completeness at faint magnitudes, but we do not use it to estimate the z ~ 7 luminosity function.

• 8) Another possible distribution for the half light radius of the simulated galaxies is the log-normal distribution with mean and dispersion adopted following the relation of Ferguson et al. (2004) for photometrically selected LBGs at z = 4, applying a (1 + z)^-1 evolution to the half light radii in physical units. We derive for the z = 7 population a mean half light radius of 0.258 arcsec and σ = 0.0358. This approach has been adopted also in Oesch et al. (2009 and 2010a). In the following it is indicated also as M1.

• 9) The typical half light radii of z = 7 candidates are of the order of 0.1−0.15 arcsec. As in Oesch et al. (2010b), we assume a log-normal distribution of the half light radius of the simulated galaxies with mean =0.153 and sigma = 0.0574 arcsec. This distribution, however, cannot explain the presence of the HUDF-1095 galaxy, with rh > 0.25 arcsec. In the following it is indicated also as M2.

• 10) A different distribution for the half light radius of the simulated galaxies at z ~ 7 has been derived starting from the observed half light radii of spectroscopically confirmed LBGs at z = 4 in GOODS from Vanzella et al. (2009). We select only galaxies with 24 ≤ I ≤ 25 magnitudes to avoid incompleteness for extended LBGs. Applying a (1 + z)^-1 evolution to the half light radii at z ~ 4 in physical units, the log-normal distribution has a mean =0.177 and sigma = 0.0751 arcsec at z ~ 7 and in apparent dimensions. In the following it is indicated also as M3.

From Fig. 3 it is clear that the observed distribution of half light radii for the z = 7 candidates peaks at  ≃0.1 arcsec and extends to 0.2−0.25 arcsec. The small area covered by the two WFC3 surveys considered here is at present too small to sample the extended tail of the half light radii distribution, and it is not possible to infer a robust size distribution. Clearly, the uniform distribution (criterion 7) is not representative of the observed z = 7 galaxy population and will not be considered in the following to estimate a reliable z ~ 7 luminosity function.

Fig. 3 The black histogram shows the observed distribution of the half light radii (in arcsec) of the z ~ 7 galaxy candidates brighter than J = 27.8 in the HUDF. At this magnitude limit the completeness of the J band image is higher than 80% even for the most extended objects. The galaxy with half light radius  ~0.25 arcsec is the HUDF-1095 one (template 4). The cyan solid line shows the uniform distribution in half light radius (template 7), while the dark pink and green solid lines indicate the two log-normal distributions in half light radius adopted by Ferguson et al. (2004) (template 8, M1) and by Oesch et al. (2010b, template 9, M2), respectively. The purple line indicates the distribution derived from Vanzella et al. (2009) (template 10, M3). All sizes are corrected for PSF. The dashed line shows the completeness of galaxy detection for different half-light radii derived by the performed simulations. This curve has been multiplied by a factor of 2 only for clarity.

The three log-normal distributions (templates 8, 9, and 10 shown in Fig. 3) cannot be confirmed or ruled out at the present stage for a number of reasons:

• a) The z ~ 7 candidate sample at magnitude J < 27.8 in the HUDF is limited to 6 galaxies. We have chosen this magnitude limit since it is known from simulations that in the HUDF we are complete at the  ~80% level for extended galaxies (see Fig. 4, blue and cyan lines). The corresponding limits for shallower surveys (ERS or the parallel HUDF observations HUDF-P1 and HUDF-P2) are 1−2 mag brighter (J < 26.0−26.8) than in the HUDF. Adopting these limits, only two additional brighter candidates have been found on the ERS survey, both having rh > 0.2   arcsec. Since the number statistic here is very low due to the limited area covered by ultradeep HST observations, at the present stage we cannot exclude the presence of more extended galaxies at J < 27.8 mag.

• b) The extrapolation to fainter magnitudes of the size distribution for z ~ 7 galaxies is not known at present. It is possible that a population of extended sources (similar to HUDF-1095) exists at high-z and at fainter magnitudes. The present HUDF observations cannot rule out such hypothesis.

• c) The spectroscopic confirmation of z ~ 7 galaxy candidates with detailed HST morphology and half light radius determination in the HUDF and ERS fields is still lacking at this stage and some of them could be lower-z contaminants. The brightest galaxy in the HUDF (HUDF-1095, which is clumpy or in a merging state) turned out to be spectroscopically confirmed at z = 6.972 (Fontana et al. 2010). Since it is the only galaxy with detailed HST morphology in the NUV rest frame (template 4) to be confirmed at z ~ 7 in the HUDF field, and it is relatively bright (J = 26.07) and extended (rh = 0.256), we cannot exclude that a similar population of extended sources exists at high-z and at fainter magnitudes.

• d) In Bouwens et al. (2011), and also in Wilkins et al. (2011), a z ~ 7 galaxy candidate with an half-light radius of 0.46 arcsec has been found (UDF092z-01191133). Moreover, other z ~ 7 candidates found in the ERS, HUDF-P1, and HUDF-P2 fields show extended half light radii (0.3 arcsec or more). If they turn out to be spectroscopically confirmed at z ~ 7, then distributions with extended half light radii, similar to template 8 (M1), would become more reliable.

Based on these considerations, it is not possible at the present stage to exclude any of the adopted log-normal half light radii distributions (M1, M2 or M3).

The simulations carried out with the different morphological templates described above have been used to estimate the completeness in the detection procedure. This is shown in Fig. 4, where we plot the fraction of detected versus input objects as a function of the J₁₂₅ input magnitude for the HUDF field.

Fig. 4 The fraction of detected galaxies versus the J125 input magnitude for the simulations in the HUDF field. Different colors indicate the various morphologies adopted: the dashed magenta lines indicate non-saturated stars on the HUDF survey (morphological template 1). The solid magenta line shows the completeness obtained using the SNa candidate outlined by Oesch et al. (2010a) (template 2) as morphological template, while the red solid curves show the three brightest z ~ 7 candidates of the HUDF survey (template 3). The blue solid line is the completeness derived by adopting as template the brightest candidate in the HUDF (HUDF-1095 or G2_1408 of Castellano et al. 2010a, template 4), and the black solid line represents the completeness obtained adopting as input template the stack of the faint candidates selected in the HUDF (template 5). The green dashed lines indicate the spectroscopically confirmed galaxies z ~ 4, shifted at z = 7 (template 6). The cyan solid line shows the completeness obtained assuming a uniform distribution in half light radius (template 7), while the dark pink and green solid lines indicates the results of two log-normal distributions in half light radius adopted by Ferguson et al. (2004), Oesch et al. (2007, 2010a) (template 8, M1) and by Oesch et al. (2010b) (template 9, M2), respectively. The purple line indicates the completeness derived from Vanzella et al. (2009) (template 10, M3).

The completeness for the HUDF is basically equal to 100% down to J ≃ 27 for all the adopted morphologies, and fades in different ways according to the different half light radii of the templates adopted in the simulation. In particular, more extended galaxies typically have a lower completeness at J > 27, and the fraction of faint z = 7 galaxies recovered depends dramatically on the morphological templates used during the simulations. Since the completeness is used to estimate the LF, we expect the faint end of the z = 7 LF to be strongly dependent on the morphological template adopted, as we will show in the next section.

5. Results

The derived LF is sensitive to the details of the simulations adopted. Here we are interested in exploring mainly two issues: the role of the color selection criteria on the estimate of the LF and the role of galaxy morphology on the measurement of α, the steepness of the faint end of the LF.

Fig. 5 The LF at z = 7, obtained using as morphological templates real bright (J < 27.4) candidates at z ~ 7 found in the HUDF (template 3). The magenta circle represents the LF found in the Hawk-I from C10a, while the other symbols indicate the LF in the HUDF and ERS derived adopting different color criteria: our color criteria (red filled circles, criterion C), Bunker et al. (2010) (small filled black triangles, criterion B) and Oesch et al. (2010a) (blue squares, criterion A). The blue points are shifted by −0.1 mag and the black points by +0.1 for clarity. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) and Bouwens et al. (2007), respectively.

5.1. Results for different color selections

We first explore three different color criteria for the selection of z ~ 7 galaxies. The LF found with these three selection criteria (A−C), described in Sect. 3, for the HUDF and ERS fields is shown in Fig. 5. The red filled circles show the LF derived with the color criterion C, compared to the one derived with the B (small filled black triangles) and A criteria (blue squares), for the HUDF data alone. We plot for comparison also the Hawk-I results (magenta open circle) of C10a. The WFC3 LF has been derived assuming as morphological inputs for our simulations real bright (J < 27.5) candidates at z ~ 7 found in the HUDF (template 3). Adopting as morphological templates spectroscopically confirmed galaxies at z ~ 4 found in the ACS i-band of the HUDF and scaled to z ~ 7 (template 6), or the stack of the fainter (J > 27.5) candidates (template 5), we obtained very similar results. We have also verified that the additional limit J₁₂₅−H₁₆₀ ≥ −1.0 adopted in criterion C does not imply significant changes in the LF in the faintest bin (M_UV ~ −18.5).

Fig. 6 The LF at z = 7, obtained using as morphological types galaxies drawn from the templates 1−6. The color code is the same of Fig. 4. The blue and green bins are shifted by −0.2 and −0.1 mag, while the red and black points are moved by +0.1 and +0.2, respectively, for clarity. The blue dotted lines indicates the range of variability for the best fit of α spanning the range −2.15 ≤ α ≤ −1.38. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) and Bouwens et al. (2007), respectively. The filled black squares are the LF presented in Oesch et al. (2010a).

It is clear that the three color selection criteria analysed here give comparable results, well within the uncertainties, despite the size of the three samples changes by a factor of  ~ 3. This is due to the fact that all the three criteria are well designed to isolate the few z ~ 7 candidates from the great number of interlopers, but each method makes different assumptions, resulting in a variable number of z-dropout galaxies. The resulting LFs, however, are similar, since the simulations required to derive them are able to take into account properly the random and systematics effects (incompleteness and contamination by photometric scatter), especially at the faintest magnitudes. In the faintest magnitude bin the sizes of the error bars depend on the criteria adopted, with the larger errors for the most strict criterion. The B criterion selects indeed a smaller number of candidates, with a resulting higher Poisson noise than the one derived adopting the A or C criteria.

The only difference between the A and C criteria is the magnitude cut adopted: we select all galaxies down to J = 29.2, while Oesch et al. (2010a) recover only galaxies with S/N greater than 5 in the Y and J bands. These color criteria are similar, and no systematic differences are found in the LF estimate. For this reason, we decide to adopt the C color criterion in the following sections.

5.2. Results for different morphologies

We then explore the issue of adopting different morphological templates in the simulations. For each morphological type shown in Fig. 2 and discussed in Sect. 4, we have simulated 2  ×  10⁴ synthetic galaxies as described in Sect. 4 and used them to compute the LF at z ~ 7. In Fig. 6 we show the LF for each template considered: at M₁₅₀₀ ~ −20, the different morphologies adopted have a small impact on the number density Φ, while at M₁₅₀₀ ~ −18 the LF is extremely dependent on the assumption of the morphological templates. In Fig. 6 the open points show the result of templates 1−6, based on single object morphologies, with the relevant uncertainties.

In Fig. 7 we considered different distributions for the simulated half light radii, namely the uniform one (template 7, cyan points), and the three log-normal distributions (template 8 using dark pink points, template 9 using dark green points, and template 10 using purple points). For comparison, we plot the LF derived by Oesch et al. (2010a) as filled black squares: these are in agreement with the result based on the simulations carried out with the distributions of half light radii which extend to 0.2−0.4 arcsec. The large variance in the galaxy number density at the faint end of the z ~ 7 LF depends on the different completeness corrections applied to the same dataset and thus reflects the uncertainties due to the poorly known half light radius distribution.

Fig. 7 The LF at z = 7, obtained using different distributions for the simulated half light radii, namely the uniform one (template 7, cyan points), and the three log-normal distributions (template 8 using dark pink points, template 9 using dark green points, and template 10 using purple points). The dark green bins are shifted by +0.1 mag, while the dark pink points are moved by −0.1 mag, respectively, for clarity. The blue dotted lines indicate the range of variability for the best fit of α spanning the range −2.15 ≤ α ≤ −1.38. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) and Bouwens et al. (2007), respectively. The filled black squares are the LF presented in Oesch et al. (2010a).

Fixing the Schechter parameters Φ ^∗  and M ^∗  to 6.9  ×  10^-4 and −20.10 (the best fit values of Ouchi et al. 2009), respectively and varying α, we derive a large range of variation for the best fit of the LF steepness, namely −2.15 ≤ α ≤ −1.38 at 68% c.l., represented by the two blue dotted lines in the plot.

The steepness of the faint end of the LF, α, depends critically on the half light radii of the synthetic galaxies used to carry out the simulations. In Fig. 8 we show the half light radii of the different morphological templates against the best fit α of the LF. For large value of the half light radius, namely  >0.2 arcsec, the best fit α is between −1.8 and −1.9, as found by Oesch et al. (2010a), Bouwens et al. (2010a, 2011) at z = 7 and beyond. For smaller values,  ~ 0.1 arcsec, we find α in the range −1.4 to −1.7, in marginal agreement with McLure et al. (2010) that have used point sources as input for their simulations. The uncertainties on this parameter are still large, of the order of  ~ 0.2−0.4. We find however a significant relation between the half light radius and α, indicating that simple conclusions based on the LF at z ~ 7 depend dramatically on the half light radius distribution assumed for the simulations.

Fig. 8 The dependence of the steepness of the LF, α, against the half light radii of the simulated galaxies. Different colors indicate various morphology adopted, using the same color code of Fig. 4.

5.3. Fitting the LF at z  ~  7

We want to derive here the LF at z ~ 7 making reliable and robust assumptions based on the considerations outlined above. We adopt the color criterion C, with a magnitude cut at J = 27.4 and J = 29.2 for the ERS and the HUDF, respectively. As shown above, the simulations are able to correct the incompleteness of this color criterion in a reliable way. The morphology templates used in the simulations, however, are critical, since the simulations cannot correct for a biased assumption on the profiles of the candidates. Since it is not possible, with the small number statistics achieved by the present WFC3 data, to derive precisely the distribution of the half light radii of the z ~ 7 LBGs, especially at the faint end and at high values of rh, we carry out the fit for the two most extreme, but reliable, variants of the LFs derived in the previous paragraph. We have adopted the simulations carried out with the stacking of the faint candidates selected in the J band of HUDF (template 5) as an extreme case (all the z = 7 galaxies are compact), while the HUDF-1095 galaxy (template 4) has been chosen as the opposite case (all high-z LBGs are very extended). We have verified that all the other options (templates 1−3 and 6−10) give results that are in between these two extremes.

We have combined the results of the ERS and HUDF fields in a single LF, to reduce the statistical uncertainties and to have an extended range in the absolute magnitude M₁₅₀₀. To cover the bright part of the LF, we also use the results of Ouchi et al. (2009) obtained from a shallow and wide Subaru survey and those of C10a and C10b derived from an intermediate depth survey with Hawk-I. In Fig. 9 we show our results (red and green filled squares) complemented with the results of Ouchi et al. (2009) (filled black circles), Bouwens et al. (2010a) (filled cyan circles) and of C10a and C10b (empty black circles).

Fig. 9 The LF at z = 7, obtained using as morphological templates the stacking of the faint candidates at z ~ 7 found in the J band image of the HUDF (template 5) and the C color criterion. The LF is computed using ERS and HUDF fields together (red filled squares, G10). The green filled squares show the LF obtained with the same color criterion C but using the HUDF-1095 candidate as input for the simulations (template 4). These data have been complemented with the results of Ouchi et al. (2009) (O09), Bouwens et al. (2010b) (B10) and C10b. The lower dashed magenta line is the best fit to all the observed points related to the template 5, the upper dashed magenta line is the best fit related to the template 4, while the black solid line shows the Hawk-I result of C10a. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) (M09) and Bouwens et al. (2007) (B07), respectively.

All the observed points shown in Fig. 9 have been used to fit the LF. We scan the parameter space in 3D (Φ^∗, M^∗, and α) adopting a Schechter parameterization (Schechter 1976) and taking into account the asymmetric errors and the presence of upper limits. The best fit is shown (long dashed magenta lines), and it is compared to the early Hawk-I result of C10a (black solid line) and the LFs at z = 6 of Bouwens et al. (2007) and McLure et al. (2009). Using our best fit, we provide a robust estimate of the LF at z = 7, excluding at  >99.73% c.l. (>3σ) the non-evolution from z = 6 to z = 7. Indeed, we are able to show that the normalization at M₁₅₀₀ ~ −19 is significantly lower than the LF at z ~ 6. Moreover, combining data with a large range in absolute magnitudes we are able to reduce the degeneracies in the three parameters of the LF. In Fig. 10 we provide the confidence regions at 68% c.l. for Φ^∗, M^∗, and α at z = 7, exploring also the constraints given on the SFR density (SFRD).

Fig. 10 The confidence regions at 68% c.l. (1σ) for the LF parameters and the SFRD. The red points show the projection of the best fit. The plot on the left shows the confidence regions based on the simulations with compact morphology (template 5), while the right plot is related to extended morphology (template 4).

Table 3 summarizes the best fit parameters for the Schechter LFs, with their 68% c.l. uncertainties. The luminosity density ρ_UV has been computed integrating the fitted LFs down to an absolute magnitude (at 1500 Å rest frame) of M_UV = −10 and to M_UV = −19. We have converted these values in a SFRD following the standard formula by Madau et al. (1998) and applying the extinction correction of Meurer et al. (1999) of A₁₆₀₀ = 0.43 (considering an average UV slope β = −2.0 for the SED of the z = 7 LBG candidates). Our estimates of the luminosity density and SFRD are consistent with other values in the literature, reducing significantly their uncertainties.

6. Discussion

We discuss the implications on reionization derived from our best fit of the LF and the uncertainties on the LF parameters implied by the assumption made during completeness simulations.

6.1. Reionization

Following Bolton & Haehnelt (2007) and Ouchi et al. (2009), the contribution of z  =  7 LBGs to the hydrogen ionizing photon can be set as (1)where α_S is the spectral index of ionizing emission and f_esc is the escape fraction of ionizing photons. The rate of ionizing photons needed to balance the recombination process of hydrogen in the IGM and hence to keep the Universe reionized is Ṅ_ion ≥ 10^47.4C_HII(1 + z)³, that translates into a constraint to the clumpiness of the IGM to be C_HII ≤ (ρ_UVf_esc)/(10²²α_S(1 + z)³). We fix these parameters to z = 7, β = −2, and α_S = 3 in the following. In Fig. 11 we show in the dashed areas the C_HII allowed by a given f_esc to keep the IGM ionized at z ~ 7: the red region is for a compact morphology (template 5), while the blue one refers to an extended morphology (template 4). The relevant UV emissivity of LBGs ρ_UV is computed by integrating the LF down to M₁₅₀₀ = −10 (Bouwens et al. 2011), assuming that the steepness of the faint end of the LF remains constant; we have no information on the number density of z ~ 7 galaxies at these faint magnitudes with the present data. Given our confidence intervals for the LF parameters, and hence the limit to ρ_UV(−10), we derive a limit C_HII ≤ 6.5 (at 68% c.l.) in order to have the Universe ionized at z = 7, assuming a maximum escape fraction of 1.0, and the simulations based on template 5 (compact morphology), as shown in Fig. 11. If we use the LF based on extended morphology (template 4), the 68% c.l. limit to the clumpiness of the ionized hydrogen is C_HII ≤ 29.6. In this reasoning we must consider C_HII ≥ 1, since galaxies at z ~ 7 are formed in biased density regions of the Universe and the IGM is plausibly not homogeneous (C_HII = 1) at these redshifts. Limiting the clumpiness to C_HII ≥ 1, we thus have f_esc ≥ 0.17 in order to have the Universe reionized by z = 7 assuming compact morphology, while the limit is f_esc ≥ 0.034 for extended morphology.

Fig. 11 The constraints on the maximum CHII and minimum fesc parameters in order to have the Universe reionized by LBGs at z = 7 given our LF. Given a value for fesc, the Universe is reionized if CHII is lower than the red line, if a compact morphology is used (template 5), or lower than the blue line if an extended one is considered (template 4). The ρUV is computed integrating the LF down to M1500 = −10 mag. The solid line at CHII = 1 indicates a homogeneous IGM, and we must consider only the CHII ≥ 1 region, since galaxies at z ~ 7 are formed in biased density regions of the Universe. The dashed line shows the value CHII = 3. In practice, only a combination of (fesc, CHII) parameters inside the shaded regions corresponds to a z ~ 7 IGM kept ionized by LBGs only.

Bolton & Haehnelt (2007) have inferred C_HII ≤ 3 at z ~ 6 from the Lyman-α forest photoionization state: since C_HII is expected to monotonically decrease towards high-z in a hierarchical Universe, an escape fraction of f_esc ~ 0.5 is enough for the reionization of the IGM (but not for the more clumpy galaxy formation regions) considering template 5 (compact morphology), while a value of f_esc ~ 0.1 is sufficient according to template 4 (extended morphology).

We neglect in our computation the contribution of AGNs, since their LF at z = 7 is still unknown and the upper limits currently available indicate that the AGN will add only 5−8% to the luminosity density of galaxies (Ouchi et al. 2009).

The reionization rate Ṅ_ion however depends largely on the contribution of faint galaxies (−19 ≤ M₁₅₀₀ ≤ −10), since ρ_UV(M₁₅₀₀ < −10) could be 10−15 times larger than ρ_UV(M₁₅₀₀ < −19) if the LF is steeper than α ≤ −2. This is eventually still possible, if a large number of faint and extended galaxies exist, similar to the object HUDF-1095 (G2_1408 of Castellano et al. 2010a). If the faint end of the LF at z = 7 is not so steep, as we have inferred here for the compact morphologies, the contribution of fainter galaxies (−19 ≤ M₁₅₀₀ ≤ −10) is only 2.4 times that of M₁₅₀₀ ≤ −19 LBGs. Thus, the present uncertainty on the steepness of the faint end of the LF translates in a large uncertainty, of a factor 2−15, on the best fit of the luminosity density, SFR density, and reionization rate at z = 7.

Following Windhorst et al. (2002), it is plausible that fainter galaxies have smaller physical sizes: in this case a combination of high escape fraction (f_esc ≥ 0.2) and small clumpiness (C_HII ≤ 3) is required for the reionization of the Universe. It is worth noticing that the recent estimates of f_esc for L ≥ L^∗ LBGs at 0 ≤ z ≤ 3 are  ≤0.15 (Bridge et al. 2010; Cowie et al. 2010; Siana et al. 2010; Vanzella et al. 2010), implying a consistent evolution of the galaxy escape fraction going to fainter luminosities or to higher redshifts.

Moreover, the knowledge of the reionization processes in the high-z Universe is currently far from clear. For example, early theoretical works (i.e. Gnedin & Ostriker 1997) have found high values for the clumpiness, C_HII ~ 10−30. More recently, Pawlik et al. (2010) studied the effect of photoionization heating by a uniform ultraviolet background and found that photo-heating can strongly reduce the clumping factor of the IGM because the increased pressure support smoothes out small-scale density fluctuations. The common expectation is that photoionization heating should provide a negative feedback on the re-ionization of the IGM because it suppresses the cosmic SFR by boiling the gas out of low-mass DM halos. However, since it also contribute to the reduction of the clumping factor, the net result is that it is easier to keep the IGM ionized. Photo-heating therefore also provides a positive feedback for the reionization processes. Neglecting the effects of photoionization heating, Pawlik et al. (2010) predict a clumpiness of the IGM of the order of 7−8 at z = 7, and in this case our limit to the LF rule out a completely ionized IGM even for f_esc = 1 if all the galaxies at z = 7 resemble our morphological template 5 (compact galaxies), while an f_esc ~ 0.3 is still sufficient assuming template 4 (extended galaxies). If instead the photoionization heating is considered, a small clumpiness C_HII ≤ 3 is expected, in agreement with measurements by Bolton & Haehnelt (2007). In this case, an escape fraction of  ~ 50% is enough to keep the Universe re-ionized, even in the most pessimistic case (of compact morphology).

6.2. The reliability of our LF estimate

The uncertainties on the estimate of the z = 7 LF are related to the method used to derive it. Errors can be due to the limited depth of the images used, the color criteria adopted, the presence of possible interlopers (cool dwarf stars, low-z red galaxies, transient objects), spurious detections, especially when going faint. All these features have been discussed in the previous sections, showing that detailed simulations are able to take into account the effects of different color criteria, and that a precise color cut can be identified with the aim of separating the few high-z galaxies from the much frequent interlopers. Moreover, the details on the steepness of the faint end of the LF are strongly influenced by the morphological templates adopted for the simulations.

There are however two issues that cannot be corrected with the simulations only: the presence of transient objects within the z-dropout candidates and the role of the cosmic variance when dealing with deep, pencil beam surveys. In the near future e.g., the HUDF field will be re-observed with WFC3 in order to go deeper than the present observations, with 34 additional orbits splitted in three filters (Y, J, and H), complementing the existing 60 orbits discussed here. With these data, it will be possible to assemble a z-dropout sample clean from variable contaminants. The discovery of a clear transient source in the HUDF data (see Fig. 1) and two possible transients in the Hawk-I data by Castellano et al. (2010a) is more than an order of magnitude higher than the 0.06 SNe per WFC3 field estimated in Oesch et al. (2010a). It is possible that some of the selected z = 7 candidates will be rejected as interlopers when a second epoch survey will be executed.

Surveys are usually bounded by a limited amount of telescope time, and thus they have to compute a trade off between the covered area and the imaging depth. This has deep implications on the number of objects recovered, on their absolute magnitudes and the sampled Universe volume. Simple observational estimates of galaxy number counts in finite volumes are thus subject to the uncertainties due to cosmic variance, arising from underlying large-scale density fluctuations. This may be particularly true at high-z, where galaxies are expected to be much more clustered than at lower-z.

We have used the Millennium simulation, using the realization of Kitzbichler & White (2007), to investigate the impact of the cosmic variance on the number of galaxies recovered in the two WFC3 fields analysed here. We extract randomly from the simulation of Kitzbichler & White (2007) 10⁴ pencil beam realizations, recording for each pointing the number count N of galaxies at 6.5 ≤ z ≤ 7.5 down to M₁₅₀₀ = −17.5 and −19.3 for the HUDF and ERS, respectively. The cosmic variance is computed as , where the first term is the measured variance and the second term is the Poissonian contribution. We found that the expected uncertainty due to cosmic variance is σ_CV/ ⟨ N ⟩   =   ~34% and  ~38% for the HUDF and ERS field, respectively. Combining the two fields, the resulting cosmic variance yields  ~28%. These results are in agreement with the predictions based on Somerville et al. (2004), extrapolated to z ~ 7.

As shown in Robertson (2010), since the errors on M^∗, Φ^∗, and α are correlated, this does not translate directly in an additional  ~30% uncertainties in the LF parameters. However, the uncertainties on the luminosity density or on the UV ionizing photons depends mainly on the total number of galaxies, so the  ±30% additional uncertainties can be applied to these quantities to take into account the effects of cosmic variance.

7. Summary

We have analysed two WFC3 fields, the HUDF and ERS surveys. More than 40 bona fide z-dropout galaxy candidates have been selected on these two areas, down to an observed magnitude of J ~ 29 (AB). The depth of ACS images observed with bluer filters (BVI) ensures a robust and well controlled selection of the high-z galaxy sample, cleaning the candidates from the large number of interlopers (i.e. low-z early galaxies, dusty starbursts or cool stars).

We have explored different color criteria for the selection of z ~ 7 galaxies, picking up three methods from the wealth of methods proposed recently in the literature. The numbers of the resulting candidates are very different, with 13, 6, and 20 candidates following the A, B, and C criteria (described in Sect. 3) in the HUDF field, while in the ERS we find 14, 6, and 22 z-dropout candidates, respectively. We have shown here that, applying these different methods both to the observed surveys and to the simulated samples, the resulting LFs are similar.

On the other hand, the output LF depends critically on the morphological templates used to simulate the completeness of the observed candidates. In Sect. 4 we have proposed 10 different options for a reliable simulation, spanning from point like objects to very extended galaxies, both observed and synthetic. The resulting LF is extremely sensitive to the adopted morphological templates, especially the parameter α governing the steepness of the faint end. We have shown that the steep LF found by Bouwens et al. (2011) is possibly due to the choice of a set of templates with typical half light radius larger than the one of the observed candidates at z = 7. If a more compact template is adopted (as used in McLure et al. 2009; and Wilkins et al. 2011), a milder slope is derived.

The resulting best fit to the z ~ 7 LF is , , and . This implies a luminosity density (down to an absolute magnitude of −10) of × 10²⁵   erg/s/Hz/Mpc³ and a star formation rate density of , if a compact morphology is used (template 5). If a more extended morphology is adopted, the resulting best fit to the z ~ 7 LF is , , and . This implies a luminosity density (down to an absolute magnitude of −10) of × 10²⁶   erg/s/Hz/Mpc³ and a star formation rate density of .

The contribution of z = 7 LBGs to the hydrogen ionizing photon density Ṅ_ion depends on the last two parameters, and the threshold to keep the Universe ionized at z ~ 7 depends also on f_esc, the escape fraction of ionizing photons, and on the clumpiness of the IGM, C_HII. Limiting the clumpiness to C_HII ≥ 1, we thus have f_esc ≥ 0.17 in order to have the Universe reionized by z = 7. Inferring a small clumping factor for the IGM at high-z (C_HII ≤ 3 at z > 6, Bolton & Haehnelt 2007), an escape fraction of f_esc ~ 0.5 is enough for the reionization of the IGM (but not for the more clumpy galaxy formation regions), assuming a LF derived from template 5 (compact morphology). If the IGM was reheated at z ≥ 9, the observed population of LBGs at z = 7 might be enough to keep the IGM ionized, given that the escape fraction is larger than 50%, since a clumpiness of C_HII ~ 2−3 is expected, according to recent theoretical predictions (Pawlik et al. 2010).

In the near future, two large surveys will provide further information and better statistics to the z ~ 7 related studies: the ground based wide NIR survey ULTRAVISTA² (Ultra Deep Survey with the VISTA telescope) on the COSMOS field and the deep space based survey CANDELS³ (Cosmic Assembly Near Infra-red Deep Extragalactic Legacy Survey). In particular, the CANDELS survey will have a deep impact on the z = 7 study since it will be able to clean out the present high-z sample from variable objects, to reduce significantly the cosmic variance, to sample in detail the break of the LF at M₁₅₀₀ ~ −20, and to provide a detailed distribution of the half light radii of the candidates, strongly needed to derive a reliable incompleteness function from the simulations.

In the meanwhile, spectroscopic confirmations of z = 7 (and beyond) candidates with optical/NIR instrumentations from the ground (Fontana et al. 2010; Vanzella et al. 2011), with JWST from space or through sub-mm emission lines with ALMA will open new frontiers in the study of the first galaxies.

Acknowledgments

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Institute. STScI is operated by the association of Universities for Research in Astronomy, Inc. under the NASA contract NAS 5-26555. Observations were also carried out using the Very Large Telescope at the ESO Paranal Observatory under Programme IDs LP181.A-0717, LP168.A-0485, ID 170.A-0788, ID 181.A-0485, ID 283.A-5052 and the ESO Science Archive under Programme IDs 67.A-0249, 71.A-0584, 73.A-0564, 68.A-0563, 69.A-0539, 70.A-0048, 64.O-0643, 66.A-0572, 68.A-0544, 164.O-0561, 163.N-0210, and 60.A-9120. We acknowledge partial financial support from ASI.

References

All Tables

All Figures

Fig. 1 The Z − Y vs. Y − J color − color diagram for the HUDF dataset limited at J < 29.2. Big points are observed galaxies, small black points are CB07 models at 0.0 < E(B − V) < 0.3, with small red points those at z > 6.5. The small blue points are the result of our simulations, with synthetic objects inserted in the FITS images and recovered using the same procedure used for the observed catalog. The solid line is the color cut of the A criterion, while the dashed line shows the B criterion (described in Sect. 3). Candidate galaxies at z ~ 7 selected with criterion C are shown in magenta. Vertical arrows are upper limits in the Z band, while the blue circle at Z − Y > 3.6 is the suspected SNa outlined by Oesch et al. (2010a).

In the text

Fig. 2 The different morphologies used for the simulations of WFC3 fields. From top left, the first four cutouts represent non-saturated stars in the HUDF survey (morphological template 1). The SNa candidate outlined by Oesch et al. (2010a) (template 2) is then shown, followed by the three brightest z ~ 7 candidates of the HUDF survey (template 3). The most extended galaxy is the brightest candidate in the HUDF (HUDF-1095 or G2_1408 of Castellano et al. 2010a, template 4). The first thumbnail on the fourth raw shows the stack of the faint candidates selected in the HUDF (template 5). All these thumbnails have been extracted from the J125 band of HUDF. The remaining objects are spectroscopically confirmed galaxies z ~ 4, extracted from the ACS-I band of the HUDF and shifted at z = 7 as described in the main text (template 6), and two artificial exponential disk galaxies with half light radius =0.3 arcsec (templates 7 and 8) and with half light radius =0.1 arcsec (templates 9 and 10).

In the text

Fig. 3 The black histogram shows the observed distribution of the half light radii (in arcsec) of the z ~ 7 galaxy candidates brighter than J = 27.8 in the HUDF. At this magnitude limit the completeness of the J band image is higher than 80% even for the most extended objects. The galaxy with half light radius  ~0.25 arcsec is the HUDF-1095 one (template 4). The cyan solid line shows the uniform distribution in half light radius (template 7), while the dark pink and green solid lines indicate the two log-normal distributions in half light radius adopted by Ferguson et al. (2004) (template 8, M1) and by Oesch et al. (2010b, template 9, M2), respectively. The purple line indicates the distribution derived from Vanzella et al. (2009) (template 10, M3). All sizes are corrected for PSF. The dashed line shows the completeness of galaxy detection for different half-light radii derived by the performed simulations. This curve has been multiplied by a factor of 2 only for clarity.

In the text

Fig. 4 The fraction of detected galaxies versus the J125 input magnitude for the simulations in the HUDF field. Different colors indicate the various morphologies adopted: the dashed magenta lines indicate non-saturated stars on the HUDF survey (morphological template 1). The solid magenta line shows the completeness obtained using the SNa candidate outlined by Oesch et al. (2010a) (template 2) as morphological template, while the red solid curves show the three brightest z ~ 7 candidates of the HUDF survey (template 3). The blue solid line is the completeness derived by adopting as template the brightest candidate in the HUDF (HUDF-1095 or G2_1408 of Castellano et al. 2010a, template 4), and the black solid line represents the completeness obtained adopting as input template the stack of the faint candidates selected in the HUDF (template 5). The green dashed lines indicate the spectroscopically confirmed galaxies z ~ 4, shifted at z = 7 (template 6). The cyan solid line shows the completeness obtained assuming a uniform distribution in half light radius (template 7), while the dark pink and green solid lines indicates the results of two log-normal distributions in half light radius adopted by Ferguson et al. (2004), Oesch et al. (2007, 2010a) (template 8, M1) and by Oesch et al. (2010b) (template 9, M2), respectively. The purple line indicates the completeness derived from Vanzella et al. (2009) (template 10, M3).

In the text

Fig. 5 The LF at z = 7, obtained using as morphological templates real bright (J < 27.4) candidates at z ~ 7 found in the HUDF (template 3). The magenta circle represents the LF found in the Hawk-I from C10a, while the other symbols indicate the LF in the HUDF and ERS derived adopting different color criteria: our color criteria (red filled circles, criterion C), Bunker et al. (2010) (small filled black triangles, criterion B) and Oesch et al. (2010a) (blue squares, criterion A). The blue points are shifted by −0.1 mag and the black points by +0.1 for clarity. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) and Bouwens et al. (2007), respectively.

In the text

Fig. 6 The LF at z = 7, obtained using as morphological types galaxies drawn from the templates 1−6. The color code is the same of Fig. 4. The blue and green bins are shifted by −0.2 and −0.1 mag, while the red and black points are moved by +0.1 and +0.2, respectively, for clarity. The blue dotted lines indicates the range of variability for the best fit of α spanning the range −2.15 ≤ α ≤ −1.38. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) and Bouwens et al. (2007), respectively. The filled black squares are the LF presented in Oesch et al. (2010a).

In the text

Fig. 7 The LF at z = 7, obtained using different distributions for the simulated half light radii, namely the uniform one (template 7, cyan points), and the three log-normal distributions (template 8 using dark pink points, template 9 using dark green points, and template 10 using purple points). The dark green bins are shifted by +0.1 mag, while the dark pink points are moved by −0.1 mag, respectively, for clarity. The blue dotted lines indicate the range of variability for the best fit of α spanning the range −2.15 ≤ α ≤ −1.38. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) and Bouwens et al. (2007), respectively. The filled black squares are the LF presented in Oesch et al. (2010a).

In the text

	Fig. 8 The dependence of the steepness of the LF, α, against the half light radii of the simulated galaxies. Different colors indicate various morphology adopted, using the same color code of Fig. 4.
In the text

Fig. 9 The LF at z = 7, obtained using as morphological templates the stacking of the faint candidates at z ~ 7 found in the J band image of the HUDF (template 5) and the C color criterion. The LF is computed using ERS and HUDF fields together (red filled squares, G10). The green filled squares show the LF obtained with the same color criterion C but using the HUDF-1095 candidate as input for the simulations (template 4). These data have been complemented with the results of Ouchi et al. (2009) (O09), Bouwens et al. (2010b) (B10) and C10b. The lower dashed magenta line is the best fit to all the observed points related to the template 5, the upper dashed magenta line is the best fit related to the template 4, while the black solid line shows the Hawk-I result of C10a. The blue dashed and green dot-dashed lines show the LF at z = 6 of McLure et al. (2009) (M09) and Bouwens et al. (2007) (B07), respectively.

In the text

	Fig. 10 The confidence regions at 68% c.l. (1σ) for the LF parameters and the SFRD. The red points show the projection of the best fit. The plot on the left shows the confidence regions based on the simulations with compact morphology (template 5), while the right plot is related to extended morphology (template 4).
In the text

Fig. 11 The constraints on the maximum CHII and minimum fesc parameters in order to have the Universe reionized by LBGs at z = 7 given our LF. Given a value for fesc, the Universe is reionized if CHII is lower than the red line, if a compact morphology is used (template 5), or lower than the blue line if an extended one is considered (template 4). The ρUV is computed integrating the LF down to M1500 = −10 mag. The solid line at CHII = 1 indicates a homogeneous IGM, and we must consider only the CHII ≥ 1 region, since galaxies at z ~ 7 are formed in biased density regions of the Universe. The dashed line shows the value CHII = 3. In practice, only a combination of (fesc, CHII) parameters inside the shaded regions corresponds to a z ~ 7 IGM kept ionized by LBGs only.

In the text