2.1. Spectral break selection

Succesful studies of the z > 3 universe have thus turned to broadband photometry. Although the spectroscopic resolution of broadband filters is extremely low (R ~ 5 for the SDSS filter set), photometry can still be used to discern strong spectral features. The intrinsic spectra of star-forming galaxies exhibit two relatively strong spectral breaks. The first is the Lyman break at 912 Å, which is the result of the hydrogen ionisation edge in massive stars, combined with the photoelectric absorption of more energetic photons by neutral gas (H i) in the ISM of galaxies. The second break is due to a combination of absorption by the higher order Balmer series lines down to the Balmer limit at 3 646 Å (strongest in A-type stars), along with absorption from metal lines in lower mass stars, primarily the Ca H and K lines (3 934 and 3 969 Å), strongest in lower-mass, G-type stars. Although this so-called ‘4 000 Å break’ can become strong in galaxies dominated by older stellar populations, it is typically much weaker than the Lyman break, which can span an order of magnitude or more in luminosity density.

The IGM adds to the amplitude of the Lyman break, as neutral gas along the line-of-sight (either in the cosmic web, or in the circum-galactic medium of galaxies) efficiently absorbs any escaping ionising radiation (with a rest-frame wavelength less than 912 Å). Additionally, the continuum of galaxy spectra between 912 and 1 216 Å will be attenuated by Lyα absorption lines in discrete systems along the line-of-sight, known as the Lyman-α forest. This latter effect is redshift-dependent, in that at higher redshift, more opacity is encountered along the line-of-sight. By z ≳ 5, the region between the Lyman continuum edge and Lyα is essentially opaque, such that no flux is received below 1 216 Å, compared to 912 Å at lower redshifts. As such, the Lyman break is occasionally referred to as the Lyman-α break at higher redshifts, although the mechanism is generically similar.

At z > 3, the Lyman break feature shifts into the optical, and can thus be accessed from large-aperture, wide-field ground-based telescopes. Building on the efforts of Tyson (Reference Tyson1988), Guhathakurta, Tyson, & Majewski (Reference Guhathakurta, Tyson and Majewski1990), and Lilly, Cowie, & Gardner (Reference Lilly, Cowie and Gardner1991), Steidel & Hamilton (Reference Steidel and Hamilton1993) was amongst the first to realise this tremendous opportunity to study the distant universe. They used a set of three filters (U_n, G, and $\mathscr{R}$ ) and devised a set of colour criteria to select galaxies at z ~ 3 in two fields around high-redshift quasars. At this redshift, the Lyman break occurs between the U_n and G bands, thus a red U_n − G colour (or, a high G/U_n flux ratio) corresponds to a strong break in the spectral energy distribution (SED) of a given galaxy between these filters. While this alone can efficiently find galaxies with strong Lyman breaks, corresponding to redshifts of z ~ 3, it may also select lower redshift galaxies with red rest-frame optical continua. Thus, a second colour is used, $G-\mathscr{R}$ , corresponding to the rest-frame UV redward of the Lyman break at z ~ 3, with a requirement that this colour is relatively blue, to exclude lower redshift passive and/or dusty galaxies from contaminating the sample. I refer the reader to the review by Giavalisco (Reference Giavalisco2002) for further details on Lyman break galaxies at z < 6.

The Lyman break is the primary spectral feature used in nearly all modern searches for z > 6 galaxies. In Section 3 below, I will discuss recent searches with this technique. Some studies continue to use a set of colour criteria, typically involving three filters, qualitatively similar to Steidel & Hamilton (Reference Steidel and Hamilton1993). However, a more recent technique is beginning to become commonplace, known as photometric redshift fitting. With this technique, one compares the colours of a photometric sample to a set of template SEDs in all available filters. The advantage here is that one uses all available information to discern a redshift. Additionally, most photometric redshift tools calculate the redshift probability distribution function, therefore giving a higher precision on the most likely redshift, typically with Δz ~ ±0.2–0.3 at z > 6, compared to Δz ~ ±0.5 for the two-colour selection technique. The disadvantage is that this technique is dependent on the set of SED templates assumed. However, at z ~ 6, this worry is alleviated as the primary spectral feature dominating the photometric redshift calculation is the Lyman break, thus the effect of different template choices is likely minimised. In Figure 1, I show model spectra at three different redshifts, highlighting how the wavelength of the Lyman break, and its strength, changes with redshift, as well as example real galaxies at each of the redshifts shown. The right panel of Figure 1 also shows a colour–colour plot, highlighting how one could use the filter set from the left panel to perform a Lyman break selection. One note for the reader on terminology—frequently authors use the term ‘Lyman break galaxy’ to denote a galaxy selected via the Lyman break technique. However, all distant galaxies, whether they are bright enough to be detected in a continuum survey or not, will exhibit a Lyman break. Thus, in this review, I will use the term ‘continuum selected star-forming galaxies’ to denote such objects, which are selected via either Lyman break colour–colour selection, or photometric redshift selection.

Figure 1. (Top-Left) Model galaxy spectra at three different redshifts, compared to the Hubble Space Telescope optical (ACS) and near-infrared (WFC3) filter set available in the GOODS/CANDELS fields. The models shown have log (M/M_⊙) = 9, an age of 10⁸ yr, and E(B − V) = 0.03. At z > 7, the Lyman break shifts into the near-infrared, rendering such distant galaxies literally ‘invisible’. (Top-Right) Colour–colour plot showing how the colours of normal star-forming galaxies (SFGs) at 4 < z < 10 change with redshift. For the SFGs, the vertical axis represents the Lyman break colour: B − V, i′ − z′, Y − J, and J − H at z = 4, 6, 8, and 10, respectively. The horizontal axis represents the rest-frame far-ultraviolet colour: V − i′, z′ − Y, J − H, and H − [3.6] at z = 4, 6, 8, and 10, respectively. The dark portion of these curves represent when the curve is within Δz ± 0.5 of the centre of the redshift bin. The dashed and dotted lines show the colours of dusty SFGs and passive galaxies from 0 < z < 5, where for these we plot z′ − Y versus i′ − z′ (i.e., showing how these galaxies would compare to z = 6 SFG colours). One can construct a box in each colour–colour combination which selects the desired high-redshift population, and excludes the low-redshift interlopers. Bottom) 3 arcsec stamp images in the seven filters shown in the top-left panel centred on example galaxies at z ~ 4, 6, 8, and 10 (the z ~ 4, 6, and 8 galaxies are spectroscopically confirmed, and come from the sample of Finkelstein et al. Reference Finkelstein2015c, while the z ~ 10 candidate galaxy comes from Bouwens et al. Reference Bouwens2015c).

2.2. Emission line selection

Another method to select distant galaxies is via strong emission lines. At z > 6, the only strong emission line accessible with current technology is the Lyman-α line, at $\lambda _{\text{rest}} =$ 1 216 Å. While a blind spectroscopic survey for this feature is likely not practical at such high redshifts, this line is strong enough that it can be discerned with imaging in narrowband filters. A galaxy with a line at a particular wavelength covered by a narrowband filter will appear brighter in that narrowband then in a broadband filter covering similar wavelengths (as it will have a greater bandpass-averaged flux in the narrowband). This particular line was noted decades ago as a possible signpost for primordial star formation in the early universe (Partridge & Peebles Reference Partridge and Peebles1967), thus a variety of studies were commissioned with the goal of selecting large samples of Lyman-α emitting galaxies (LAEs). Although one of the first narrowband-selected LAEs was discovered by Djorgovski et al. (Reference Djorgovski, Spinrad, McCarthy and Strauss1985), it wasn’t until the advent of large aperture telescopes and/or wide-field optical imagers that the first large samples of LAEs were discovered (e.g., Cowie & Hu Reference Cowie and Hu1998; Rhoads et al. Reference Rhoads2000).

LAEs form a complementary population of galaxies to continuum-selected star-forming galaxies. The study of Steidel & Hamilton (Reference Steidel and Hamilton1993, and subsequent studies from that group) typically restricted their analyses to galaxies with observed optical AB magnitudes brighter than approximately 25, due to the limited depth available from ground-based broadband imaging. As narrowband selection techniques require only evidence of a large flux ratio between the narrowband and encompassing broadband filters, continuum detections are not required, and high equivalent width (EW) Lyα emission from galaxies with much fainter continuum levels can be detected. Deep broadband imaging of ground-based narrowband-selected LAEs shows that they are indeed on average fainter than continuum-selected galaxies, with continuum magnitudes as faint as ~ 28 (for z = 3–5 LAEs, e.g., Ouchi et al. Reference Ouchi2008; Finkelstein et al. Reference Finkelstein, Cohen, Malhotra and Rhoads2009). Whether LAEs form a completely separate population, or are simply the low-mass extension of the more massive continuum-selected galaxy population, remains an active area of study (e.g., Hashimoto et al. Reference Hashimoto2013; Nakajima et al. Reference Nakajima2013; Song et al. Reference Song2014).

2.3. Infrared selection

The previously discussed selection methods relied either on the detection of stellar continuum emission, or of nebular gas emission (due to photoionisation from predominantly stellar-produced ionising photons). An alternative method is also available, selecting galaxies based on the far-infrared emission from UV radiation re-processed by dust grains in the ISM (e.g., Smail, Ivison, & Blain Reference Smail, Ivison and Blain1997; Hughes et al. Reference Hughes1998; Barger et al. Reference Barger1998). The advent of the Herschel Space Observatory as well as large-dish ground-based telescopes such as the James Clerk Maxwell Telescope (JCMT) with fast survey capabilities have allowed searches for rare, highly luminous dusty star-forming galaxies. The redshift distribution of such dusty star-forming galaxies peaks at z ~ 2–3 (with the exact peak redshift depending on the selection-wavelength, with redder wavelengths selecting on average higher redshift galaxies due to the redshifting of the dust-emission SED peak).

Surprisingly, the redshift distribution extends out to z > 5 (see Casey, Narayanan, & Cooray Reference Casey, Narayanan and Cooray2014, and references therein), implying that vast quantities of dust are produced in the early universe. Obtaining spectroscopic redshifts for such sources is difficult, as Lyα is easily attenuated by dust (though see, e.g., Barger et al. Reference Barger1999; Chapman et al. Reference Chapman, Blain, Smail and Ivison2005; Capak et al. Reference Capak2008). However, the advent of the Atacama Large Millimeter Array (ALMA) now allows spectroscopic confirmation via a number of sub/millimetre lines, such as those on the CO ladder, or the [C ii] 158 μm fine structure line (e.g., Riechers et al. Reference Riechers2013). Although dust emission has been detected from galaxies as distant as z ≈ 7.5 (Watson et al. Reference Watson2015), the vast majority of known dusty star-forming galaxies lie at z < 6; thus, I will not discuss these studies further in this review. However, ALMA, combined with a recent update to the Plateau du Bure interferometer (NOEMA), and a potential future update to the Jansky Very Large Array (NGVLA; Carilli et al. Reference Carilli2015), will further enable z > 6 mm studies, and will soon provide key insight into such distant galaxies.

3.1. Broadband searches for star-forming galaxies at z = 6

Selecting galaxies at z = 6 via the Lyman break technique requires imaging at the extreme red end of the optical, in the z-band at 0.9 μm, as galaxies at such redshifts are not visible at bluer wavelengths. Additionally, the imaging must be deep enough to detect these galaxies, as their increased luminosity distance results in an observed magnitude at z = 6 that is ~ 2 mag fainter than a comparably intrinsically bright galaxy at z = 2. Deep imaging in the z-band is difficult both due to the early lack of red-sensitive detectors, as well as (from the ground) the numerous night sky OH emission lines, which set a high background level. Although there were earlier examples of z ≈ 6 galaxies from, e.g., the Hubble Deep Field with WFPC2 (e.g., Weymann et al. Reference Weymann1998), larger samples of z = 6 star-forming galaxies were not compiled until the installation of the red-sensitive optical Advanced Camera for Surveys (ACS) on HST in 2002. While initial studies using early observations found a handful of z ~ 6 candidates (Bouwens et al. Reference Bouwens2003; Stanway, Bunker, & McMahon Reference Stanway, Bunker and McMahon2003; Yan, Windhorst, & Cohen Reference Yan, Windhorst and Cohen2003, including spectroscopic confirmations, e.g., Bunker et al. Reference Bunker, Stanway, Ellis, McMahon and McCarthy2003), larger samples of more than 100 candidate z = 6 galaxies were soon compiled using both the Great Observatories Origins Deep Survey (GOODS; Giavalisco et al. Reference Giavalisco2004a) and Hubble Ultra Deep Field (HUDF; Beckwith et al. Reference Beckwith2006) datasets (Dickinson et al. Reference Dickinson2004; Bunker et al. Reference Bunker, Stanway, Ellis and McMahon2004; Giavalisco et al. Reference Giavalisco2004b; Bouwens et al. Reference Bouwens, Illingworth, Blakeslee and Franx2006).

Significant progress at z = 6 has also been made from the ground. Although difficult due to the bright night sky emission, ground-based surveys can cover much larger areas, and thus provide complementary results on the bright-end of the galaxy luminosity function, which is difficult to constrain with HST due to the small field-of-view of HST’s cameras. Surveys such as the Subaru/XMM-Newton Survey (SXDS), the UKIRT Infrared Deep Sky Survey (UKIDDS)/Ultra Deep Survey (UDS), the Canada–France Hawaii Telescope Legacy Survey (CFHTLS), and the UltraVISTA Survey have searched areas of the sky from 1–4 deg² for z = 6 galaxies (e.g., Kashikawa et al. Reference Kashikawa2006; McLure et al. Reference McLure, Cirasuolo, Dunlop, Foucaud and Almaini2009; Willott et al. Reference Willott2013; Bowler et al. Reference Bowler2015). As discussed in Section 5, these wide-field surveys are necessary to probe luminosities much brighter than the characteristic luminosity at z = 6.

One of the major conclusions from these early studies was that the galaxy population at z = 6 could maintain an ionised IGM only if faint galaxies dominate the ionising budget, which required that the luminosity function must maintain a steep faint-end slope well below L* (e.g., Bunker et al. Reference Bunker, Stanway, Ellis and McMahon2004; Yan & Windhorst Reference Yan and Windhorst2004). Turning to the evolution of the cosmic star-formation rate (SFR) density, Giavalisco et al. (Reference Giavalisco2004b) found that the evolution was remarkably flat out to z = 6, such that the rate of cosmic star formation was similar at z = 6 as at z = 2. However, in a combined analysis using data from multiple HST surveys, Bouwens et al. (Reference Bouwens, Illingworth, Franx and Ford2007) found that there was a steep drop in the SFR density, by more than 0.5 dex from z = 2 to 6. Some of this discrepancy may be due to the fact that many of these early z = 6 galaxies were only detected in a single band, making robust samples (and their completeness corrections) difficult to construct, as well as hampering the ability to derive a robust dust-correction, which is necessary for an accurate measure of the SFR density.

This issue was alleviated with the installation of the Wide Field Camera 3 (WFC3) on HST in 2009, which contains both an ultraviolet (UV)/optical camera (WFC3/UVIS), and a near-infrared camera (WFC3/IR). Three major surveys were initiated with the infrared camera to probe high-redshift galaxies. The Hubble Ultra Deep Field 2009 survey (HUDF09; PI Illingworth) obtained deep imaging in three near-infrared filters (centred at 1, 1.25, and 1.6 μm) on the HUDF as well as two nearby parallel fields, while the Ultra Deep Field 2012 (UDF12; PI Ellis; see Ellis et al. Reference Ellis2013 and Koekemoer et al. Reference Koekemoer2013) survey increased the depth in these filters in the HUDF, and added a fourth filter at 1.4 μm. At the same time, the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; PIs Faber & Ferguson; see Grogin et al. Reference Grogin2011 and Koekemoer et al. Reference Koekemoer2011) was one of the three HST Multi-cycle Treasury Programs awarded in Cycle 18. CANDELS observed both GOODS fields in the same three WFC3/IR filters as the HUDF09 surveyFootnote ¹ , as well as three additional fields (COSMOS; Extended Groth Strip/EGS, and UDS) in the 1.25 and 1.6 μm filters. CANDELS also obtained optical imaging with ACS in parallel, which was particularly useful in the COSMOS, EGS, and UDS fields, which had less archival ACS imaging than the GOODS fields. Using the combination of these new near-infrared data with the previously available ACS optical data, these data now allow full two-colour selection of z = 6 galaxies (Bouwens et al. Reference Bouwens2015b), as well as accurate photometric redshifts (Finkelstein et al. Reference Finkelstein2015c), for a sample of ~ 700–800 robust candidates for z = 6 galaxies.

3.2. Broadband searches for star-forming galaxies at z ⩾ 7

At z = 7, the Lyman break redshifts into the near-infrared, making deep near-infrared imaging a requirement. Prior to the advent of WFC3, the availability of the necessary imaging was scarce. Efforts were made with NICMOS imaging (e.g., Bouwens et al. Reference Bouwens2010a), but with a survey efficiency > 40 × worse than WFC3 (a combination of depth and field size; Illingworth et al. Reference Illingworth2013), significant progress in understanding early galaxy formation was difficult. The first large, robust samples of z > 7 galaxies thus did not come about until the acquisition of the first year of the HUDF09 dataset. Several papers were published in the first few months after the initial Year 1 HUDF09 imaging (e.g., Oesch et al. Reference Oesch2010; Bouwens et al. Reference Bouwens2010b; Bunker et al. Reference Bunker2010; McLure et al. Reference McLure2010; Finkelstein et al. Reference Finkelstein2010), finding ~ 10–20 robust z = 7 candidate galaxies, as well as 5–10 z = 8 candidate galaxies. Sample sizes of faint galaxies were increased with the completed 2-yr HUDF09 dataset, as well as the added depth from the UDF12 programme (e.g., McLure et al. Reference McLure2013; Schenker et al. Reference Schenker2013). The largest addition in sample size of z = 7–8 galaxies was made possible by the CANDELS programme, which, in the two GOODS fields alone, provided 5.4 and 3.7 × more galaxies than the HUDF alone at z = 7 and z = 8, respectively (Finkelstein et al. Reference Finkelstein2015c). A large-area search for z = 8 galaxies was also made possible by the Brightest of Reionizing Galaxies (BoRG; PI Trenti) programme, which used HST pure parallel imaging to find an additional ~ 40 z = 8 galaxy candidates at random positions in the sky. Recently, even larger samples of plausible z = 7 and 8 galaxies were obtained by Bouwens et al. (Reference Bouwens2015b), who extended the search to all five CANDELS fields (as well as BoRG), making use of ground-based Y-band imaging in those regions without HST Y-band (i.e., the COSMOS, EGS, and UDS CANDELS fields). The HST samples of z = 7 and 8 galaxies now number ~ 300–500 and ~ 100–200, respectively.

Even at such great distances, ground-based searches have provided valuable data, particularly at z = 7, though the bright night sky background still limits these studies to be restricted to relatively bright galaxies. Some of the first robust z = 7 candidates discovered from the ground were published by Ouchi et al. (Reference Ouchi2009), who used ground-based near-infrared imaging over the Subaru Deep Field (SDF) and GOODS-N to select 22 z = 7 candidate galaxies, all brighter than 26th magnitude, including some which are spectroscopically confirmed (Ono et al. Reference Ono2012). Castellano et al. (Reference Castellano2010) also used deep ground-based near-infrared imaging, here from the VLT/HAWK-I instrument, to find 20 candidate z = 7 galaxies brighter than 26.7. Tilvi et al. (Reference Tilvi2013) took a complementary approach, using ground-based medium-band imaging to select three candidate z = 7 galaxies from the zFourGE survey. Although the numbers discovered in this latter study were small, the higher spectral resolution afforded by the medium bands allows much more robust rejection of stellar contaminants, particularly brown dwarfs, which can mimic the broadband colours of z ⩾ 6 galaxies (Figure 2). The most recent, and most constraining, ground-based results come from Bowler et al. (Reference Bowler2012) and Bowler et al. (Reference Bowler2014), who used deep, very wide-area imaging 1.65 deg² from the UltraVISTA COSMOS and the UKIDSS UDS surveys to discover 34 bright z = 7 candidates. The combination of the very large area with the depth allowed Bowler et al. (Reference Bowler2014) to have some overlap in luminosity dynamic range with the HST studies, which allows more robust joint constraints on the luminosity function.

Figure 2. A comparison of the SEDs of star-forming galaxies at high redshift with possible lower redshift contaminants. The blue shaded region shows model spectra of high-redshift star-forming galaxies (SFGs), with z = 6 as the upper bound, and z = 8 as the lower bound (both models have log(M/M_⊙) = 9.7, an age of 300 Myr, and A_V = 0.4). The purple and red lines show a dusty star-forming and a passive galaxy, respectively, both at z = 1.3. Cyan, green, and yellow curves denote M, L, and T dwarf star empirical near-infrared spectra, taken as the weighted mean of M, L, and T dwarf standards from the SpeX Prism Spectral Libraries. The gray-shaded regions denote the wavelengths covered by the HST ACS, WFC3/IR, and Spitzer/IRAC imaging used in space-based searches for z > 6 galaxies, with the lower bound denoting the magnitude depths at these wavelengths in the CANDELS Deep and S-CANDELS surveys. All contaminants shown would likely satisfy a Lyman break criterion for a z > 6 galaxy, as they would not be detected in typical optical imaging. However, colours at redder observed wavelengths can begin to distinguish between true high-redshift galaxies and low-redshift contaminants, though this can be difficult when working with low signal-to-noise data.

Searches at even higher redshifts have been performed, with a number of studies now publishing candidates for galaxies at z = 9 and 10. This is exceedingly difficult with HST alone, as at z ≳ 8.8, galaxies will be detected in only the reddest two WFC3 filters (at 1.4 and 1.6 μm), while at z ≳ 9.3, the Lyman break is already halfway through the 1.4 μm filter, rendering many higher redshift galaxies one-band (1.6 μm) detections. However, initial surveys did not include observations in the 1.4 μm band, thus only one-band detections were possible. These can be problematic, as one-band detections can pick up spurious sources such as noise spikes or oversplit regions of bright galaxies; the possibility of such a spurious source being detected in two independant images at the same location is extremely low (see discussion in Schmidt et al. Reference Schmidt2014 and Finkelstein et al. Reference Finkelstein2015c).

The first z ⩾ 9 candidate galaxy published was a single z ~ 10 object found in the HUDF by Bouwens et al. (Reference Bouwens2011a). The addition of 1.4 μm imaging in the HUDF by the UDF12 programme led to further progress, with a handful of two-band detected z ~ 9 candidate galaxies being discovered (Ellis et al. Reference Ellis2013; McLure et al. Reference McLure2013; Oesch et al. Reference Oesch2013). Interestingly, the initial z ~ 10 galaxy from Bouwens et al. (Reference Bouwens2011a) was not detected in this new 1.4 μm imaging, implying that if it is truly at high redshift, it must be at z ~ 12, although there is slight evidence that it may truly be an emission-line galaxy at z ~ 2 (Brammer et al. Reference Brammer2013). Very high redshift galaxies have also been found via lensing from the CLASH and Hubble Frontier Fields programmes, with candidates as high as z ~ 11, although none have been spectroscopically confirmed (Coe et al. Reference Coe2013; Zheng et al. Reference Zheng2012; Zitrin et al. Reference Zitrin2014; Ishigaki et al. Reference Ishigaki2015; McLeod et al. Reference McLeod2015). More recent work by Oesch et al. (Reference Oesch2014) and Bouwens et al. (Reference Bouwens2015b) have increased the sample sizes of plausible z = 9 and 10 candidate galaxies by probing the full CANDELS area. Although extremely shallow 1.4 μm imaging is available (from pre-imaging for the 3D-HST programme), these studies leverage the deep available Spitzer Space Telescope Infrared Array Camera (IRAC) imaging in these fields. These data cover 3.6 and 4.5 μm, which encompasses rest-frame 0.3–0.4 μm at z = 9 and 10, and thus can potentially provide a second detection filter (though this is limited by the shallower depth and much broader point-spread function of the IRAC imaging). The latest results come from Bouwens et al. (Reference Bouwens2015c), which combine the results from Oesch et al. (Reference Oesch2014) and Bouwens et al. (Reference Bouwens2015b) with new candidates from additional 1 μm imaging over selected galaxy candidates, finding a total sample of ~ 15 and 6 robust candidate galaxies at z = 9 and 10, respectively.

3.3. Narrowband searches for star-forming galaxies at z ⩾ 6

There has also been an intensive effort to discover galaxies on the basis of strong Lyα emission with narrowband imaging surveys at z > 6. These have been primarily ground-based, as the narrow redshift window probed combined with the small-area HST cameras renders space-based narrowband imaging inefficient. The narrowband technique has proven highly efficient at discovering large samples of LAEs at z = 3–6 (e.g., Cowie & Hu Reference Cowie and Hu1998; Rhoads et al. Reference Rhoads2000; Gawiser et al. Reference Gawiser2006; Ouchi et al. Reference Ouchi2008; Finkelstein et al. Reference Finkelstein, Cohen, Malhotra and Rhoads2009); thus, clearly an extension to higher redshift is warranted, though surveys at z > 6 are restricted in redshift to wavelengths clear of night sky emission lines. The most complete survey for z > 6 LAEs comes from Ouchi et al. (Reference Ouchi2010), who used the wide-area SuprimeCam instrument on the Subaru Telescope to discover > 200 LAEs at z = 6.6 over a square degree in the SXDS field. Matthee et al. (Reference Matthee2015) have recently increased the area searched for LAEs at z = 6.6 to 5 deg² over the UDS, SSA22, and COSMOS fields, finding 135 relatively bright LAEs.

Moving to higher redshift has proven difficult, as the quantum efficiency of even red-sensitive CCDs is declining. Nonetheless, Ota et al. (Reference Ota2010) imaged 0.25 deg² of the SXDS with SuprimeCam with a filter centred at 9730 Å, finding three candidate LAEs at z = 7.0. Hibon et al. (Reference Hibon, Malhotra, Rhoads and Willott2011) used the IMACS optical camera on the Magellan telescope to find six candidate z = 6.96 LAEs in the COSMOS field, while Hibon et al. (Reference Hibon, Kashikawa, Willott, Iye and Shibuya2012) found eight candidate LAEs at z = 7.02 in the SXDS with a 9 755 Å narrowband filter on SuprimeCam. To observe LAEs at z > 7 requires moving to the near-infrared, which has only recently been possible due to the advent of wide-format near-infrared cameras, such as NEWFIRM on the Kitt Peak 4-m Mayall telescope. An additional complication is the increasing presence of night sky emission lines, which leaves few open wavelength windows, and drives many to use even narrower filters to mitigate the night sky emission as much as possible. One such window is at 1.06 μm, which corresponds to Lyα redshifted to z = 7.7. At z = 7.7, Hibon et al. (Reference Hibon2010) used WIRCam on the CFHT to discover seven candidate LAEs, Krug et al. (Reference Krug2012) used NEWFIRM on the Kitt Peak 4 m to discover four candidate LAEs, and Tilvi et al. (Reference Tilvi2010), also using NEWFIRM, found four additional candidate LAEs. However, the majority of these candidate LAEs remain undetected in accompanying broadband imaging, due primarily to the difficulty of obtaining deep broadband imaging in the near-infrared from the ground, and most also remain spectroscopically unconfirmed (e.g., Faisst et al. Reference Faisst, Capak, Carollo, Scarlata and Scoville2014), though see Rhoads et al. (Reference Rhoads, Hibon, Malhotra, Cooper and Weiner2012) for one exception. Thus, the validity of the bulk of these sources is in question, and requires either deep (likely space-based) broadband imaging, or spectroscopic follow-up. Although the depths of these z ≳ 7 studies vary, a relatively common conclusion is that the LAE luminosity function is likely evolving very strongly at z > 7 compared to that at lower redshifts. We will discuss the physical implications of this perceived lack of strong Lyα emission in Section 7.

3.4. Searches for non-star-forming galaxies at high redshift

The previous sub-sections focussed on searches for star-forming galaxies at high redshift, via either rest-frame UV continuum emission from massive stars, or Lyα emission from H ii regions surrounding such stars. In the local universe, such a search technique would be extremely biased, as it would miss passively evolving galaxies. An ongoing debate is whether such a bias exists at very high redshift. It is conceivable that so close to the Big Bang, galaxies have not had time to quench and stop forming stars, and thus current surveys are highly complete. However, observational evidence for this is lacking, as the detection of passive galaxies with only optical and near-infrared imaging is difficult. Although no robust passive galaxies have yet been discovered at z > 6 (Mobasher et al. Reference Mobasher2005, but see Chary et al. Reference Chary2007), the first robust samples of handfuls of passive galaxies at z > 3 have only recently been compiled with state-of-the-art near-infrared imaging surveys, relying either on photometric selection via the Balmer break, or full photometric redshift analyses (e.g., Muzzin et al. Reference Muzzin2013; Spitler et al. Reference Spitler2014; Nayyeri et al. Reference Nayyeri2014; Stefanon et al. Reference Stefanon, Marchesini, Rudnick, Brammer and Whitaker2013). If such galaxies exist at z > 6, their discovery should be possible with very deep infrared imaging with JWST, allowing selection based on the rest-frame optical emission from lower mass stars.

A large population of such galaxies at z > 6 is not likely, as they would exist at a time < 1 Gyr removed from the Big Bang. For example, a z = 6 galaxy which formed log (M/M_⊙) = 10 in a single burst at z = 20 would have a magnitude of 29 and 26 at 1.6 and 3.6 μm, respectively. Such a galaxy would be detectable in the HUDF presently. The lack of such galaxies places an upper limit on the abundance, although one needs to be cautious as these types of objects may not be selected by some selection techniques, and it is possible that they are presently mis-identified as foreground interlopers.

3.5. Contamination

All of the studies discussed above select galaxy candidates, meaning that their derived SEDs are consistent with them lying at a high redshift, but the vast majority have not had their precise redshifts measured with spectroscopy. I will discuss spectroscopic efforts in the following section, but here I discuss the possible sources of contamination. In Figure 2, I show the SEDs of true high-redshift galaxies, along with the plausible contaminating sources discussed below.

For continuum-selected galaxies, the most common contaminants are lower redshift dusty galaxies, lower redshift passively evolving galaxies, and stars. Low-redshift dusty galaxies can contaminate as they would be observed to have very red colours near the anticipated Lyman break of a true high-redshift galaxy. Similarly, a lower redshift passively evolving galaxy can contaminate if the 4 000 Å break is mistaken for the Lyman break (at z = 6 and 8 the redshifts of such contaminants would be z ~ 1.1 and 1.7, respectively). Both types of contamination can happen, as depending on the depth of the imaging used, these contaminants may not be detected in the bluer of the filters used to constrain the Lyman break, and detected in the redder of the filters. However, both types of galaxies should be rejected as they will also have red colours in the filter combination just redward of the desired Lyman break, while true high-redshift galaxies are likely bluer. Additionally, extremely dusty galaxies may be detected in mid/far-infrared imaging, which is typically much too shallow to detect true high-redshift galaxies. Photometric redshift analysis techniques typically show that the redder the galaxy, the more probability density shifts from a high-redshift solution to a low-redshift solution, reflecting the decreased likelihood that the object in question is a truly red high-redshift galaxy. For galaxies that are very blue, it is trivial to rule out any possibility of either a dusty or passive low-redshift interloper, but there is usually a non-zero chance of such contamination amongst the redder galaxies in a given sample. Contamination estimates from such objects are typically low at < 10% (e.g., Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c), though the difficulty of spectroscopically identifying such interlopers makes it difficult to empirically measure this contamination rate.

Stellar contamination is typically handled differently, as many studies search for and remove stellar contaminants after the construction of the initial galaxy sample (e.g., McLure et al. Reference McLure2006; Bowler et al. Reference Bowler2012, Reference Bowler2014; Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c). At z > 6, the colours of M, L, and T (brown) dwarf stars can match the colours of candidate galaxies due to the cool surface temperatures of these objects. With HST imaging, it can be straightforward to remove the brighter stellar contaminants as the brighter candidate galaxies are all resolved, while stars remain point sources. However, this works less well for fainter galaxies, as near the detection limit, it can be difficult to robustly tell whether a given object is resolved. This is not a major problem for HST studies, as at J > 27, the expected surface density of such contaminating stars in the observed fields is low (Finkelstein et al. Reference Finkelstein2015c; Ryan & Reid Reference Ryan and Reid2016).

The more major concern is at intermediate magnitudes, J = 25–26, where the numbers of candidates are small, yet it can be difficult to robustly discern if an object is spatially resolved. To alleviate this issue, for any objects which may be unresolved one can examine whether its observed colours are consistent with any potential contaminating stellar sources. For this to be possible, one needs to ensure that the photometric bands available can robustly delineate between stellar sources and true high-redshift galaxies; as discussed in Finkelstein et al. (Reference Finkelstein2015c), this requires imaging at 1 μm when working at z = 6–8 (see also Tilvi et al. Reference Tilvi2013, for a discussion of the utility of medium bands). Using a combination of object colours and spatial extent, it is likely that space-based studies are relatively free of stellar contamination. This may be more of a problem with ground-based studies, though with excellent seeing even bright z > 6 galaxies can be resolved from the ground (e.g., Bowler et al. Reference Bowler2014). Future surveys must be cognizant of the possibility of stellar contamination, and choose their filter set wisely to enable rejection of such contaminants.

4.1. Spectroscopy at z = 6–6.5

With the advent of ACS on HST, the frontier for spectroscopic confirmations of galaxies moved to z ~ 6. At these redshifts, Lyα is still accessible with optical spectrographs and was thus an attractive choice for spectroscopic confirmation. However, the extreme distances means that this line will be extremely faint, thus 8–10 m class telescopes were needed to follow them up spectroscopically. One of the first studies to spectroscopically observe Lyα from continuum-selected star-forming galaxies at z ≈ 6 was that of Dickinson et al. (Reference Dickinson2004), who used serendipitous SNe follow-up ACS grism spectroscopy to detect the Lyman continuum break from one galaxy, following it up with LRIS on Keck to discover Lyα emission at z = 5.8. At that same time, Stanway et al. (Reference Stanway2004) used the GMOS optical spectrograph on the Gemini 8.2-m telescope to measure the redshifts to three galaxies discovered in the ACS imaging of the HUDF, at z = 5.8–5.9 [one of these, originally published by Bunker et al. (Reference Bunker, Stanway, Ellis, McMahon and McCarthy2003), is bright enough to spectroscopically detect the Lyman break]. Stanway et al. (Reference Stanway2007) continued this survey, confirming the redshifts to two additional galaxies, at z = 5.9–6.1. Dow-Hygelund et al. (Reference Dow-Hygelund2007) added another six redshifts via Lyα at z = 5.5–6.1. In a series of papers by Vanzella et al., a larger sample of confirmed z ≈ 6 galaxies was obtained with the FORS2 spectrograph on the 8.2-m VLT, culminating in the spectroscopic confirmation of a total of 32 z ~ 6 galaxies (Vanzella et al. Reference Vanzella2006, Reference Vanzella2008, Reference Vanzella2009).

Another effort for z ~ 6 spectroscopic follow-up comes from Stark et al. (Reference Stark, Ellis, Chiu, Ouchi and Bunker2010, Reference Stark, Ellis and Ouchi2011), who used the DEIMOS optical spectrograph on Keck to spectroscopically observe continuum selected star-forming galaxies at 4 < z < 6. In particular, Stark et al. (Reference Stark, Ellis and Ouchi2011) obtained a very deep 12.5 h single mask observation with DEIMOS, measuring the redshifts for 11 galaxies at 5.7 < z < 6.0 via Lyα emission. Stark et al. (Reference Stark, Ellis and Ouchi2011) examined the fraction of galaxies with strong (here defined as EW > 25 Å) Lyα emission, finding that for fainter galaxies (M $_{\text{UV}} > -$ 20.25) it rises from ~ 35% at z = 4, to ~ 55% at z ~ 6 (for brighter galaxies, the fraction rises from ~ 10% at z = 4 to ~ 20% at z = 6). These results imply that galaxies at higher redshifts have a higher escape fraction of Lyα photons, potentially due to reduced dust attenuation. In addition to measuring the redshifts of many galaxies at z ~ 6, the rising fraction of galaxies with detectable Lyα emission with increasing redshift implied that Lyα should continue to be a very useful tool at z > 6.5.

HST provides an alternative to ground-based spectroscopy, as ACS has grism spectroscopic capabilities. In this mode, one obtains very low-resolution spectra for every object in the camera’s field. The main advantage is in the multiplexing. The primary disadvantage is the contamination from overlapping sources, though this can be mitigated by splitting the observations over multiple roll angles. Due to the very low resolution (the G800L grism on ACS has R ~ 100), only very strong emission lines can be detected. However, the very low sky background affords this mode much greater continuum sensitivity, particularly when searching for galaxies at z > 6, where the night sky emission makes continuum detections from the ground problematic. There have thus been a number of HST surveys seeking to confirm galaxy redshifts via a detection of the Lyman continuum break. This provides somewhat less precision than an emission line detection, but if the sharpness of the break can be measured, one can confirm that the break seen in photometry is indeed the Lyman break (which is sharp, in contrast to the 4 000 Å break which is more extended in wavelength; see Figure 2).

The GRAPES survey (PI Malhotra) obtained deep ACS grism observations over the HUDF. Malhotra et al. (Reference Malhotra2005) presented the spectroscopic confirmation of 22 galaxies at 5.5 < z < 6.7 from this survey, detecting the continuum break from galaxies as faint as z _AB ~ 27.5. The PEARS survey (PI Malhotra) extended these observations to cover eight additional pointings in the GOODS fields, culminating in the spectroscopic detection of a Lyman break at z = 6.6 (with z _AB = 26.1). Ground-based follow-up with the Keck 10-m telescope showed a Lyα emission line at z = 6.57 for this galaxy, confirming its high-redshift nature (Rhoads et al. Reference Rhoads2013). The WFC3/IR camera also has grism capability, and there are have been efforts (though none successful at this time; e.g., Pirzkal et al. Reference Pirzkal2015) to confirm redshifts at z > 7 with HST. The in-progress FIGS (PI Malhotra), CLEAR (PI Papovich), and GLASS (PI Treu; Schmidt et al. Reference Schmidt2016) surveys may change this, as the very deep spectroscopy should detect Lyα emission or possibly continuum breaks for galaxies out to z ~ 8.

4.2. Pushing to z > 6.5: A deficit of Lyα emission?

As the first z ~ 7 galaxy samples began to be compiled with the initial WFC3/IR surveys and ground-based surveys, confirmation via Lyα was an obvious next step. However, this proved more difficult than thought. One of the first hints that not all was as expected came from Fontana et al. (Reference Fontana2010), who observed seven candidate z > 6.5 galaxies with FORS2 on the VLT. Given the expected Lyα EW distribution and the magnitudes of their targeted sample, they expected to detect three Lyα emission lines at ⩾ 10σ significance, yet they found none (Lyα emission from one galaxy was detected at 7σ significance at z = 6.972). Progress was still made, as Pentericci et al. (Reference Pentericci2011) and Vanzella et al. (Reference Vanzella2011) each reported the confirmation of two galaxies via Lyα, at z ~ 6.7 in the former, and z = 7.0–7.1 in the latter. Yet, as discussed in Pentericci et al. (Reference Pentericci2011), the fraction of confirmed galaxies was only 25%, much less than the ≳ 50% predicted by Stark et al. (Reference Stark, Ellis and Ouchi2011).

Observations of galaxies in this epoch were also performed by Ono et al. (Reference Ono2012) and Schenker et al. (Reference Schenker2012), confirming Lyα-based redshifts at z = 7–7.2, yet still finding less galaxies that would be the case if the Lyα EW distribution was unchanged from z = 6. Pentericci et al. (Reference Pentericci2014) published a larger sample, reporting Lyα emission from only 12 of 68 targeted sources at z ≳ 6.5. After accounting for the depth of observations and accurate modelling of night sky emission, Pentericci et al. (Reference Pentericci2014) found the fraction of faint galaxies with Lyα EW > 25 Å to be only ~ 30%. This deficit is unlikely to be due to significant contamination, as Pentericci et al. (Reference Pentericci2011) showed a much higher fraction of detectable Lyα emission from galaxies at z ~ 6 selected in a similar way.

Part of these difficulties may be technological in nature, as at z > 6.5, these observations were working at the extreme red end of optical spectrographs, where the sensitivity begins to be dramatically reduced. Until recently, similarly efficient multi-object spectrographs operating in the near-infrared were not available. This changed with the installation of MOSFIRE on the Keck I telescope in 2012. Finkelstein et al. (Reference Finkelstein2013) used MOSFIRE to observe > 40 galaxies at z = 7–8 from the CANDELS survey in GOODS-N, obtaining very deep 5 h integrations over two configurations (of ~ 20 galaxies each). A single emission line was detected, which was found to be Lyα from a galaxy at z = 7.51, the most distant spectroscopic detection of Lyα at that time. Accounting for incompleteness due to wavelength coverage and spectroscopic depth, Finkelstein et al. (Reference Finkelstein2013) found that they should have detected Lyα from ~ six galaxies, finding that the Lyα ‘deficit’ continues well beyond z = 7. Other observations have been performed with MOSFIRE, yet most have only achieved relatively short exposure times, resulting in primarily non-detections (e.g., Treu et al. Reference Treu, Schmidt, Trenti, Bradley and Stiavelli2013; Schenker et al. Reference Schenker, Ellis, Konidaris and Stark2014). Recently, two new record holders for the most distant spectroscopically confirmed galaxy at the time of this writing have been found, at z = 7.73 by Oesch et al. (Reference Oesch2015b) and at z = 8.68 by Zitrin et al. (Reference Zitrin2015), both detected with MOSFIRE. Interestingly, the four highest redshift galaxies known, at z = 7.5–8.7 (including a recent detection at z = 7.66 by Song et al. Reference Song2016b), all appear to have very low Lyα EWs (of < 30 Å, respectively) and thus are not similar to the much higher EW sources frequently seen at z ≲ 6. In Table 1, I summarize the currently known spectroscopically confirmed galaxies at z > 7. Of note here are again the relatively low EWs, especially at z > 7.2, as well as the bright UV magnitudes of the confirmed sources.

Table 1. Spectroscopically confirmed galaxies at z > 7.

The upper portion of the table contains published redshifts based on significantly detected (> 5σ) Lyα emission at z > 7. We include published uncertainties on the equivalent width when available. Not listed are two additional sources which fall in the 4–5σ significance range, from Schenker et al. (Reference Schenker, Ellis, Konidaris and Stark2014) at z = 7.62, and Roberts-Borsani et al. (Reference Roberts-Borsani2016) at z = 7.47. The bottom portion contains the highest redshift spectroscopically confirmed quasar (via several emission lines, including Lyα), and gamma-ray burst (via spectroscopic observations of the Lyman break), respectively.

While there have been some notable successes in the search for Lyα emission at z > 7, in general all studies report Lyα detections from fewer objects than expected, as well as weak Lyα emission from any detected objects. It thus appears that something, either in the physical conditions within the galaxies, or in the universe around them, is causing either less Lyα photons to be produced, or preventing most of them from making their way to our telescopes. I will discuss physical possibilities for this apparent lack of strong Lyα emission in Section 7.1.3.

4.3. Alternatives to Lyα

Given the apparent difficulties with detecting Lyα at z > 6.5, it is prudent to examine whether other emission lines may be useful as spectroscopic tracers. While photometric colours imply these galaxies likely have strong rest-frame optical emission (e.g., Finkelstein et al. Reference Finkelstein2013; Smit et al. Reference Smit2014; Oesch et al. Reference Oesch2015b), spectroscopic observations of for example [O iii] requires JWST. However, there may be weaker rest-frame UV emission lines that can be observed. Erb et al. (Reference Erb2010) published a spectrum of a blue, low-mass star-forming galaxy at z ~ 2 (called BX418) which possessed physical characteristics similar to typical galaxies at z > 6. Amongst the interesting features in the spectrum of this object was detectable emission lines of He ii λ1640 and C iii] λλ1907,1909. Stark et al. (Reference Stark2014) obtained deep optical spectroscopy of 17 similarly low-mass z ~ 2 galaxies, finding nearly ubiquitous detections of C iii]. The strength of this emission was on average 10% that of Lyα. However, at z > 6.5, most of the Lyα is being attenuated (or scattered); for example, the strength of Lyα at z = 7.51 observed by Finkelstein et al. (Reference Finkelstein2013) was only ~ 10% of that expected from the stellar population of the galaxy. Thus, exposures deep enough to detect Lyα may also detect C iii] at very high redshift. Stark et al. (Reference Stark2015) searched for C iii] from two galaxies with known Lyα redshifts at z > 6, obtaining tenuous ~ 3σ detections of C iii] at z = 6.029 and z = 7.213. Part of the difficulty at higher redshift is because the C iii] doublet becomes resolved (splitting the line flux over more pixels), yet these possible detections imply this may be a promising line for future study. Some progress may be made with MOSFIRE on Keck, though multi-object spectrographs on the next generation of telescopes, such as the 25-m Giant Magellan Telescope (GMT) or the Thirty Meter Telescope (TMT), will have the capability to probe these alternative UV emission features to very faint levels, probing the redshifts of galaxies out to z > 10.

5.1. The significance of the UV luminosity function

Distribution functions are an immensely useful quantity to measure as they are relatively straightforward to compute in both observations and theory, and thus provide a direct means to compare the two. Distribution functions of galaxy luminosities, stellar masses, and even velocity dispersions have been measured at a variety of redshifts, leading to detailed insight into the physical processes inherent in galaxy evolution. At z > 6, however, we are limited in what we can measure. The rest-frame UV is the wavelength regime which can be observed very deeply from the ground and with HST, thus the galaxy rest-frame UV luminosity function is the best-studied distribution function at such redshifts. While stellar mass functions are also useful (and will be discussed in the next section), it is much more direct to correct the UV luminosity function for incompleteness, as the luminosity is a direct observable, while the stellar mass is a derived quantity. There is a downside to the UV luminosity, in that it is highly susceptible to dust attenuation, thus to compare observations to theory, simulations must add dust attenuation or observations must attempt to correct for this attenuation.

Observations of luminosity functions at lower redshifts have shown that it typically follows a characteristic shape with a power-law slope at the faint-end and an exponential decline at the bright-end, transitioning at a characteristic magnitude or luminosity typically referred to as the ‘knee’ of the luminosity function. Parameterised by Schechter (Reference Schechter1976), the ‘Schechter function’ requires three parameters to describe this shape: the characteristic magnitude or luminosity at the knee (M* or L*), the power-law slope at the faint-end (α), and the characteristic number density (ϕ*) which is a normalisation factor which defines the overal number density of galaxies. Schechter function parameterisations for luminosity and magnitudes are given in Equations (1) and (2), respectively (in units of number per magnitude or luminosity bin, per volume).

(1) $$\begin{equation} \phi (L) = \phi ^{\ast }\left(\frac{L}{L^{\ast }}\right)^{\alpha } \mathrm{exp} \left(-\frac{L}{L^{\ast }}\right), \end{equation}$$

(2) $$\begin{eqnarray} \phi (M) = 0.4\,\mathrm{ln}\,(10)\,\phi ^{\ast }\,10^{-0.4(M-M^{\ast })(\alpha + 1)} e^{-10^{-0.4(M-M^{\ast })}}. \end{eqnarray}$$

A comparison between the shape of the luminosity function and that of the underlying halo mass function can provide insight into the mechanisms driving galaxy evolution. A simple toy model may assume that the shape of the luminosity function is similar to that of the halo mass function, scaling for some constant baryon conversion efficiency. However, as shown in Figure 4, this is not the case. In this figure, I show the z = 7 luminosity function from Finkelstein et al. (Reference Finkelstein2015c), along with the halo mass function at z = 7 from the Bolshoi ΛCDM simulation (Klypin, Trujillo-Gomez, & Primack Reference Klypin, Trujillo-Gomez and Primack2011), measured by Behroozi et al. (Reference Behroozi, Wechsler and Conroy2013b). I place the luminosity function on this figure by converting from luminosity to stellar mass via the liner relation derived by Song et al. (Reference Song2016a), and scaling vertically such that the two distribution functions touch at the knee. Assuming in this case a ratio of halo mass to stellar mass of 80, the number densities match at log (M _halo/M_⊙) ~ 11.5 (approximately the halo mass of a L* galaxy at this redshift; Finkelstein et al. Reference Finkelstein2015b), yet the number of both more and less massive halos is higher than that of galaxies. To phrase this another way, the conversion of gas into stars in galaxies in both more and less massive halos is less efficient than at the knee.

Figure 4. The cumulative halo mass function from the Bolshoi simulations at z = 6 and 7, shown in red. In blue, I show the cumulative luminosity functions from Finkelstein et al. (Reference Finkelstein2015c), using the relation between stellar mass and UV absolute magnitude from Song et al. (Reference Song2016a), and scaling by a stellar mass-to-halo mass ratio such that the z = 7 functions match at the knee. Even after this scaling, there is still a discrepancy, which is commonly attributed to feedback due to supernovae at the faint-end, and AGN feedback at the bright-end (image of the Crab Nebula from Loll et al. Reference Loll, Desch, Scowen and Foy2013).

Such differences have been observed at all redshifts where robust luminosity functions exist, and a number of physical mechanisms have been proposed for this observation. One mechanism that is currently actively debated is that of feedback. Models which invoke feedback, typically due to supernovae (primarily at the faint-end), stellar radiative feedback, and (primarily at the bright-end) accreting supermassive black hole/active galactic nucleus (AGN) (see discussion of these processes in the review of Somerville & Davé Reference Somerville and Davé2015, and references therein) can more successfully match observations than those which do not include such effects, in which case too many stars are frequently formed. This feedback can heat and/or expel gas from galaxies, effectively reducing, or even quenching further star-formation. Such feedback can explain a variety of observations. For example, the mass-metallicity relation observed at lower redshift (e.g., Tremonti et al. Reference Tremonti2004; Erb et al. Reference Erb2006) can be explained by supernova-driven winds preferentially removing metals from lower-mass galaxies, while the increased potentials from higher mass galaxies allows retention of these metals (e.g., Davé, Finlator, & Oppenheimer Reference Davé, Finlator and Oppenheimer2011).

Given that these physical processes affect the shape of the luminosity function, studying the evolution of this shape with redshift can therefore provide information on the evolution of these processes. Observations have shown that the abundance of bright quasars decreases steeply at z > 3 (e.g., Richards et al. Reference Richards2006). Although bright quasars do exist at z > 6 (Fan et al. Reference Fan2006), they are exceedingly rare. Thus, if AGN feedback was the primary factor regulating the bright-end of the luminosity function, one may expect a decrease in the difference between the luminosity function and the halo mass function at high redshift. Likewise, if supernova feedback became less efficient at higher redshift, one would expect the faint-end slope to steepen at high redshift, approaching that of the halo mass function (α ~ −2). In practice, this is more complicated, as other effects are in play, such as luminosity-differential dust attenuation (e.g., Bouwens et al. Reference Bouwens2014), and perhaps a changing star-formation efficiency (e.g., Finkelstein et al. Reference Finkelstein2015b). Nonetheless, the shape of the luminosity function is one of the key probes that we can now measure which can begin to constrain these processes.

The integral of the rest-frame UV luminosity function is also a physically constraining quantity. As the UV luminosity probes recent star-formation (on timescales of ~ 100 Myr; Kennicutt Reference Kennicutt1998), the integral of the UV luminosity function, corrected for dust, measures the cosmic SFR density in units of solar masses per year per unit volume. One can measure this quantity as a function of redshift, and such a figure shows that from the present day this quantity rises steeply into the past (e.g., Lilly et al. Reference Lilly, Le Fevre, Hammer and Crampton1996; Schiminovich et al. Reference Schiminovich2005), peaking at z ~ 2–3 (e.g., Reddy & Steidel Reference Reddy and Steidel2009), and declining at z > 3 (e.g., Madau et al. Reference Madau1996; Steidel et al. Reference Steidel, Adelberger, Giavalisco, Dickinson and Pettini1999; Bouwens et al. Reference Bouwens, Illingworth, Franx and Ford2007). A recent review of this topic can be found in Madau & Dickinson (Reference Madau and Dickinson2014). The extension of the cosmic SFR density to z > 6 can provide detailed constraints on the build-up of galaxies at early times. If it continues in a smooth trend, as observed from z ~ 3 to 6, it implies a smooth build-up of galaxies from very early times. Alternatively, if the SFR density exhibits a steep dropoff at some redshift, it may imply that we have reached the epoch of initial galaxy formation.

Finally, another use of the integral of the rest-frame UV luminosity function is as a constraint on reionisation. Galaxies are the leading candidate for the bulk of the necessary ionising photons for reionisation. By assuming a (stellar-population dependent) conversion between non-ionising and ionising UV luminosity, one can convert the integral of the UV luminosity function (the specific luminosity density, in units of erg s⁻¹ Hz⁻¹ Mpc⁻³) to an ionising emissivity, in units of photons s⁻¹ Mpc⁻³. This can then be compared to models of the needed ionising emissivity to reionise the IGM, to assess the contribution of galaxies to reionisation. As shown over a decade ago, galaxies much fainter than the detection limit of HST are likely needed to complete reionisation (e.g., Bunker et al. Reference Bunker, Stanway, Ellis and McMahon2004; Yan & Windhorst Reference Yan and Windhorst2004). Thus, measuring an accurate faint-end slope is crucial to allow a robust measurement of the total UV luminosity density, and thus the total ionising emissivity. I will cover this issue in Section 7.

5.2. Observations at z = 6–10

In Sections 3.1 and 3.2, I covered recent surveys for star-forming galaxies at z ⩾ 6. Here, I will discuss the measurements of the rest-frame UV luminosity function from these surveys. A number of papers have studied this quantity at z = 6 (Bouwens et al. Reference Bouwens, Illingworth, Franx and Ford2007; McLure et al. Reference McLure, Cirasuolo, Dunlop, Foucaud and Almaini2009; Willott et al. Reference Willott2013; Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c; and Bowler et al. Reference Bowler2015; Atek et al. Reference Atek2015; Livermore, Finkelstein, & Lotz Reference Livermore, Finkelstein and Lotz2016), z = 7 (Ouchi et al. Reference Ouchi2009; Castellano et al. Reference Castellano2010; Bouwens et al. Reference Bouwens2011b; Tilvi et al. Reference Tilvi2013; McLure et al. Reference McLure2013; Schenker et al. Reference Schenker2013; Bowler et al. Reference Bowler2014; Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c; Atek et al. Reference Atek2015; Livermore, Finkelstein, & Lotz Reference Livermore, Finkelstein and Lotz2016), z = 8 (Bouwens et al. Reference Bouwens2011b; McLure et al. Reference McLure2013; Schenker et al. Reference Schenker2013; Schmidt et al. Reference Schmidt2014; Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c; Atek et al. Reference Atek2015; Livermore et al. Reference Livermore, Finkelstein and Lotz2016), z = 9 (McLure et al. Reference McLure2013; Oesch et al. Reference Oesch2013, Reference Oesch2014; McLeod et al. Reference McLeod2015; Bouwens et al. Reference Bouwens2015c; Ishigaki et al. Reference Ishigaki2015), and z = 10 (Oesch et al. Reference Oesch2015a; Bouwens et al. Reference Bouwens2015b, Reference Bouwens2015c). In the interest of presenting constraints from the most recent studies in a concise manner, I will focus on the studies of Bowler et al. (Reference Bowler2014, Reference Bowler2015), Finkelstein et al. (Reference Finkelstein2015c), and Bouwens et al. (Reference Bouwens2015b) at z = 6–8, and Bouwens et al. (Reference Bouwens2015c) and McLeod et al. (Reference McLeod2015) at z = 9–10.

At z = 6–8, Finkelstein et al. (Reference Finkelstein2015c) and Bouwens et al. (Reference Bouwens2015b) used data from the CANDELS and HUDF surveys, while Bowler et al. (Reference Bowler2014, Reference Bowler2015) used ground-based imaging from the UltraVISTA and UKIDSS UDS surveys to discover brighter galaxies. Finkelstein et al. (Reference Finkelstein2015c) used only data from the CANDELS GOODS-N and GOODS-S fields, which have deep HST imaging in seven optical and near-infrared filters, versus only four filters in the other three fields. Specifically, only the CANDELS GOODS fields have deep space-based Y-band imaging, which is necessary for robust removal of stellar contaminants (Finkelstein et al. Reference Finkelstein2015c). Bouwens et al. (Reference Bouwens2015b) used all five CANDELS fields, making use of ground-based optical and Y-band imaging to fill in the missing wavelengths from HST. Both studies use the HUDF and associated parallels, while Finkelstein et al. (Reference Finkelstein2015c) also used the parallels from the first year of the Hubble Frontier Fields programme, and Bouwens et al. (Reference Bouwens2015b) used the BoRG/HIPPIES pure parallel programme data (Schmidt et al. Reference Schmidt2014).

In spite of different data reduction schemes, selection techniques, data used, and completeness simulations, the results of Finkelstein et al. (Reference Finkelstein2015c) and Bouwens et al. (Reference Bouwens2015b) are broadly similar (Figure 6). Both studies find a characteristic magnitude M* which is constant or only weakly evolving from z = 6–8, and both find a significantly evolving faint-end slope (to steeper values at higher redshift), and characteristic number density ϕ* (to lower values at higher redshift). This is a change from initial studies at z > 6 (most probing smaller volumes), which found that M* significantly evolved to fainter values from z = 4 to 8, with less evolution in ϕ* (e.g., Bouwens et al. Reference Bouwens, Illingworth, Franx and Ford2007, Reference Bouwens2011b; McLure et al. Reference McLure2013). The faint-end slope α now appears to match (or even exceed) the value from the halo mass function (α ~ −2) at the highest redshifts (α < −2 can be possible due to effects of baryonic physics). The primary difference between these studies appears to be in the normalisation, as the Bouwens et al. (Reference Bouwens2015b) data points are systematically slightly higher/brighter than those of Finkelstein et al. (Reference Finkelstein2015c), which can be attributed to differences in the assumed cosmology (~ 5%), aperture corrections utilised to calculate the total fluxes in the photometry, and differences in contamination.

As Bowler et al. (Reference Bowler2014, Reference Bowler2015) used ground-based data to probe larger volumes, they were thus sensitive to brighter galaxies than either Finkelstein et al. (Reference Finkelstein2015c) or Bouwens et al. (Reference Bouwens2015b). Broadly speaking their results are consistent with the HST-based studies where there is overlap. However, there does seem to be a modest tension, in that the Bowler et al. (Reference Bowler2014, Reference Bowler2015) ground-based results exhibit slightly lower number densities than those from either of the space-based studies, though the tension only exceeds 1σ significance at the faintest ground-based magnitude (M _UV = −21.5). This is true compared to both HST studies at z = 7, while Bowler et al. (Reference Bowler2015) and Finkelstein et al. (Reference Finkelstein2015c) are in agreement at z = 6. To fit a Schechter function, Bowler et al. (Reference Bowler2015) combine their data with data at fainter luminosities from Bouwens et al. (Reference Bouwens, Illingworth, Franx and Ford2007) at z = 6, while Bowler et al. (Reference Bowler2014) combine with fainter data from McLure et al. (Reference McLure2013) at z = 7. The combination of the deeper HST imaging with their much larger volumes allows the ground-based studies to perhaps place tighter constraints on M*. They do find more significant evolution in M* than that found by Finkelstein et al. (Reference Finkelstein2015c) or Bouwens et al. (Reference Bouwens2015b), with M* _{z = 7} = −20.56 ± 0.17, compared to − 21.03^+0.37 _−0.50 from Finkelstein et al. (Reference Finkelstein2015c) and − 20.87 ± 0.26 from Bouwens et al. (Reference Bouwens2015b). Given these uncertainties, some evolution in M* towards fainter luminosities at higher redshift is certainly plausible (including the modest $\text{d}M/\text{d}z =$ 0.2 proposed by Bowler et al.), and there appears no strong disagreement between these complementary studies. The data points from these three studies, as well as a number of other recent works, are shown in Figure 5, and the fiducial Schechter function parameters are shown in Figure 6.

Figure 5. A compilation of luminosity function data from the literature. Data from space-based surveys are shown in blue, and ground-based surveys in green. In each panel, I show the reference Schechter function fit (Section 5.3) to all available data points as the red curves. The lower right panel overplots the fiducial Schechter functions together at all five redshifts with darker colors denoting higher redshifts. These Schechter function values and associated uncertainties are given in Table 2. The studies used in the fitting are: Bouwens et al. (Reference Bouwens2015b) at z = 4–10; Finkelstein et al. (Reference Finkelstein2015c) at z = 4–8; van der Burg, Hildebrandt, & Erben (Reference van der Burg, Hildebrandt and Erben2010) at z = 4–5; McLure et al. (Reference McLure, Cirasuolo, Dunlop, Foucaud and Almaini2009) at z = 5–6; McLure et al. (Reference McLure2013) at z = 7–9; Schenker et al. (Reference Schenker2013) at z = 7–8; Bouwens et al. (Reference Bouwens2015c) at z = 9–10; Bowler et al. (Reference Bowler2015) at z = 6; Castellano et al. (Reference Castellano2010), Tilvi et al. (Reference Tilvi2013) and Bowler et al. (Reference Bowler2014) at z = 7; Schmidt et al. (Reference Schmidt2014) at z = 8; Oesch et al. (Reference Oesch2013) and McLeod et al. (Reference McLeod2015) at z = 9; and Oesch et al. (Reference Oesch2014) and McLeod, McLure, & Dunlop (Reference McLeod, McLure and Dunlop2016) at z = 10.

Figure 6. The derived evolution of the three Schechter function parameters with redshift, derived by fitting all redshifts simultaneously. The shaded blue regions denote the 68% confidence range of the linear evolution of these parameters with redshift, while the circles denote the reference value at a given redshift. Evolution with increasing redshift in M* (to fainter values), α (to steeper values), and ϕ* (to lower values) is detected at > 10σ significance.

Given the wide dynamic range now probed in luminosity, each of the aforementioned studies understandably pay careful attention to the shape of the luminosity function. Specifically, they investigate whether a Schechter function shape is required by the data, or whether another function, such as a double power law (where the bright-end exponential cut-off is replaced by a second power law), or even a single power law (where the exponential cut-off disappears altogether) is demanded. Finkelstein et al. (Reference Finkelstein2015c) considered all three functions. While they found significant evidence that a Schechter function was required at z = 6 and 7, at z = 8 the luminosity function was equally well fit by a single power law as by a Schechter function. This is intriguing, as this is exactly the signature one might expect were AGN feedback to stop suppressing the bright-end (changing dust attenuation is likely not a dominant factor, as bright galaxies at z = 6–8 have similar levels of attenuation; Finkelstein et al. e.g., Reference Finkelstein2012b; Bouwens et al. e.g., Reference Bouwens2014; Finkelstein et al. e.g., Reference Finkelstein2015b). Bouwens et al. (Reference Bouwens2015b) found a similar result, claiming that there was no overwhelming evidence to support a departure from a Schechter function, but that the available data at z > 6 made it difficult to constrain the exact functional form.

While Bowler et al. (Reference Bowler2015) reach a similar conclusion at z = 6, in that either a double power law or a Schechter functional form could explain the bright-end shape, at z = 7 Bowler et al. (Reference Bowler2014) find significant evidence that a double power law is a better fit than Schechter. With this result, they find that the shape of the bright-end closely matches that of the underlying halo mass function, implying little quenching in bright galaxies at z = 7. Bouwens et al. (Reference Bouwens2015b) discuss this difference in conclusions, and attribute it to the differences in the measured number density at M = −21.5; the higher number densities found by the HST studies allow a Schechter function to be fit equally well. In any case, at z > 6, there no longer remains overwhelming evidence to support a Schechter function parameterisation. To robustly constrain the shape of the luminosity function requires further data, in particular in the overlap range between the ground-based and HST studies.

At higher redshifts of z = 9–10, the data do not currently permit such detailed investigations into the shape of the luminosity function. However, we have begun to gain our first glimpse into the evolution of the Schechter parameters to the highest redshifts (if, in fact, this function ultimately does describe the shape of the luminosity functions). Studies of galaxies at these redshifts are difficult with the data used in the above studies, as at these redshifts the Lyman break passes through the HST/WFC3 F125W filter, rendering most objects only detectable in the F160W filter in the CANDELS fields. Nonetheless, several groups have attempted to select F125W-dropout galaxies, with some using Spitzer/IRAC data as a potential secondary detection band (e.g., Oesch et al. Reference Oesch2014). Additionally, recent surveys have started adding observations in the HST/WFC3 F140W filter, which should show a detection for true z ≈ 9 galaxies, though galaxies at z > 10 may still have most of their flux attenuated by the IGM in that band.

McLeod et al. (Reference McLeod2015) probed the faint-end at z ~ 9 using data from the first year of the Hubble Frontier Fields (though they explicitly do not include regions with high magnifications), which contains F140W imaging, allowing for robust two-band detections. They found 12 new z ~ 9 galaxy candidates in these data, which they combined with previously discovered z ~ 9 candidates in the HUDF from McLure et al. (Reference McLure2013) to constrain the luminosity function. They do not fit all Schechter parameters, choosing to leave the faint-end slope fixed to α = −2.02, and explore possible values of the characteristic magnitude and normalisation under the scenarios of pure density or pure luminosity evolution. They find only modest evolution from z = 8, consistent with a continued smooth decline in the UV luminosity density. Bouwens et al. (Reference Bouwens2015c) recently searched the CANDELS fields for bright z ~ 9 and 10 galaxies, finding nine z ~ 9 galaxies, and five z ~ 10 galaxies (all with H < 27). They do not attempt to derive Schechter function parameters, rather they use the redshift evolution of said parameters from Bouwens et al. (Reference Bouwens2015b) to show that it is consistent with their updated results at the bright-end, as well as results from the literature at the faint-end. They however conclude that the UV luminosity density at z ~ 9 and 10 is ~ 2 × lower than would be expected when extrapolating the observed trend from z ~ 4–8.

These studies highlight the remarkable capability of the modestly-sized Hubble to study galaxies to within 500 Myr after the Big Bang. However, the sizable uncertainties remaining, especially at z > 8, lead to fundamental disagreements about the evolution of the UV luminosity density (Section 5.4) which will likely not be resolved until the first datasets come in from JWST.

5.3. A reference luminosity function

As discussed at the beginning of the previous subsection, many groups have published measurements of rest-frame UV luminosity functions at z ⩾ 6. There are a variety of datasets used, both deep space-based, and relatively shallow ground-based, thus one study alone may not have the dynamic range to fully constrain the full luminosity function. As a perhaps crude, yet illustrative attempt to shed light on the evolution of this useful tool, in this subsection I combine published results from a number of studies in an attempt to derive a set of ‘reference’ rest-frame UV luminosity functions at z = 4–10Footnote ² . I use the data from all studies listed in the first paragraph of Section 5.2 where the luminosity function data was available, with the exception of the z = 6 results from Willott et al. (Reference Willott2013) as Bowler et al. (Reference Bowler2015) found that the Willott et al. sample did not include many true z ~ 6 galaxies, possibly due to the shallowness of the earlier data. Likewise, I do not include the z = 7 results from Ouchi et al. (Reference Ouchi2009), as they applied a very high contamination correction of 50%; see discussion in Appendix F.3 of Bouwens et al. (Reference Bouwens2015b). I also do not include the data from Bouwens et al. (Reference Bouwens, Illingworth, Franx and Ford2007, Reference Bouwens2011b), as they are superseded by Bouwens et al. (Reference Bouwens2015b). Finally, I do not include the recent lensed luminosity functions of Ishigaki et al. (Reference Ishigaki2015), Atek et al. (Reference Atek2015), or Livermore et al. (Reference Livermore, Finkelstein and Lotz2016) due to the likely strong presence of Eddington bias (see discussion in Livermore et al. Reference Livermore, Finkelstein and Lotz2016).

I acknowledge that many of these studies use similar datasets in the same survey fields; thus, the same galaxies may contribute to multiple data points, and this analysis does not include possible systematic effects which could broaden the error budget. However, even studies which utilise the same surveys use a wide range of data reduction, photometry, and sample selection techniques which can and does result in differences in the measured number densities at a given magnitude (Figure 5). The results of this analysis are illustrative of the constraints possible when marginalising over these differences, as well as when combining ground with space-based datasets. However, I do acknowledge that the faintest data points used typically all come from the same field: the HUDF. Therefore, during the MCMC analysis (described below), in magnitudes bins at − 18 or fainter each step of the MCMC chain randomly chooses one data point per magnitude bin. This should suppress any ‘artificial’ damping of the uncertainties on the faint-end slope, as during each step, we use only a single measurement from the HUDF per magnitude bin.

From each of the studies used, I extracted the measured, completeness corrected, number densities (and associated uncertainties) as a function of UV absolute magnitude, shown in Figure 5. I then calculated the likelihood that the data were represented by a given Schechter function via a Markov Chain Monte Carlo (MCMC) algorithm. I used an IDL implementation (R. Ryan, private communication) of the Python ${\tt Emcee}$ package (Foreman-Mackey et al. Reference Foreman-Mackey, Hogg, Lang and Goodman2013), which utilises an affine-invariant ensemble sampler to sample the parameter space. However, rather than fitting a Schechter function to each redshift independently, I made the assumption that the Schechter function parameters linearly vary with redshift. This is an assumption often assumed in the literature, though typically after the fact (e.g., Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c). The advantage of this method is that it allows the data from all redshifts to be fit simultaneously, solving for the posterior distribution of the coefficients of the linear function which describes M*(z), ϕ*(z), and α(z). In effect, this replaces 21 parameters (the three Schechter parameters at seven redshifts) with six parameters: The slope and intercept of the three redshift-dependent Schechter parameter trends. An additional advantage is that this method allows convergence of the chain even at z = 9 and 10, where there are few datapoints. Without the assumption of linear evolution, a prior would have to be placed on the individual Schechter parameters at these highest redshifts to obtain a result.

I ran the MCMC algorithm for 10⁶ total steps, and checked each parameter for convergence using the Geweke diagnostic, which compares the average and the variance of the first 10% to the last half of the steps (Geweke Reference Geweke1992). For each of the six parameters of the fit, I calculated the median value and the 68% confidence range from the median and central 68% range of the posterior distribution, respectively. The results are:

$$\begin{eqnarray*} M^{\ast }(z) &=& -20.79_{-0.04}^{+0.05} + 0.13_{-0.01}^{+0.01}\, (z-6),\\ \alpha (z) &=& \,\!\!-1.91_{-0.03}^{+0.04} - 0.11_{-0.01}^{+0.01}\, (z-6),\\ log\,\phi ^{\ast }(z) &=& \,\!\!-3.37_{-0.04}^{+0.05} - 0.19_{-0.01}^{+0.01}\, (z-6). \end{eqnarray*}$$

The fiducial Schechter functions from this exercise derived by evaluating the above equations at a given redshift are shown alongside the data points in Figure 5 and the values are tabulated with their uncertainties in Table 2. In the lower right-hand panel of Figure 5, I show the fiducial Schechter functions from this analysis at z = 4–10 together, highlighting a relatively smooth decline from z = 4 to 10. In Figure 6, I show the redshift-dependent Schechter parameter evolution with the median values at each redshift, along with the results of a variety of studies listed above. The reference results are for the most part consistent with results in the literature within the uncertainties [with the exception of M* at z = 4, which is brighter than the value found by Bouwens et al. (Reference Bouwens2015b) and Finkelstein et al. (Reference Finkelstein2015c)].

Table 2. Reference Luminosity Function Schechter Parameters.

Median Schechter function values to a compilation of data from the literature. These values are derived from a MCMC fit with a prior that the evolution of these parameters are linear with redshift. The value of each parameter is the median value of that parameter’s linear function evaluated at the given redshift. The uncertainty listed is the 68% confidence range.

In the discussion in Section 5.2, I highlighted recent investigations into the evolution of the faint-end slope α and the characteristic magnitude M*. As shown in the left-hand panel of Figure 6, the results at 4 < z < 7 are consistent with a significant dimming in M* with increasing redshift, with $\text{d}M^{\ast }/\text{d}z =$ 0.13. The sign of this effect is similar to that found in Bowler et al. (Reference Bowler2015), though they found a more accelerated evolution, with $\text{d}M^{\ast }/\text{d}z =$ 0.20–0.25 over 5 < z < 7. Investigating α and ϕ*, the reference luminosity functions confirm the previous results of a steepening faint-end slope and decreasing characteristic number density with increasing redshift, with $\text{d}\alpha /\text{d}z =-$ 0.11 and $\text{d}{\log}\phi ^{\ast }/\text{d}z =-$ 0.19. These trends are qualitatively similar to those from Bouwens et al. (Reference Bouwens2015b) of dα/dz = −0.10 and dlogϕ*/dz = −0.27, and from Finkelstein et al. (Reference Finkelstein2015c) of dα/dz = −0.19 and dlogϕ*/dz = −0.31. However, due to the higher fidelity of the reference Schechter fits, the significance of the trends derived here is higher, such that a > 10σ significance evolution is detected in both α and ϕ*, respectively (though the partial correlation between some of the data points may be responsible for some of the apparent improvement). While I did not include recent lensed-galaxy data from the Hubble Frontier Fields, the faint-end slopes I derive at z = 6–8 are in excellent agreement with those from the Frontier Fields of α ~ −2 (Atek et al. Reference Atek2015; Livermore et al. Reference Livermore, Finkelstein and Lotz2016).

A general conclusion which can be made from this exercise is that the current data can reasonably constrain all three Schechter function parameters out to z = 10 when all redshifts are fit simultaneously, under the assumption that the parameters evolve linearly with redshift. By inspecting Figure 5, one can see that the resultant parameterisations at each redshift appear to excellently describe the data; thus, it would appear that this assumption is valid, at least at the limit of the data presently in hand. By fitting all data simultaneously, I arrive at potentially more precise values of the evolution of the rest-frame UV luminosity function with redshift, which increases the constraints on reionisation which can be derived by integrating the UV luminosity function. I explore this further in Section 7.

5.4. The evolution of the cosmic star-formation rate density

A final quantity one can study with the UV luminosity function is that of the evolution of the cosmic SFR density. This is derived by integrating the UV luminosity function down to a specified threshold magnitude, and then converting to a SFR using a conversion based on stellar population modelling (e.g., Madau & Dickinson Reference Madau and Dickinson2014). As this is derived from the UV which is susceptible to dust reddening, an extinction correction also needs to be estimated and applied. As shown in a variety of previous studies, the evolution of this quantity has been shown to smoothly vary with redshift at 3 < z < 8, such that the data are consistent with a single power-law ∝(1 + z)^{−(3 − 4)} (e.g., Madau & Dickinson Reference Madau and Dickinson2014; Oesch et al. Reference Oesch2014; Finkelstein et al. Reference Finkelstein2015c; Bouwens et al. Reference Bouwens2015b). The evolution to higher redshift has been less well constrained, with some observations supporting a continued smooth decline with increasing redshift at z > 8 (e.g., Ellis et al. Reference Ellis2013; Coe et al. Reference Coe2013; McLeod et al. Reference McLeod2015), while others report that their data require a steeper decline (e.g., Oesch et al. Reference Oesch2014; Bouwens et al. Reference Bouwens2015c, Reference Bouwens2015b).

To explore whether the reference luminosity functions derived here can distinguish between these two scenarios, in Figure 7, I show the values of the SFR density as a function of redshift from this reference analysis, alongside values at z < 4 from the compilation of Madau & Dickinson (Reference Madau and Dickinson2014). These values were derived by integrating the luminosity functions at each redshift to − 17 (a value typically used; Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c), to calculate the observed UV luminosity density. This was then converted into a SFR density using the conversion factor of 1.15 × 10⁻²⁸ M_⊙ yr⁻¹/(erg s⁻¹ Hz⁻¹) from Madau & Dickinson (Reference Madau and Dickinson2014). Dust corrections were derived using the redshift-dependent M _UV – β relation from Bouwens et al. (Reference Bouwens2014, see Section 6.1, including at 0.35 dex scatter) coupled with the relation between dust extinction and β from Meurer, Heckman, & Calzetti (Reference Meurer, Heckman and Calzetti1999). Formally, the combination of these two relations can result in negative values of A _UV; these were set to zero. The total dust correction factors (L _intrinsic/L _observed) were 3.0, 2.2, 1.7, 1.4, and 1.4 at z = 4, 5, 6, 7, and 8, respectively. Zero dust attenuation was assumed at z = 9 and 10.

Figure 7. The evolution of the cosmic star-formation rate density, comparing the values from the integral of the reference luminosity function to those from the literature. All points have been corrected to represent a lower limit on the luminosity function integration of M _UV < − 17, and have been corrected for dust attenuation (with the exception being the low-redshift far-infrared datapoints from Madau & Dickinson Reference Madau and Dickinson2014). The solid blue curve shows a power-law fit to the reference data at 4 < z < 8 (∝[1 + z]^{−4.17 ± 0.27}), extrapolated to higher redshift, while the dashed line shows a fit only to the data at z ⩾ 8 (∝[1 + z]^{−5.10 ± 0.69}). The results from the reference luminosity function are consistent with a smooth decline in the SFR density at 4 < z < 10, with no significant evidence for an accelerated evolution at z > 8. However, the 68% confidence range on the total SFR density (blue-shaded region; derived from integrating the reference luminosity functions to M _UV = − 13) is consistent with an even shallower decline in the SFR density over 4 < z < 10. The light purple region denotes constraints on the total luminosity density from unresolved background fluctuations (Mitchell-Wynne et al. Reference Mitchell-Wynne2015), which also imply a relatively shallow evolution of the total SFR density. The potential for a surprisingly high SFR density at z > 8 will soon be settled by JWST.

I show my results in Figure 7 alongside several results from the most recent studies at each redshift in the literature. Those literature values that were integrated to brighter limits (typically − 17.7) were corrected down to − 17 using the luminosity functions from each paper. Those which did not perform a dust correction had a correction applied, using the value for a given redshift derived here as discussed in the previous paragraph. We note that the reference results are consistent with those from Finkelstein et al. (Reference Finkelstein2015c) and McLeod et al. (Reference McLeod2015), while those from Bouwens et al. (Reference Bouwens2015b) are elevated at z ⩽ 7. The latter discrepancy can be explained as the Bouwens et al. (Reference Bouwens2015b) luminosity functions are slightly elevated over those from other studies.

Madau & Dickinson (Reference Madau and Dickinson2014) found that the low-redshift evolution of the SFR density was proportional to (1 + z)^2.7, while at high redshift, it was proportional to (1 + z)^−2.9 (fitting only to the previously available data at z < 8). To explore the evolution of the SFR density implied by the reference values derived here, I fit a single power law to the reference data at 4 < z < 10, finding that the SFR density is ∝ (1 + z)^{−4.55 ± 0.14}. Fitting to the data only at 4 < z < 8, the slope is slightly shallower, (1 + z)^{−4.17 ± 0.27}, in excellent agreement with the slope of − 4.3 found by Finkelstein et al. (Reference Finkelstein2015c) over that same redshift range. The fit to the 4 < z < 8 data is shown as the solid blue line in Figure 7. We extrapolate this line out to z = 10, and find that it is generally in good agreement with the reference SFR density values at those redshifts, although they do appear to be slightly lower.

To explore this further, I fit a separate power law only to the data at 8 < z < 10, finding a power-law slope of − 5.10 ± 0.69. The difference between this slope and that derived from the data at 4 < z < 8 is only 1.3σ. Therefore, I conclude there is no significant evidence for an accelerated evolution of the SFR density at z > 8 given the currently available data. This is contrary to the evolutionary trend proposed by Oesch et al. (Reference Oesch2014) of (1 + z)^−10.9, shown by the dashed purple line, which becomes progressively more discrepant with the reference shallower evolution with increasing redshift at z > 8.

I note that by forcing the assumption that each of the Schechter function parameters evolve linearly with redshift, I effectively bias the reference SFR density results to evolve smoothly with redshift. However, as shown by the excellent agreement between the data and the Schechter fits in Figure 5, the reference results supporting a smooth evolution of the SFR density are fully consistent with the data out to those high redshifts. To examine this quantitatively, I compared the fiducial results at z = 9 and 10 to those from a single Schechter function fit at each redshift (holding M* fixed at the z = 8 value of − 20.5 to allow the fits to converge). Comparing the goodness-of-fit via the Bayesian Information Criterion (Liddle Reference Liddle2004), I find no evidence that the fiducial fit is worse than when fitting the z = 9 and 10 luminosity functions without the assumption of linear evolution. Therefore, I conclude that this assumption is not biasing these results.

As a final point here, I emphasise that a decline in this quantity—the SFR density for galaxies with M _UV ⩽ −17—is not unexpected, as we are comparing SFR densities above a fixed absolute magnitude. As the luminosity function evolves with increasing redshift, a progressively larger fraction of the total UV luminosity density will come from galaxies with UV absolute magnitudes fainter than − 17 due to the steepening of the faint-end slope. Assuming that the luminosity function extends with the measured faint-end slope down to M _UV = −13, the SFR density derived at M _UV < −17 is ~ 60% of the total SFR density at z = 6, yet only ~ 20% at z = 10. Thus, in the regime of an evolving luminosity function, the choice of the limiting magnitude directly affects the inferred evolution in the SFR density.

To illustrate this point, the shaded blue region in Figure 7 shows the 68% confidence range on the total SFR density from the reference UV luminosity functions (integrated down to M = −13). This is observed to evolve only very shallowly from z ~ 4–8, leaving open the intriguing possibility that there is significant star formation hiding below our current detection limits at z = 10. A resolution of this issue will require a more robust measurement of the total SFR density at z = 10, which we will have in just a few years with the launch of JWST. A JWST deep field should be able to detect galaxies down to M _UV = −16 at z = 10; more than 1.5 mag fainter than the currently available sample of galaxies. However, the recent detection of near-infrared background fluctuations in the CANDELS GOODS-S field (Mitchell-Wynne et al. Reference Mitchell-Wynne2015) imply a total UV luminosity density at z > 8 consistent with our estimate of the total SFR density, thus substantial star-formation activity at z > 8 may yet be found.

6.1. Dust attenuation and chemical enrichment

One of the most direct, and straightforward measures of galaxy evolution is the evolution of galaxy colours. Although the interpretation of these colours requires assumptions and/or modelling, the direct measurement is unambiguous, and can be done with a simple photometric catalogue. At lower redshift, simple colour–colour or colour–magnitude plots are immensely useful to probe the star-forming properties of galaxies, as galaxy populations typically separate out into the star-forming ‘blue cloud,’ and the relatively quiescent ‘red sequence.’ At z > 6, when the universe was less than 1 Gyr old, our HST observations probe only the rest-frame UV, thus our colour diagnostics are somewhat more limited than at lower redshift (though we don’t expect a large fraction of quiescent galaxies, e.g., Muzzin et al. Reference Muzzin2013; Nayyeri et al. Reference Nayyeri2014).

However, the rest-frame UV can provide robust constraints on the dust content in distant galaxies. While one would prefer to directly probe the metallicities of such systems to track chemical enrichment, in practice, this is difficult with current technology. However, as dust grains are made from the heavy elements which form in stars, the build-up of dust in the universe tracks the build-up of metals. Additionally, although dust is not unique in reddening the rest-frame UV colours of galaxies, the change in UV colour with increased dust content is larger than that with a comparable increase in metallicity or stellar population age (particularly at z > 6, as stars must be younger than the age of the universe). Thus, tracking the evolution of the UV colours of distant galaxies provides an excellent path for following the evolution of the dust content in these systems, therefore providing a proxy to track the build-up of heavy metals in the universe.

Rather than using a single pair of filters, the rest-frame UV colour is typically parameterised by the UV spectral slope β, which is defined as f _λ∝λ^β (Calzetti, Kinney, & Storchi-Bergmann Reference Calzetti, Kinney and Storchi-Bergmann1994). The quantity β is an excellent tracer of dust extinction, as it has been found to be well correlated with the ratio of the far-infrared to UV flux at both low and moderate redshift (e.g., Meurer et al. Reference Meurer, Heckman and Calzetti1999; Reddy et al. Reference Reddy2012), though the exact shape of this correlation depends on the dust attenuation curve (e.g., Capak et al. Reference Capak2015) and other factors (e.g., Wilkins et al. Reference Wilkins, Gonzalez-Perez, Lacey and Baugh2012, Reference Wilkins2013). Earlier work showed that β became substantially more negative (i.e., bluer UV colours) with increasing redshift, from β ~ −1.5 at z ~ 2 to β ~ −2 at z ~ 6 (e.g., Bouwens et al. Reference Bouwens2009).

6.2. The growth of galaxies

6.2.1. Galaxy stellar masses

The distribution function of galaxy stellar masses is also a useful tool to probe the galaxy physics discussed in Section 5.1. In particular, comparing the slope of the low-mass-end of the stellar mass function to that of the underlying dark matter halo mass function can provide an impression of the impact of feedback (supernovae, or stellar radiative) on the suppression of star formation. This has the advantage over the UV luminosity function in that it directly probes an intrinsic galaxy property without the troubling effect of dust. However, the measurement of galaxy stellar masses is an inferred quantity, therefore the galaxy stellar mass function is a less direct quantity to measure than the UV luminosity function. I will not discuss here the detailed methods for the measurements of galaxy stellar masses, but a recent review is available in Conroy (Reference Conroy2013). In general, with decent photometry, galaxy stellar masses can be measured accurately to within a factor of a few, typically a higher accuracy than other photometrically determined properties (though a changing initial mass function can result in a much larger systematic uncertainty in the stellar mass).

The optical and near-infrared imaging typically used to discover distant galaxies probes only the rest-frame UV, and is thus dominated by the most massive stars present in these galaxies. Stellar masses measured using only these filters can be subject to the ‘outshining’ phenomenon, where an older, perhaps more massive population is ‘outshone’ by a more recent episode of star formation at the observed wavelengths. Robust stellar mass measurements thus require longer wavelength imaging. While rest-frame near-infrared imaging would be desirable, this falls longward of 10 μm at z > 6, and is thus presently inaccessible. However, rest-frame optical imaging can still constrain moderately older populations, and is supplied by Spitzer/IRAC imaging at 3–5 μm at these redshifts. Direct detections at these wavelengths are not necessarily required, as if the IRAC upper limits occur at similar flux levels as the optical/near-infrared detections, useful constraints on the stellar mass can still be obtained (by, e.g., ruling out a massive older population, which would exceed the IRAC upper limits).

Individual galaxy stellar masses were determined for galaxies at z = 6 from a variety of studies. In particular, Yan et al. (Reference Yan2005, Reference Yan2006) and Eyles et al. (Reference Eyles2007) studied i-band dropout galaxies in the HUDF and/or GOODS fields which had clear IRAC detections, and found that they had stellar masses of ≳ 10¹⁰ M_⊙, surprisingly high for only ~ 1 Gyr after the Big Bang. Stark et al. (Reference Stark2009) also measured stellar masses for a sample of a few dozen z ~ 6 galaxies. They were not all detected in IRAC, thus the typical stellar masses were lower, ~ 10⁹ M_⊙ for galaxies with M _UV = −20 to − 21, though there was a significant scatter with some galaxies as massive as 10¹⁰ M_⊙. These massive galaxies with IRAC detections were fit with strong 4 000 Å breaks, indicative of a dominant population of older stars (~ 500 Myr; e.g., Eyles et al. Reference Eyles2007). At higher redshift, Finkelstein et al. (Reference Finkelstein2010) measured the stellar masses for ~ 30 z = 7–8 galaxies in the HUDF (though see Egami et al. Reference Egami2005 for an earlier mass measurement of an individual lensed z ~ 7 candidate galaxy). As these galaxies were typically faint, they were mostly not detected in IRAC, and had correspondingly low stellar masses of ~ 10⁹ M_⊙, an order of magnitude less massive than that of an L* galaxy at z = 3, though intriguingly similar to the masses of narrowband-selected LAEs at z = 3–6. Labbé et al. (Reference Labbé2010) found similar stellar mass results in a stack of z = 7–8 galaxies in the HUDF, but through a stacked analysis was able to measure a significant average IRAC flux for these galaxies, inferring a somewhat surprising age of 300 Myr, implying that the majority of stars in these galaxies formed prior to z = 10.

One complicating issue to the inclusion of IRAC photometry in these fits is the potential for a significant contribution from nebular emission lines. In particular, if the [O iii] λλ4959,5007 and/or Hα emission lines were strong enough, they could add significantly to the bandpass-averaged flux in the IRAC bands. This was initially not thought to be a problem, as the majority of low-redshift star-forming galaxies do not show extremely strong lines, and thus these lines were not accounted for in these initial studies. However, a bevy of observational results have since shown that high-EW [O iii] and Hα emission is common, if not ubiquitous amongst distant star-forming galaxies, possibly indicating changes in the physical environments of typical star-forming regions (e.g., Schaerer & de Barros Reference Schaerer and de Barros2009, Reference Schaerer and de Barros2010; Finkelstein et al. Reference Finkelstein, Cohen, Malhotra and Rhoads2009, Reference Finkelstein2011; Shim et al. Reference Shim2011; Stark et al. Reference Stark2013; Labbé et al. Reference Labbé2013; Song et al. Reference Song2014, Reference Song2016a; Salmon et al. Reference Salmon2015; Smit et al. Reference Smit2015; Finkelstein et al. Reference Finkelstein, Finkelstein, Malhotra and J.2015a). These strong lines boost the IRAC flux, mimicking the presence of a strong 4 000 Å break, causing one to infer a higher mass and/or older population than may be correct. Schaerer & de Barros (Reference Schaerer and de Barros2009, Reference Schaerer and de Barros2010) showed that by including nebular lines in the stellar population modelling process much lower stellar masses and stellar population ages can be obtained, with ages as young as 4 Myr possible for the stacked photometry of faint z ~ 7 galaxies from Labbé et al. (Reference Labbé2010). Salmon et al. (Reference Salmon2015) quantified this further, finding that when including nebular emission lines, the stellar masses of z ~ 6 galaxies in the CANDELS/GOODS-S field were systematically measured to be ~ 20–30% lower than when emission lines were ignored.

Turning to the distribution of stellar masses, the first full, completeness-corrected stellar mass functions at z ⩾ 6 were published by González et al. (Reference González2011). They found low-mass-end slopes which were surprisingly shallow, with α = −1.44 at z = 6 and − 1.55 at z = 7, compared to simulations at the time, which often predicted α < −2 (e.g., Jaacks et al. Reference Jaacks, Choi, Nagamine, Thompson and Varghese2012a). The observations thus implied that the implementation of feedback in the simulations was possibly too weak. However, these initial results were based on small samples of only ~ 60 galaxies total at z = 6 and 7 combined, and also did not include the impact of nebular lines. Recently, three separate studies on the stellar mass function have been completed using a combination of different subsets of the CANDELS data with the HUDF. Duncan et al. (Reference Duncan2014) used the CANDELS GOODS-S field to find much steeper faint-end slopes of ~ −1.9 at these redshifts. Grazian et al. (Reference Grazian2015) used both the CANDELS GOODS-S and UDS fields to find a similarly steep slope at z = 7, though they found α = −1.55 at z = 6. The most recent study, by Song et al. (Reference Song2016a), also used two CANDELS fields (GOODS-N and GOODS-S), but made use of deeper IRAC data newly available from the S-CANDELS survey (PI Fazio). By using simulations to verify that their methods minimised systematic offsets in their results, Song et al. (Reference Song2016a) found that the low-mass-end slope α steadily becomes steeper, from − 1.53 at z = 4 (similar to González et al. Reference González2011) to − 2.05 at z = 7. At the highest redshifts, this slope is similar to that of the dark matter halo mass function, implying that feedback may be less efficient at suppressing star formation towards higher redshift.

The integral of the stellar mass function provides the stellar mass density of the universe at a given epoch. As measurements at very low-redshift from, e.g., the SDSS, are quite robust, it is interesting to compare the stellar mass density at high redshift to low, to investigate what fraction of the total stellar mass of the universe exists at a given epoch. This has been done in a variety of studies, (e.g., Dickinson, Giavalisco, & The Goods Team Reference Dickinson, Giavalisco, Bender and Renzini2003; Rudnick et al. Reference Rudnick2003; Duncan et al. Reference Duncan2014; Grazian et al. Reference Grazian2015; Oesch et al. Reference Oesch2015b), which are highlighted in Figure 8. In this figure, I also show the stellar mass density obtained by integrating the evolution of the SFR density from Figure 7, showing both the updated reference trend from this work, as well as the values from Madau & Dickinson (Reference Madau and Dickinson2014). At z < 2, the stellar mass density derived in this way is systematically higher by a small amount than that obtained directly from stellar mass measurements. Madau & Dickinson (Reference Madau and Dickinson2014) discussed a number of possible causes for this offset, including overestimation of the SFRs (which could be the case if, for example, galaxies have a dust attenuation curve similar to the Small Magellanic Cloud), or underestimation of the stellar masses (due to lack of constraints on older stars, and/or an evolving IMF). However, at z > 3, our integrated SFR density shows excellent agreement with the observations of Grazian et al. (Reference Grazian2015), Song et al. (Reference Song2016a), and Oesch et al. (Reference Oesch2014). The values from Duncan et al. (Reference Duncan2014) show a modest significant difference at z ⩾ 5 due to their much steeper low-mass-end slope (see discussion in Song et al. Reference Song2016a, and also Graus et al. Reference Graus, Bullock, Boylan-Kolchin and Weisz2015). By comparing the values of my model at a given redshift to that at z = 0, I find that the universe formed 50% of its stellar mass by z ~ 1.3, with ~ 10%, 1%, 0.1%, and 0.01% of the stellar mass in place by z = 2.9, 5.4, 8.0, and 10.2, respectively. Observations have therefore tentatively inferred the presence of 99.99% of all the stellar mass which has ever formed in the universe, though much of it remains to be directly detected.

Figure 8. The evolution of the total stellar mass density in the universe, all derived assuming a Salpeter IMF. The low redshift measurements have a range of definitions, but I note that all high redshift measurements were obtained by integrating stellar mass functions from 8 $< \log\, M/\text{M}$ _⊙ < 13 (the flatter slope of the low-mass end of the stellar mass function at z < 4 implies that the lower limit of the integration is less important). The gray points show the data from the compilation of Madau & Dickinson (Reference Madau and Dickinson2014), while the other symbols come from recent estimates of the stellar mass function at high redshift, defined in the legend. The blue and dashed black curves show the stellar mass density obtained by integrating the SFR density evolution from Figure 7. The left side of the shaded regions denote the redshift at which the total stellar mass density formed is equal to the listed percentage of the z = 0 value.

6.2.3. Evolution of galaxy specific star-formation rates

The specific star-formation rates (sSFR; SFR/stellar mass) of galaxies are an excellent tracer of galaxy growth, as this quantity normalises the growth of new stars by the current stellar mass. Theoretical simulations predict that the sSFR at fixed stellar mass should be rising with increasing redshift, as such models predict that star formation is governed by cold gas accretion, and should therefore increase with increasing redshift (e.g., Neistein & Dekel Reference Neistein and Dekel2008; Davé et al. Reference Davé, Finlator and Oppenheimer2011; Krumholz & Dekel Reference Krumholz and Dekel2012). It was thus surprising when the first measurements of the evolution of the sSFR at high redshift showed it to be roughly flat from z ~ 3 to 7 (González et al. Reference González2010). Subsequent work by Stark et al. (Reference Stark2013) considered the impact of nebular emission, which was not accounted for in González et al. (Reference González2010), and found hints that the sSFR may be rising from z = 5 to 7. González et al. (Reference González2014) performed an updated analysis with a larger sample, including nebular lines and varying SFHs (including those which increase with time) and also found evidence for a rising sSFR with increasing redshift, though with a redshift dependance much weaker than simulations.

Most recently, Salmon et al. (Reference Salmon2015) found a steeper increase in sSFR with increasing redshift from z = 4–6. Their fiducial analysis tracking galaxies at constant mass was not a great match to simulations, but when they instead considered tracking galaxies with an evolving cumulative number density to more accurately track progenitors and descendants (Behroozi et al. Reference Behroozi2013a), they found an excellent match to both Neistein & Dekel (Reference Neistein and Dekel2008) and Davé et al. (Reference Davé, Katz, Oppenheimer, Kollmeier and Weinberg2013) in both normalisation and slope. Therefore, at present, there is no significant discrepancy with models at high redshift. However, current results at z > 4 have large uncertainties, and, at the highest redshifts, are based on small samples, thus there is room for improvement. Additionally, there is a significant discrepancy between models and observations at 0.5 < z < 2 (see Figure 9 of González et al. Reference González2014), where observed sSFRs evolve shallowly down to z ~ 1 then fall off quickly towards lower redshift, while simulations predict a smoother evolution. Thus, the evolution of the sSFR will likely remain an active area of inquiry at all redshifts.

Figure 9. The volume-averaged neutral fraction as a function of redshift, with the 68% confidence range on constraints from the integral of the reference UV luminosity functions shown by the shaded blue region. An upper limit from the Lyα narrowband survey by Ouchi et al. (Reference Ouchi2010) is shown by the blue circle, while constraints via Lyα spectroscopy at z = 6.5–7 by Pentericci et al. (Reference Pentericci2014) and at z = 7.5–8 by Tilvi et al. (Reference Tilvi2014) are shown by the blue arrows. Constraints via quasars are shown by the magenta regions (from Fan et al. Reference Fan2006 at z < 6, and Bolton et al. Reference Bolton2011 at z = 7). The Thomson optical depth to electron scattering derived from the reference UV luminosity function constraints is shown by the hatched blue region, and is consistent with constraints from Planck, shown by the hatched red region (both corresponding to values on the right-hand vertical axis), as well as the recent study by Robertson et al. (Reference Robertson, Ellis, Furlanetto and Dunlop2015) shown in green. This fiducial reionisation history constrains reionisation to be < 10% complete by z ~ 10, and > 85% complete by z ~ 7.

7.1. The timeline of reionisation

7.2. Constraints from galaxies

The declining quasar luminosity function at z > 3 (e.g., Richards et al. Reference Richards2006) has led to the conclusion that star-formation in galaxies produced the bulk of ionising photons necessary for reionisation. A key measure is thus a cosmic ‘census’ of distant galaxies, measuring their rest-UV light, then inferring the likely ionising photon emissivity to see if the observed star formation at a given epoch can sustain an ionised IGM. If one assumes that galaxies do indeed provide all the needed ionising photons, this analysis provides an independent measure into the timeline of reionisation.

There are a number of assumptions which need to be made to convert rest-frame UV observations of galaxies into a constraint on reionisation.

• • Stellar population: As one observes non-ionising UV photons, that observable needs to be converted to the emergent ionising radiation, which requires assuming a stellar population. The primary considerations here are the initial mass function and the metallicity. The conversion factors do not vary wildly if one assumes a Salpeter IMF and non-zero metallicity (~ 20–30%; e.g., Finkelstein et al. Reference Finkelstein2012a). However, a changing IMF, particularly one which results in more massive stars, could produce a much higher number of ionising photons per unit UV luminosity making it easier for galaxies to reionise the IGM.

• • Escape fraction: Only ionising photons which escape a galaxy are available to contribute to reionisation, thus one needs to assume an escape fraction of ionising photons (f $_{\text{esc}}$ ). This quantity cannot be directly observationally ascertained in the epoch of reionisation due to attenuation along the line-of-sight. Studies at z < 4 have found that the majority of galaxies have little-to-no escaping ionising radiation (e.g., Siana et al. Reference Siana2010), though examples of escape fractions on the order of 20% have been found (e.g., Nestor et al. Reference Nestor, Shapley, Steidel and Siana2011). Finkelstein et al. (Reference Finkelstein2012a) showed that values greater than 13% at z = 6 would violate the ionising photon emissivity obtained from observations of the Lyα forest. Most simulations predict that galaxies massive enough to be currently observed have very low escape fractions, and values > 20% only occur in the very lowest mass halos (e.g., Paardekooper et al. Reference Paardekooper, Khochfar and Dalla Vecchia2015). Values of ~ 20% are frequently assumed in the literature when obtaining reionisation constraints from galaxies, which may thus be optimistically too high (e.g. Finkelstein et al. Reference Finkelstein2010; Robertson et al. Reference Robertson2013, Reference Robertson, Ellis, Furlanetto and Dunlop2015).

• • Limiting Magnitude: When integrating the UV luminosity function, one needs to assume a limiting magnitude. An absolute limit would be the magnitude of a single O star (~ −5); however, galaxies likely have their star formation suppressed in more massive halos, due to photoionisation of their gas from star formation in nearby halos. Galaxies are commonly observed at z ~ 6–8 down to magnitudes of M _UV = −17, thus it may be reasonable to assume fainter values of − 13 or − 10, as is often done in the literature. Recent simulations have found that the luminosity function may break at M = −14 to − 16 due to a combination of physical processes (Jaacks et al. Reference Jaacks, Thompson and Nagamine2013; O’Shea et al. Reference O’Shea, Wise, Xu and Norman2015), or at − 13 by simulating the progenitors of galaxies in the Local Group (e.g., Boylan-Kolchin et al. Reference Boylan-Kolchin2015). However, initial results from the Hubble Frontier Fields lensing survey show that the luminosity function continues steeply to M _UV = ~13 at z = 6 (Livermore et al. Reference Livermore, Finkelstein and Lotz2016), and future studies with the full Frontier Fields dataset should produce even better constraints.

• • Clumping Factor of IGM: The density of gas in the IGM is a key parameter, as it is more difficult to keep dense gas ionised than diffuse gas, as the recombination time is shorter. This is typically parameterised by the clumping factor (C), where higher values of C require more ionising photons to sustain reionisation. While C was originally thought to be quite high (~ 30; Gnedin & Ostriker Reference Gnedin and Ostriker1997), more recent models, with a better understanding of the interface between galaxies and the IGM, find lower values (< 5; e.g., Pawlik, Schaye, & van Scherpenzeel Reference Pawlik, Schaye and van Scherpenzeel2009; Finlator et al. Reference Finlator, Oh, Özel and Davé2012; Pawlik, Schaye, & Dalla Vecchia Reference Pawlik, Schaye and Dalla Vecchia2015). The most recent models of Pawlik et al. (Reference Pawlik, Schaye and Dalla Vecchia2015) find a clumping factor which is < 3 at z > 8, rising to C = 4 at z = 6. Most studies below assume C = 3.

With these assumptions, one can use a model of reionisation to determine whether the galaxy population at a given redshift can sustain an ionised IGM. The model of Madau, Haardt, & Rees (Reference Madau, Haardt and Rees1999) is commonly used, which calculates a redshift-dependent limiting UV luminosity density to sustain an ionised IGM, assuming a stellar population term, and values for f $_{\text{esc}}$ and C. One can then integrate an observed luminosity function (to an assumed limiting magnitude), and compare to this model to explore the contribution of galaxies to reionisation. This has been done in a variety of studies; here, I will review a few recent results.

Using new samples of z = 6, 7, and 8 galaxies from the CANDELS and HUDF surveys, Finkelstein et al. (Reference Finkelstein2012a) examined the contribution to reionisation solely from observed galaxies (i.e., with no extrapolation to fainter magnitudes). They found that if the typical galaxy had an escape fraction of ~ 30%, then the observed population of galaxies could sustain an ionised IGM by z = 6, in contrast to previous studies which found that fainter galaxies were required (e.g., Bunker et al. Reference Bunker, Stanway, Ellis and McMahon2004; Yan & Windhorst Reference Yan and Windhorst2004). By using limits on the ionising emissivity from observations of the z ≈ 6 Lyα forest, Finkelstein et al. (Reference Finkelstein2012a) found that if the UV luminosity function continues with the observed steep slope to M _UV = −13, then the typical escape fraction could be no higher than ~ 13%. However, even with this modest escape fraction, galaxies could still sustain an ionised IGM by z = 6 if the luminosity function extended down to such faint magnitudes. Robertson et al. (Reference Robertson2013) performed a complementary analysis, using not only the UV luminosity function, but also measurements of the stellar mass density and the electron scattering optical depth from the WMAP measurements of the CMB as constraints. They found similar results, in that by assuming f $_{\text{esc}} =$ 20%, galaxies can sustain an ionised IGM by just after z = 6, but that fainter galaxies would be required to complete reionisation at earlier times.

Over the next few years, the advent of Planck coupled with the large improvement on the precision of the UV luminosity function from the completed CANDELS survey allowed these constraints to be tightened. Finkelstein et al. (Reference Finkelstein2015c) used their updated luminosity functions to revisit the constraints placed on reionisation from galaxies. Assuming a fiducial model with C = 3 and f $_{\text{esc}} =$ 13%, they found that galaxies suggest that reionisation was largely complete by z = 6, and that the universe was primarily neutral by z > 10. Their observations constrained a midpoint of reionisation (defined as the epoch when the volume-averaged ionised fraction reached 50%) of 6.7 < z < 9.4 (at 68% confidence). Although this analysis did not use the recent observations from Planck as constraints, the electron scattering optical depth derived from this reionisation history is in excellent agreement with the Planck observations, implying that a significant contribution to reionisation from galaxies at z > 10 is not needed. Robertson et al. (Reference Robertson, Ellis, Furlanetto and Dunlop2015) updated their previous analysis with the Planck constraints and the latest luminosity function results, and found a similar result in that a reionisation process which starts at z ≈ 10 and completes by z ≈ 6 is consistent with both the observations of galaxies and the Planck observations of the CMB. Finally, Bouwens et al. (Reference Bouwens2015a) recently performed a complementary analysis, eschewing constraints from the luminosity functions of galaxies in their analysis to derive independant constraints on the ionising sources. As constraints, they combined observations of the Lyα forest from quasar spectra with the Planck optical depth measurements, as well as observations of the evolving Lyα emission from galaxies. They derived an evolution in the ionising background which was consistent with the ionising emissivity derived from observations of galaxies (with $f_{\text{esc}} \sim$ 10–20%) even though galaxies were not used as constraints. This again suggests that galaxies are the major producers of the ionising photons which drove the reionisation of the IGM, provided such ionising escape fractions are achievable.

7.3. Reference reionisation constraints

The studies discussed above uniformly arrive at a similar conclusions: galaxies were responsible for reionisation, and reionisation completed by z ~ 6. However, the time-evolution of the neutral fraction during reionisation is less well constrained. For example, the constraints on the mid-point of reionisation from Finkelstein et al. (Reference Finkelstein2015c) of Δz = 2.7 are largely due to the faint-end slope uncertainties from the luminosity function. As discussed in Section 5.3, I can improve upon these constraints with my derived reference UV luminosity function. For example, the fractional uncertainty on the faint-end slope (σ_α/α) from Finkelstein et al. (Reference Finkelstein2015c) and Bouwens et al. (Reference Bouwens2015b) was 10 and 6%, respectively, while from my reference UV luminosity function at z = 7, I find σ_α/α = 2%. Therefore, to conclude this discussion of reionisation, I explore the improvements in the constraints on the contribution of galaxies to reionisation available when I use luminosity functions derived in Section 5.

To derive this reionisation history, I follow the methodology of Finkelstein et al. (Reference Finkelstein2015c), assuming C = 3 and f $_{\text{esc}} =$ 13% (which does not violate constraints on the Lyα forest at z = 6). I then integrate the fiducial reference UV luminosity functions from Section 5.3 at each redshift to a limiting magnitude of M _UV = −13 to derive the UV luminosity density, using the results from the MCMC chains to derive uncertainties on these values. Using the model of Madau et al. (Reference Madau, Haardt and Rees1999) with the above assumptions, I then calculate the volume-averaged neutral fraction (x _HI) at each redshift, propagating through the uncertainties on the UV luminosity density to derive a 68% confidence range on x_HI. The resultant 68% confidence reionisation history is shown in Figure 9.

The tighter constraints on the luminosity functions available when using all available data lead to more constraining results on the evolution of the volume averaged neutral fraction x $_{\text{HI}}$ than previous studies utilising only galaxy data. This analysis constrains x _HI = 0% at z = 6 and x _HI > 80% at z = 10, with the 68% confidence constraint on the midpoint of reionisation being 7.3 < z < 7.7. I calculated the Thomson optical depth to electron scattering from both my upper and lower 68% confidence constraints on the neutral fraction, finding $\tau _{\text{es}} =$ 0.057 – 0.060 (blue hatched region in Figure 9). This is in consistent with the low end of the observational constraints from Planck of $\tau _{\text{es}} =$ 0.066 ± 0.012. This is particularly striking as the Planck measurements were not used as an input constraint on the galaxy luminosity function-based reionisation history. Rather, the galaxy observations, combined with assumptions on the limiting magnitude, clumping factor, and escape fraction, produce a result which is consistent with the Planck observations. The results found here are highly consistent with those from Robertson et al. (Reference Robertson, Ellis, Furlanetto and Dunlop2015), as shown by the green-shaded region in Figure 9. If the true value of $\tau _{\text{es}}$ is closer to the higher end of the Planck constraints, it could leave room for modest star-formation activity at z > 10, while if it is closer to the value implied by our observations at the lower end of the Planck constraints, it leaves little room for a significant number of free electrons at higher redshift.

Comparisons to other complementary observations are also shown in Figure 9. Current results from quasars imply a fully ionised universe at z < 6 (e.g., Fan et al. Reference Fan2006), while the single known quasar at z = 7 implies a volume averaged neutral fraction > 10%. On the other hand, Ouchi et al. (Reference Ouchi2010) found x_HI < 20% at z = 6.6, using the observed evolution of the Lyα luminosity function via narrowband imaging. Using the incidence of Lyα emission in spectroscopic follow-up of known high-redshift sources, Pentericci et al. (Reference Pentericci2014) and Tilvi et al. (Reference Tilvi2014) constrained 30% < x _HI < 90% at z = 7, and x _HI > 30% at z = 8, respectively. While both the quasar and Lyα luminosity function-derived constraints are fully consistent with my inferred reionisation history, taken at face value the z = 7 Lyα results are in modest tension. These studies imply a significant neutral fraction of > 30% by z = 7, while the luminosity functions seemingly allow at most 20% of the IGM to be neutral.

However, there is presently enough uncertainty in all of these measures such that the tension can be easily resolved. For example, reducing the assumed escape fraction to 10% increases the neutral fraction at z = 7 to be consistent with the lower limit from Pentericci et al. (Reference Pentericci2014), while still generating enough free electrons to satisfy the Planck constraint. Likewise, the constraints via Lyα spectroscopy are relatively poor, due both to small sample sizes as well as modelling uncertainties, both which should improve significantly in the coming years as observers build up more robust samples of z > 7 spectroscopic observations. The escape fraction is likely the biggest wild card, and while here I impose an upper limit set by Lyα forest constraints at z = 6, we have no reasonable constraint on the lower limit of this quantity. Future studies at z = 2–3, as well as more advanced modelling, should shed light on this issue in the coming years. Should the escape fractions be universally lower than currently assumed, it may be that galaxies play less of a role in reionisation than is currently thought, opening the door for a larger contribution from low-luminosity AGN (e.g., Giallongo et al. Reference Giallongo2015; Madau & Haardt Reference Madau and Haardt2015).

8.1. Probing the dark side with ALMA

Much of the observational studies above have focussed on observations of the starlight from galaxies. While this is a direct probe into galaxy evolution, it leaves out a key component, which is the ubiquitous gas which is fueling the star formation we observe. Observations at lower redshifts imply that galaxies in the distant universe are substantially more gas-rich than at low redshift (e.g., Tacconi et al. Reference Tacconi2010; Papovich et al. Reference Papovich, Finkelstein, Ferguson, Lotz and Giavalisco2011, Reference Papovich2015). This could have strong implications on the evolution of such distant galaxies, and indeed there are hints, discussed above, that the star-forming properties of distant galaxies may be different than expected, perhaps due to large gas reservoirs. The advent of ALMA now enables direct studies into the gas properties of distant (z > 6) galaxies. While semi-direct probes of molecular gas in normal galaxies such as CO emission may be observationally highly expensive (e.g., Tan et al. Reference Tan2013), more indirect probes such as the dust emission (e.g., Scoville et al. Reference Scoville2016) and/or dynamical modelling via resolved [C ii] emission can shed light on gas at high redshift. While ALMA observations have only just begun, recently improved programmes such as the mapping of the HUDF (PI Dunlop), and of a portion of the GOODS-S field (PI Elbaz) should begin to allow us to probe the dark side of the distant universe. A better understanding of the gas properties of these distant galaxies will shed light on the physical mechanisms fueling galaxy evolution, perhaps explaining observed trends in galaxy sSFRs and the possibly evolving SMHM relation.

8.2. Action items for JWST

The launch of JWST is now less than 3 yr away, thus the time is now for identifying the key issues for study immediately following the launch. While there are clearly many important issues that JWST will tackle at z > 6, here I focus on a few key issues. The first issue is the faint-end slope of the galaxy UV luminosity function. Although the reference luminosity functions derived here have improved the constraints on the faint-end slope at z = 7, there remains a large uncertainty in the shape (and extent) of the faint-end at 8 < z < 10, which leaves possible a wide range of possible IGM neutral fractions, as shown in Figure 9. A deep JWST survey with a time investment similar to the HUDF, will reach UV absolute magnitudes as faint as − 16 (− 15.5) at z = 10 (z = 7), two magnitudes fainter than the current HST observations (without the addition of magnification uncertainties). Not only will this allow much better constraints on the slope of the faint-end of the luminosity function, but it will also probe whether the single power law slope is valid down to such faint magnitudes, which is not the case in some models (Jaacks et al. Reference Jaacks, Thompson and Nagamine2013; O’Shea et al. Reference O’Shea, Wise, Xu and Norman2015).

Such a survey will also be capable of detecting galaxies out to z > 12, and perhaps z = 15 (depending on their brightness, and on the filters used). A uniform sample of galaxies from 7 < z ≲ 12 will allow a robust investigation into the evolution of the cosmic SFR density. This analysis will answer the question of whether there is a steep decline at z = 8, as not only will it discover larger numbers of moderately bright galaxies at 8 < z < 10, but it will probe much further down the luminosity function, allowing us to gain a much greater confidence on the amount of total star formation occurring, constraining galaxy formation at the earliest epochs.

The infrared sensitivity of JWST will allow a more robust investigation into the stellar mass density of the distant universe. The high-redshift data shown in Figure 5 are all based on rest-frame UV selected samples, thus the stellar mass measured is only for galaxies which have been forming stars at some point in the ~ 100 Myr prior to observation. If the duty cycle of star formation is not near unity, it is possible that there are populations of galaxies between bursts of star formation which are not accounted for in the current observations. Such galaxies, however, still host lower mass stars, and thus would be visible in a longer wavelength survey. A deep survey at observed 5 μm would measure rest-frame 0.6 μm light from z = 7 galaxies, allowing the detection of such massive UV-faint systems. The stellar mass function computed from such a survey would represent a more robust census of the total stellar mass formed at such redshifts, and will also allow a truly direct investigation into the duty cycle of star formation at such redshifts.

Finally, spectroscopy with JWST will clearly be an enormous discovery space. JWST will observe Lyα to much fainter line flux levels than is currently possible, both increasing the number of spectroscopically confirmed distant galaxies, and also providing tighter constraints on reionisation via Lyα emission. Second, observations of rest-optical emission lines such as [O iii], which are already implied to be quite bright from strong Spitzer colours, will provide an easier route to spectroscopic confirmation, and also allow probes into the chemical enrichment and ionisation conditions inside these galaxies.

8.3. Future space-based telescopes

The next decade will bring the Wide-Field Infrared Survey Telescope (WFIRST), which will sport a 0.28 deg² camera housed on a 2.4-m space telescope. While such a telescope will have a similar flux sensitivity as HST, the wide field will allow a number of unique investigations into the distant universe. The planned 2 227 deg² High Latitude Survey (HLS) to a near-infrared depth of 26.5 AB mag will allow a truly unique view into the evolution of the bright-end of the galaxy UV luminosity function at very high redshifts. When paired with optical data from the Large Synoptic Survey Telescope (LSST), this survey will enable the selection of galaxies at 4 < z < 10 over this enormous area, approximately 10⁴ times the size of the CANDELS survey. WFIRST will also have a General Observer component, and as this telescope is similar in size to HST, a unique programme would be to perform an extremely deep survey over a single square degree, which would encompass an area ~ 720 × that of the HUDF. While deep JWST surveys will likely provide tight constraints on the faint-end slopes at 4 < z < 10 by this time, they will still be subject to significant cosmic variance uncertainties due to the small field-of-view of the JWST instruments. By probing to 29 AB mag over a very wide field, these uncertainties can be mitigated.

In Table 3, I list the expected number counts and source densities based on the reference UV luminosity functions, for both the WFIRST HLS and a hypothetical 1 deg² m _AB = 29 WFIRST survey. At z = 7, over one hundred thousand m < 27 galaxies will be discovered, > 1 000 × the number of bright z = 7 galaxies currently known. This will allow the relatively rough constraints on the bright-end of the high-redshift UV luminosity functions to be dramatically improved, allowing much more detailed investigations into the issues discussed in Section 5.1, such as the evolution of the impact of feedback on distant galaxies. In particular, the cosmic variance which currently plagues the bright-end of the galaxy UV luminosity function (Bouwens et al. Reference Bouwens2015b; Finkelstein et al. Reference Finkelstein2015c) will be mitigated. Combining the cosmic variance estimator of Newman & Davis (Reference Newman and Davis2002) for dark matter halos with the measured bias at z = 7 from Barone-Nugent et al. (Reference Barone-Nugent2014), Finkelstein et al. (Reference Finkelstein2015c) found an uncertainty in the number of bright z = 7 galaxies to be ~ 20%. Using these same methods, the WFIRST HLS will have a fractional uncertainty due to cosmic variance of only 0.2%. For faint galaxies, Finkelstein et al. (Reference Finkelstein2015c) found a cosmic variance uncertainty of ~ 40%, while the analysis here shows that a deep 1 deg² survey (discovering galaxies to m = 28.5) would reduce this uncertainty to only 8%. We will get a taste of this science early in the next decade, when the smaller Euclid launches, which will perform a 40 deg² deep survey ~ 1 mag shallower than the WFIRST HLS.

Table 3. Predicted source counts and densities for WFIRST HLS and deep GO surveys.

The reference luminosity functions derived in this review were used for these predictions. Columns 3 and 4 refer to the planned 2 227 deg² WFIRST High Latitude Survey, while columns 5 and 6 refers to a hypothetical deep 1 deg² survey. To derive these predictions, I assumed a depth of the WFIRST HLS of 26.5 AB in the Y, J, and H filters. Assuming a desired detectable Lyman break of at least 0.5 mag this implies a limiting magnitude for discoverable z = 8, 9, and 10 galaxies of m _AB = 26 for the HLS. Galaxies at z = 7 will be limited by the depth of the LSST z′-band, which is estimated to be 26.2³, limiting z = 7 galaxy studies in the WFIRST HLS to m < 25.7. The bluer LSST bands will go much deeper, thus I use a limiting magnitude of 26.5 for z = 4–6, such that these galaxies are still detected in the WFIRST HLS.

The somewhat more distant future (likely in the 2030s) should see the advent of a 12–14-m class space telescope. Several science case and design studies are underway for concept such as the High Definition Space Telescope (HDST; Dalcanton et al. Reference Dalcanton2015) and the Advanced Technology Large Aperture Space Telescope (ATLAST; Postman et al. Reference Postman, Foresto, Gelino and Ribas2010). Such a telescope will likely be designed for the discovery of habitable worlds around other stars, but the technical requirements will result in an observatory with direct applications to the distant universe. One particular goal of such an observatory would be to account for all of the ionising photons needed for reionisation. As discussed above, the current observations likely only account for a tiny fraction of the galaxies which power reionisation. As shown in Figure 10, even with a steep luminosity function, JWST deep surveys will still only account for 63% of the necessary photons. A 12–14-m space telescope, reaching a limiting UV magnitude of − 14 at z = 7, will account for 85% of the total ionising photon budget (albeit indirectly, through observations of non-ionising UV light). More crucially, these studies will allow a direct investigation into the shape of the UV luminosity function, removing the need for the currently applied extrapolations to levels well below that which we can observe.

Figure 10. The reference z = 7 luminosity function, highlighting the absolute magnitudes reached in deep surveys by HST, and future deep surveys by JWST and a 12–14-m space telescope. Assuming that the z = 7 luminosity function integrated down to M = −13 (with f $_{\text{esc}} =$ 13%) can complete reionisation, current HST observations only account for 30% of the needed ionising photons. While an Ultra Deep Survey with JWST would only double this, a similar allocation of resources with a 12–14-m class space observatory, such as HDST or ATLAST, would account for 85% of the needed photons.

8.4. The next generation of ground-based telescopes

Much as todays 8–10-m class telescopes provide detailed spectroscopic follow-up for sources discovered by HST, the next generation of 25–40-m class ground-based telescopes will allow us to study in detail the new discoveries made by JWST and future space-based telescopes. There are presently three such telescopes in development: the 25 m Giant Magellan Telescope (GMT), the Thirty Meter Telescope (TMT) and the 39-m European Extremely Large Telescope (E-ELT). These telescopes will be complementary in many ways. For example, the E-ELT will clearly have the largest light-gathering power, while the GMT will have the largest simultaneous field-of-view. These facilities will enable a number of new observational analyses for the z > 6 universe; here, I will focus on two promising future lines of inquiry directly relevant to the topics in this review.

First, high resolution spectroscopic follow-up of very high redshift gamma ray bursts (GRBs) offers the potential to study both galaxy ISMs and the state of the IGM at 6 < z < 10. Amongst the most intense explosions in the universe, GRBs are predicted to occur early in the Universe’s history, and have been observed at z = 8.2 (Tanvir et al. Reference Tanvir2009) and likely at z = 9.4 (Cucchiara et al. Reference Cucchiara2011), thus GRBs provide extremely bright pencil-beam flashlights into gas in the distant universe. By examining metal absorption line features in the otherwise featureless spectra of GRBs, one can directly study the abundance of heavy metals in the host galaxy of the GRB (e.g., Prochaska et al. Reference Prochaska, Chen, Dessauges-Zavadsky and Bloom2007). Likewise, GRBs offer a chance to probe the ionisation state and temperature of the IGM at z > 7. While quasars have been utilised for these analyses previously (e.g., Fan et al. Reference Fan2006), there is only a single quasar currently known at z = 7, and none at higher redshifts. Should a space-based gamma ray observatory capable of discovering such sources be in place by the next decade, instruments such as GMTNIRS (Jaffe et al. Reference Jaffe2014) available at first light on the GMT, or NIRES on the TMT and SIMPLE (Origlia et al. Reference Origlia2010) on the E-ELT will allow high-resolution spectroscopic follow-up into the state of the ISMs of galaxies and the IGM at early times.

A second promising future analysis is the study of Lyα emission at z > 6 to much fainter flux levels than currently possible. Current studies (primarily with MOSFIRE on Keck) can only reach low limiting Lyα EWs for the brightest known z > 7 galaxies, rendering it presently very difficult to place constraints on Lyα in galaxies fainter than 27 AB mag, which comprise the majority of the known z > 7 galaxy population. The 5–15 × gain in light gathering power by these future telescopes will allow much deeper line flux depths to be probed. For example, deep spectroscopy on the GMT is expected to reach 5σ limiting line fluxes of a few × 10⁻¹⁹ erg s⁻¹ cm⁻², with the TMT and E-ELT reaching even deeper.

A direct future application of such a facility would be to map Lyα across both ionised and neutral regions in the 6 < z < 10 universe. A future multi-object spectrograph on the GMT should have a field-of-view of ~ 8′, comparable to the angular size of ionised bubbles during reionisation, and > 5 × larger in area than the equivalent spectrograph on JWST. While the imaging coverage needed to select galaxies over an area covering multiple such bubbles is likely beyond the realm of HST or JWST galaxy surveys, the WFIRST deep survey source densities in Table 3 show that such an instrument can simultaneously observe hundreds of 6 < z < 10 galaxies down to m = 28.5, allowing the use of Lyα emission to probe a wide range of IGM states during the epoch of reionisation. This will be highly complementary with 21-cm line surveys by, e.g., the Square Kilometer Array, as the combination of both surveys will allow a full picture of both the bright and dark side of reionisation.

8.5. Final words

The future is very bright for observational studies of galaxies in the first billion years after the Big Bang. As a community, we have made tremendous progress in not only discovering sources at such redshifts, but beginning to probe the physical processes dominating their evolution and gaining a better picture of how they affect the universe around them during reionisation. These tantalising glimpses into the nascent stages of our Universe leave us yearning for answers to the most pressing questions. When did the first stars form, and what were their characteristics? When did the true first galaxies form, and were their stellar populations substantially different than expected? Is reionisation truly over at z = 6, and was the neutral fraction substantial at z = 7, or is there some other effect hindering our view of Lyα? Finally, what is the interplay between gas and stars in such distant galaxies, and does star formation proceed differently? Although we may never truly be able to answer all of these questions, the next generation of astronomical observatories will allow us to make great progress in beginning to tease out the truth of the physics of the z > 6 universe. Combining the pursuit of these studies with the new questions we will uncover with these future facilities promises to make the coming decades an extremely exciting time for observational studies of the early universe.