Probing cosmic dawn: Ages and star formation histories of candidate z ≥ 9 galaxies

ABSTRACT

We discuss the spectral energy distributions and physical properties of six galaxies whose photometric redshifts suggest they lie beyond a redshift z ≃ 9. Each was selected on account of a prominent excess seen in the Spitzer/IRAC 4.5 |$\mu$|m band which, for a redshift above z = 9.0, likely indicates the presence of a rest-frame Balmer break and a stellar component that formed earlier than a redshift z ≃ 10. In addition to constraining the earlier star formation activity on the basis of fits using stellar population models with BAGPIPES, we have undertaken the necessary, but challenging, follow-up spectroscopy for each candidate using various combinations of Keck/MOSFIRE, VLT/X-shooter, Gemini/FLAMINGOS2, and ALMA. Based on either Lyman-α or [O iii] 88 |$\mu$|m emission, we determine a convincing redshift of z = 8.78 for GN-z-10-3 and a likely redshift of z = 9.28 for the lensed galaxy MACS0416-JD. For GN-z9-1, we conclude the case remains promising for a source beyond z ≃ 9. Together with earlier spectroscopic data for MACS1149-JD1, our analysis of this enlarged sample provides further support for a cosmic star formation history extending beyond redshifts z ≃ 10. We use our best-fitting stellar population models to reconstruct the past rest-frame UV luminosities of our sources and discuss the implications for tracing earlier progenitors of such systems with the James Webb Space Telescope.

galaxies: formation, galaxies: high-redshift, galaxies: stellar content

1 INTRODUCTION

An important goal of contemporary cosmology is to locate the epoch when the Universe was first bathed in starlight, popularly referred to as ‘cosmic dawn’. According to numerical simulations, dark matter haloes sufficiently massive to induce star formation could develop as early as 150–250 Myr after the big bang, corresponding to the redshift range z ≃ 15–20 (Wise et al. 2014; Villanueva-Domingo, Gnedin & Mena 2018). With some assumptions about their ionizing capabilities, studies based on the demographics of early galaxies (Robertson et al. 2015) have claimed to be compatible with independent constraints on early reionization based on analyses of the optical depth of electron scattering to the cosmic microwave background (CMB; Planck Collaboration VI 2020), as well as tentative measures of 21 cm absorption at a redshift z ≃ 16–19 in the CMB possibly associated with UV photons produced by early star formation (Bowman et al. 2018).

Recent deep imaging with the Hubble and Spitzer Space Telescopes has extended our redshift horizon to z ≃ 9–11 (Ellis et al. 2013; Oesch et al. 2018). At this redshift, the Universe is less than 600 Myr old and the presence of main sequence stars older than 250 Myr would imply galaxy formation originating before a redshift z ≃ 14. Such a case of ‘age-dating’ a high redshift galaxy was presented by Zheng et al. (2012) for the lensed galaxy MACS1149-JD1 based upon a prominent photometric excess seen in the Spitzer/IRAC 4.5 |$\mu$|m band. For the assumed photometric redshift of z = 9.6 ± 0.2, the IRAC excess would imply the presence of a Balmer break indicative of a burst of star formation some 200–300 Myr earlier. However, central to this conclusion is excluding the possibility of a redshift below z ≃ 9 where such an IRAC excess could arise from contamination of the 4.5 |$\mu$|m band by intense rest-frame optical [O iii] 5007 Å emission consistent with a much younger stellar population (Labbé et al. 2013; Roberts-Borsani et al. 2016). An early attempt to secure a spectroscopic redshift of MACS1149-JD1 using grism capabilities onboard the Hubble Space Telescope (HST) (Hoag et al. 2018) provided a tantalizing hint of a Lyman break at z > 9, however the lack of an emission line prevented a secure redshift constraint. A convincing redshift of z = 9.11 was eventually determined through Lyman-α and [O iii] 88 micron emission by Hashimoto et al. (2018), using ALMA and VLT/X-Shooter. The subsequent re-analysis of the object using the secure spectroscopic redshift and improved IRAC photometry resulted in an estimated age of 290 ± 150 Myr for the dominant stellar component, consistent with a formation redshift of 15.4 ± 2.3.

Various studies have discussed whether MACS1149-JD1 might be representative of the galaxy population at early times and hence provide a realistic first estimate of when cosmic dawn occurred. However, accurate interpretation of the associated photometry is not without its challenges. For instance, although Hashimoto et al. (2018) derived a stellar mass of ≃ 10⁹M_⊙ for MACS1149-JD1, this relied on an assumed and uncertain lensing magnification, thus introducing important uncertainties when comparing its properties to theoretical predictions. Furthermore, a fundamental difficulty has been that a prominent Balmer break requires a declining star formation rate with cosmic time for at least some portion of the early history of a galaxy, something that is difficult to reconcile with contemporary numerical simulations (Binggeli et al. 2019).

Observationally, several 7.5 < z < 9 objects (both lensed and unlensed) with secure spectroscopic redshifts display a significant IRAC excess which could be indicative of the presence of a Balmer break and early episodes of star formation (e.g. Tamura et al. 2019; Bakx et al. 2020; Strait et al. 2021). However, here the interpretation is confused by the likely prominence of [O iii] 5007 Å emission. To break the degeneracy of line emission and starlight as contributors to the IRAC excess, Roberts-Borsani, Ellis & Laporte (2020) combined HST and Spitzer/IRAC photometry with ALMA far-infrared [O iii] 88 |$\mu$|m emission and dust continuum constraints for four z > 7 sources (including MACS1149-JD1). Although their analysis gave a lower age for MACS1149-JD1 than Hashimoto et al. (2018), they none the less concluded that, collectively for these four sources, nearly half the stellar mass was produced prior to z ≃ 10.

Even for sources believed to be beyond z ≃ 9, it is possible that an IRAC excess may not entirely arise via a Balmer break. For MACS1149-JD1 Hashimoto et al. (2018) ruled out strong contamination of the 4.5 |$\mu$|m band from intense [O iii] 4959 Å and H β, which would have indicated a much younger system. Likewise, very strong [O ii] 3727 Å emission might contaminate the band, but only if the true redshift approaches z ≃ 10. Using hydrodynamical simulations and a comoving volume of 70 Mpc³, Katz et al. (2019) predicted the presence of three massive z > 9 galaxies with similarly strong IRAC colours and matched photometry to MACS1149-JD1. While their results also point to likely episodes of extended star formation, they conclude that dust may contribute to the IRAC excess, thereby reducing the strength of the Balmer break.

Given the implications for the birth of galaxies, determining the amount of star formation prior to z ≃ 9 is an important, but challenging endeavour. While ambiguities of interpretation may remain given the signal/noise of the available photometric data, our motivation in this paper is to further explore the Balmer break interpretation for six promising z ≥ 9 candidates (including MACS1149-JD1) each with a 4.5 |$\mu$|m IRAC excess. A major advance is a concerted attempt to secure their spectroscopic redshifts. Our goal is to explore the possible extent of star formation activity prior to a redshift z ≃ 10. We consider this timely given it may have important consequences for exploring this early period further with the James Webb Space Telescope.

A plan of the paper follows. In Section 2, we describe the selection criteria we applied to construct our sample of z ≥ 9 objects. In Section 3, we present the relevant spectral energy distributions (SEDs) and use these to infer the likely photometric redshifts and their star formation histories. This analysis provides a first assessment of the early star formation history for an enlarged population of z ≃ 9 galaxies. The results of our spectroscopic follow-up campaign are described in Section 4. Given the challenges of working at these limits, we are not able to convincingly secure a redshift for all of our candidates. Section 5 returns to a discussion of the physical properties of the sample in the light of the limited spectroscopic data. We also discuss the extent to which our IRAC excess sample may be representative of the broader population of z ≃ 9 galaxies and the frequency with which Lyman-α emission is being located deep in the reionization era. Finally, we discuss the implication of our findings in the context of direct searches for earlier activity with the James Webb Space Telescope in Section 6.

Throughout this paper, we use a concordance cosmology with H₀ = 70, Ω_M = 0.3, Ω_L = 0.7. Magnitudes are in the AB system (Oke & Gunn 1983).

2 SELECTION OF Z > 9 SOURCES

To make further progress in charting the cosmic star formation history beyond the limits directly accessible to HST, we have selected further candidates for detailed study from three near-infrared surveys: CANDELS (Grogin et al. 2011) and the Frontier Fields (Lotz et al. 2017) from Hubble Space Telescope imaging, and the UltraVISTA ground-based campaign (McCracken et al. 2012).

We applied the following selection criteria to likely z ≃ 9 candidates with m_H-band < 26.5 mag (see Section 3): (i) a 2σ non-detection in all bands bluewards of the Lyman break at the putative redshift; (ii) a 5σ detection in two consecutive bands redwards of the Lyman break; (iii) red Spitzer/IRAC colours (3.6–4.5 |$\mu$|m > 0.5); (iv) a photometric redshift probability distribution permitting z ≥ 9 at 1σ. We used publicly available catalogues for the CANDELS and UltraVISTA fields. AstroDeep catalogues were used for four Frontier Field clusters (A2744, MACS0416, MACS0717, and MACS1149) and we constructed our own catalogues for the clusters A370 and AS1063 following the method described in Laporte et al. (2016). The Spitzer/IRAC data used in this paper are drawn from several surveys and therefore not uniform in depth. CANDELS IRAC data have a 5σ depth of 25.5 mag; the 5σ depth of the UltraVISTA IRAC data is 25.8 mag and the Frontier Fields data reach a 5σ depth of 26.5 mag. We extracted the Spitzer photometry using a 1.2 arcsec radius aperture centred at the position of the source on the F160W image, and applied the usual aperture correction. However noting that, for blended sources, this may overestimate the IRAC flux, in those cases we extracted the photometry more carefully, taking into account the shape of the galaxies with GALFIT (Peng et al. 2010) assuming a Sersic profile.

Including MACS1149-JD1, these criteria led to six suitable sources distributed as follows: two (presumably lensed) sources viewed behind a Frontier Fields cluster, three in the CANDELS fields GOODS-N and GOODS-S, and one in UltraVISTA survey (see Table 1). The most restrictive criterion in their selection is that relating to the red IRAC colour as this can only be applied to relatively bright targets given the depth difference between ∼1 |$\mu$|m data and the IRAC images. All these sources have already been identified as promising high-redshift candidates by previous studies (Ishigaki et al. 2018; Oesch et al. 2018; Bowler et al. 2020).

3 PHOTOMETRIC ANALYSIS

The spectral energy distributions (SEDs) and photometric redshift likelihood distributions for our selected six z > 9 candidates are presented in Fig. 1. The SEDs utilize photometry drawn from publicly available imaging data taken with Hubble and Spitzer Space Telescopes and, where available, the VISTA deep surveys. We matched the Point-Spread Function (PSF) of HST images using Tiny Tim models (Krist, Hook & Stoehr 2011). We extracted the total SED following the method described in Finkelstein et al. (2013) with an aperture correction computed from the Sextractor MAG_AUTO measured on the F160W image from WFC3/IR or the H-band from VISTA. The IRAC/Spitzer photometry following method described in the previous section . Two sources in our sample, MACS0416-JD and MACS1149-JD1, are located behind Frontier Fields clusters and thus are gravitationally lensed. We estimated the magnification for these using the online magnification calculator for the Frontier Field survey¹ which is based on lensing models developed by Hoag et al. (2016), Caminha et al. (2017), Jauzac et al. (2014), Richard et al. (2014), Kawamata et al. (2016, 2018), Ishigaki et al. (2015), Johnson et al. (2014), and Grillo et al. (2015). This yields a magnification of 1.90|$^{+0.03}_{-0.02}$| and 28|$^{+44}_{-11}$|⁠, respectively for MACS0416-JD and MACS1149-JD1. The latter estimate represents a revision of the magnification adopted by Hashimoto et al. (2018).

Figure 1.

Spectral energy distributions and best-fitting BAGPIPES (Carnall et al. 2018) models for the six galaxies in Table 1. The blue points represent photometric data or limits, and the red points represent the predicted fluxes. All models incorporate contributions from nebular emission (see the text). The right-hand inset panel shows the redshift probability distribution, while the left-hand inset panel shows the posterior distribution of the stellar mass as a function of the age of the stellar population. HST upper limits are showed at 1σ, IRAC limits are at 2σ.

We simultaneously determine the optimal photometric redshift and a range of physical properties for the six selected galaxies using BAGPIPES (Carnall et al. 2018). We allow the photometric redshift to lie in the range 0.0 < z_phot < 10.0 which, for z ≤ 9.0, permits the IRAC excess to be formed, in part, from [O iii] line emission (Roberts-Borsani et al. 2020). BAGPIPES generates H ii regions from CLOUDY photoionization code (Ferland et al. 2017) and assumes that the nebular emission is the sum of emission from H ii regions of different ages following Byler et al. (2017). In particular, it is capable of reproducing intense [O iii] emission with equivalent widths EW[O iii + H β] >1000 Å, sufficient for consideration of dominant contributions to an IRAC excess at high redshift. BAGPIPES assumes that the gas-phase metallicity is similar to that of the stars producing the ionizing radiation (with a range of Z ∈[0.0:2.5] Z_⊙). Finally, BAGPIPES accounts for dust emission via the Draine & Li (2007) model, with a reddening law following the Calzetti (2001) relation. We fit the SED of our candidates assuming a single-component model with four possible star-formation histories (SFH): delayed, exponential, constant, or burst-like. The best SED-fit is selected to be that with the largest likelihood and is presented as the grey curve for each candidate in Fig. 1, with the associated physical parameters listed in Table 1. The best fit is always obtained for the delayed or constant SFH. Although we tabulate the age of the best-fitting stellar population, a more meaningful measure is the fraction (expressed as a percentage) of the observed stellar mass that was assembled prior to a redshift z = 10. For the two gravitationally lensed galaxies, the luminosities and stellar masses have been corrected for the magnifications assumed above. Table 2 shows the range of parameters adopted in BAGPIPES for each particular SFH.

Table 1.

Bright z_phot ≥9 galaxies with red IRAC colours (3.6–4.5 |$\mu$|m > 0.5) with physical properties determined using BAGPIPES (Carnall et al. 2018). Photometric redshift uncertainties represent 1σ values and absolute magnitudes and stellar masses are corrected for lensing magnification where appropriate (see the text). The final column displays the inferred percentage of the presently observed stellar mass formed before z = 10.

Target	RA	DEC	zphot	|$\rm {\mathit{ M}_{uv}}$|	|$\rm {\mathit{ M}_{\star }}(\times 10^9$|M⊙)	Age (Myr)	fM(z> 10) (per cent)
MACS0416-JDa	04:16:11.52	−24:04:54.0	9.25|$^{+0.08}_{-0.09}$|	−20.83 ± 0.22	1.50|$^{+0.84}_{-0.61}$|	360|$^{+108}_{-157}$|	74.1|$^{+8.1}_{-25.6}$|
MACS1149-JD1b	11:49:33.59	22:24:45.76	9.44|$^{+0.02}_{-0.03}$|	−19.17 ± 0.04	0.44|$^{+0.05}_{-0.04}$|	484|$^{+17}_{-36}$|	86.5|$^{+2.0}_{-4.3}$|
			9.11	−19.12 ± 0.04	0.66|$^{+0.09}_{-0.04}$|	512|$^{+10}_{-12}$|	93.3|$^{+4.2}_{-5.1}$|
GN-z10-3c	12:36:04.09	+62:14:29.6	9.57|$^{+0.23}_{-0.27}$|	−20.72 ± 0.12	1.57|$^{+0.81}_{-0.63}$|	265|$^{+153}_{-145}$|	38.0|$^{+23.1}_{-38.0}$|
GN-z9-1d	12:36:52.25	+62:18:42.4	9.22|$^{+0.32}_{-0.33}$|	−20.81 ± 0.18	2.18|$^{+1.24}_{-0.85}$|	323|$^{+137}_{-168}$|	65.1|$^{+18.7}_{-32.2}$|
GS-z9-1d	03:32:32.05	−27:50:41.7	9.26|$^{+0.41}_{-0.42}$|	−20.38 ± 0.20	2.47|$^{+1.62}_{-1.03}$|	326|$^{+128}_{-176}$|	70.6|$^{+12.5}_{-25.4}$|
UVISTA-1212e	10:02:31.81	02:31:17.10	8.88|$^{+0.27}_{-0.46}$|	−22.93 ± 0.20	9.7|$^{+5.10}_{-4.72}$|	280|$^{+172}_{-174}$|	32.4|$^{+47.2}_{-32.4}$|

Note.^a Laporte et al. (2016), ^b Zheng et al. (2012), ^cOesch et al. (2014), ^dOesch et al. (2018), ^eBowler et al. (2020)

Table 2.

The BAGPIPES parameter ranges used to fit the SEDs of the six targets studied in this paper. The age (or maximum age) is chosen to be less than 1 Gyr given the selection criteria applied to build our sample. The ionization parameter range we used is large to account for extreme values (log U>−2) as observed in previous high-redshift galaxies (Stark et al. 2017).

Parameters	Burst	Exponential	Constant	Delayed
tSF begin (Gyr)	0.0–0.1	0.0 ; 1.0
tSF end (Gyr)	–	–	0.0 ; 1.0	–
τ (Gyr)	–	0.1 ; 10.0	–	0.1 ; 10.0
log M⋆ (M⊙)	0.0 ; 12.0
Z [Z⊙]	0.0 ; 2.5
Av (mag)	0.0 ; 3.0
log U	−4.0 ; 0.0

Despite the different HST and ground-based data sets from which our six IRAC excess sources were selected, they span a limited range in luminosity and stellar mass. As a consequence of its detection to the shallower magnitude limit with a ground-based telescope, UVISTA-1212 is the most luminous and massive galaxy in the sample. Later in Section 5, we discuss how representative is our sample with respect to the larger population of sources thought to be at z ≃ 9. For all of the sources except UVISTA-1212, the most probable photometric redshift derived is above z ≃ 9 and thus a Balmer break may be the most likely explanation for the IRAC excess. None the less, as discussed by Hashimoto et al. (2018) and Roberts-Borsani et al. (2020), alternative explanations based on contributions from strong nebular emission lines or dust reddening can be considered.

Although [O iii] 5007 Å is redshifted out of the IRAC 4.5 |$\mu$|m band for z > 9.1, Hashimoto et al. (2018) considered the contribution of H β and [O iii] 4959 Å. Even for a very young starburst with a hard radiation field (log U = −2.0) and low metallicity (20 per cent solar), it is not possible to reproduce a significant IRAC excess. Furthermore, for the youngest ages, the IRAC 3.6 |$\mu$|m flux is raised by the presence of a significant nebular continuum. Likewise, [O ii] 3727 Å will enter the 4.5 |$\mu$|m band for redshifts above z ≃ 9.7. The photometric likelihood distribution does extend to such a high redshift for some of our sources (Fig. 1). However, the IRAC excess would only be matched if [O ii] had a rest-frame equivalent width of over 1000 Å which seems improbable. Finally, it is possible to match the SEDs by combining a delayed SFH with dust extinction. In this case, to match the IRAC excess, the extinction would typically need to exceed A_V >1.7 mag, which is larger than that inferred in previous measures of reionization era galaxies (e.g. Schaerer & de Barros 2010; Strait et al. 2021).

In Fig. 1, we show the posterior distribution of the age and stellar mass to indicate the reliability of concluding there was star formation prior to z = 10 as listed in Table 1. Despite uncertainties in the available photometry, the shape of the posterior is reasonably well-defined for ages greater than 200 Myr, indicating the likelihood of significant earlier star formation. The BAGPIPES solution for MACS1149-JD1 can be compared with the solutions derived by Hashimoto et al. (2018) and Roberts-Borsani et al. (2020), both of which included ALMA [O iii] 88 |$\mu$|m fluxes and upper limits on the dust continuum in Band 7 as additional constraints. Although Roberts-Borsani et al advocated the use of ALMA to break the degeneracy between optical [O iii] emission and a Balmer break as contributors to the IRAC excess, it is not a key issue in this case given the spectroscopic redshift of z = 9.11 precludes the former case. Although an academic exercise given its redshift is known, both our photometric redshift (⁠|$z_{\text{phot}}=9.44^{+0.02}_{-0.03}$|⁠) and age (⁠|$484 ^{+17}_{-36}$| Myr) are consistent with the range of early estimates. We include a separate entry in Table 1 for the case where the redshift is fixed at its spectroscopic value of z = 9.11 (Hashimoto et al. 2018).

We emphasize that the inferred contribution of early star formation is uncertain in many cases, particularly for GN-z10-3 whose IRAC excess is the least prominent. None the less, overall our sample is consistent with having significant star formation prior to z ≃ 10, as was concluded for MACS1149-JD1. We will return to a further discussion of this in Section 6.

4 SPECTROSCOPIC FOLLOW-UP

Although Lyman-α is the most prominent rest-frame UV emission line in star-forming galaxies at intermediate redshift (Stark et al. 2010; Stark, Ellis & Ouchi 2011), its visibility deep in the reionization era at z ≃ 9 is expected to be significantly diminished via resonant scattering by neutral hydrogen. For this reason, various efforts have been made to target metallic UV lines such as [C iii]1909 Å and C IV1548 Å with limited success (Laporte et al. 2017; Mainali et al. 2020). Surprisingly, for 11 galaxies with spectroscopic redshifts above z = 7.5, Lyman-α was detected in 7 with additional UV metallic lines in only four cases. In all cases, detecting the intrinsically fainter metallic lines has been more challenging than searching for Lyman-α.

Far infrared emission lines, such as [O iii] 88 |$\mu$|m arising from ionized gas and [C ii] 158 |$\mu$|m from the neutral gas, have been pivotal in spectroscopic confirmation at high redshifts with the ALMA interferometer (Inoue et al. 2014). For the same 11 z > 7.5 galaxies of which seven are accessible with ALMA, [O iii], and/or [C ii] emission have been detected in five cases. Moreover, in cases where both ALMA and near-infrared spectrographs on ground-based telescopes were both used in similar exposure times (e.g. Laporte et al. 2017; Hashimoto et al. 2018), the ALMA detections were more significant. For this reason, we adopted a multifacility approach in attempting to secure spectroscopic redshifts for the targets introduced in Table 1.

We discuss the results of our spectroscopic campaigns for each target in turn below. A summary of the results of our spectroscopic campaign is given in Table 2.

4.1 MACS0416-JD

We obtained 12.6 h (8 h on source) ALMA band 7 observation in Cycle 6 (ID: 2019.1.00061.S, PI: R. Ellis) to search for [O iii] 88 |$\mu$|m in MACS0416-JD. The spectral windows were setup to cover the 1-σ redshift probability distribution estimated from SED fitting (Fig. 1, viz. |$\nu \, \in$| [321.8 GHz : 325.7 GHz], ∪ [328.6 GHz : 335.6 GHz ], and ∪ [340.6 GHz : 347.6 GHz]. The data were obtained between 2019 October and 2020 January and reduced using the ALMA pipeline (v. 5.6.1–8) with a natural weighting scheme, a uv-tapper of 0.35 arcsec with a channel width of 16 MHz and a beam size of 1.00 × 0.76 arcsec in order to maximize the signal-to-noise ratio of any emission line (Tamura et al. 2019). At the HST position of MACS0416-JD, we identified a clear feature at ν = 329.69 GHz corresponding to a redshift z = 9.28 (Fig. 2).

Figure 2.

Top: ALMA contours overplotted on the HST F160W image (left-hand panel) and ALMA band 7 data (right-hand panel) from a 2σ level upwards for MACS0416-JD. The beam size is indicated at the bottom right-hand panel of the ALMA panel. Centre: Extracted 1D spectrum over the full frequency range sampled showing a clear [O iii] 88 |$\mu$|m emission line at 329.69 GHz corresponding to a redshift of z = 9.28. Bottom: Distribution of the pixel signal/noise ratio (SNR) along the line of sight at the position of MACS0416-JD. In the case of non-detection, the histogram shape should be consistent with a Gaussian. The yellow region highlights a deviation from a Gaussian shape at SNR≃6.

The offset between the centroid of the UV continuum and the peak of [O iii]88 |$\mu$|m emission is 0.40 arcsec, corresponding to an observed physical separation of ≃1.8 kpc at z ∼ 9.28. Such a small displacement is consistent with previous [O iii]88 |$\mu$|m observations at high-z (e.g. Carniani et al. 2017; Harikane et al. 2020). We measured the integrated flux on a map created using the immoments casa recipe. We determine f_int = 0.298 ± 0.051 Jy km s⁻¹ which is within the range values found previously in z ≥ 8 galaxies (Hashimoto et al. 2018; Tamura et al. 2019). Assuming a Gaussian profile, we measure a full width at half-maximum (FWHM) of 313 ± 22 km s⁻¹. The significance of the line can be estimated via a histogram of the pixel-by-pixel signal-to-noise ratio ratio (SNR) at the spatial position of the target. In the absence of emission, this will approximate a Gaussian distribution. If emission is present, a second peak will be seen indicating its significance. The lower panel of Fig. 2 reveals such a second peak at SNR ≃ 6 consistent with the analysis above. No dust continuum was detected to a 1-σ level of 18 |$\mu$|Jy beam⁻¹.

We also observed this target in service mode with X-Shooter/VLT between 2019 August and September (ID: 0104.A-0028(A), P.I. : R. Ellis). The target was centred using a small blind offset (<20 arcsec) in the 0.9 arcsec JH slit to reduce the background and a nod length of 4 arcsec was used. The total on source exposure time in the near-infrared arm was 10 h in excellent seeing condition (<0.8 arcsec). The data were reduced using the latest version of the X-shooter pipeline (3.5.0) and inspected visually by two authors (NL and RAM). No emission line was found in any of the three X-Shooter arms. Given the ALMA-based redshift, the absence of Lyman-α implies a 1σ upper limit of 1.8 × 10⁻¹⁸ cgs assuming an FWHM = 200 km s⁻¹ (Fig. 3).

Figure 3.

Portion of the X-shooter spectra of MACS0416-JD indicating no convincing detection of Lyman-α at the redshift of [O iii] 88 |$\mu$|m (red arrow). The grey rectangles show the position of sky lines. The grey line shows the 1σ noise spectrum and blue line shows the extracted spectrum at the position of MACS0416-JD.

4.2 GN-z10-3 and GN-z9-1

We observed GN-z10-3 with MOSFIRE/Keck undertaking a 5 h 20 min on source exposure in 2019 April (Keck 2019A_U128, PI: B.E. Robertson). In order to monitor the seeing and transparency, two stars were included in the MOSFIRE mask. For the final analysis, only frames with seeing less than 0.8 arcsec FWHM were selected representing a total exposure of 4.9 h. Data reduction was undertaken using the 2018 version of the MOSFIRE DRP. As the target was centred in the MOSFIRE slit, emission lines will feature a positive signal accompanied by two negative counterparts separated by the nodding scale of 1.25 arcsec. We identified a convincing (>5σ) emission line at λ = 11891.4 ± 1.3 Å with an FWHM = 133 ± 36 km s⁻¹ and an integrated flux of (1.37 ± 0.24) × 10⁻¹⁸ cgs (Fig. 4). No other line is present in the entire J-band spectra. Assuming this line is Lyman-α, it corresponds to a redshift of z = 8.78. In order to verify its reliability, we examined data comprising independent halves of the total exposure time. The line was retrieved with suitably reduced signal/noise on both spectra. The line width is comparable with that measured for other z > 7 Lyman-α detections with MOSFIRE (e.g. Finkelstein et al. 2013; Zitrin et al. 2015; Hoag et al. 2017).

Figure 4.

Top: 2D MOSFIRE spectrum of GN-z10-3. The red circles show the location of a 5σ emission line detected at λ = 11891.4 corresponding to Lyman-α at z = 8.78 in a pattern consistent with the telescope nodding. Bottom: 1D extracted spectrum (blue) with the emission line highlighted in yellow and the 1σ noise level in grey. The grey rectangles indicate the position of sky lines.

The spectroscopic data refines the redshift value within the photometric redshift likelihood distribution (Fig. 1). Although the IRAC 4.5 |$\mu$|m excess at this improved redshift below z = 9 implies a likely contribution from [O iii] 5007 Å emission, we will return to the analysis of its SED in Section 5. Further sub-mm observations are currently on-going with NOEMA to detect [O i] 145 |$\mu$|m, one of the brightest FIR emission lines (Katz et al. 2019).

During the same observing run, a separate mask was designed for GN-z9-1. A total on source exposure time of 5.6 h was secured in good seeing conditions. The data wer ereduced and analysed in the manner described above. However, in this case, no clear (>5σ) emission line over the entire J-band. Reduced the significance threshold to 3σ, a faint feature is present at λ = 13241.2 Å with a positive signal at the expected position of the target in the slit and two negative counterparts on each side. We estimated a flux of (2.2 ± 0.7) × 10⁻¹⁸cgs. Assuming this line is Lyman-α, the redshift of GN-z9-1 would be z = 9.89 (Fig. 5). However, although consistent with the photometric redshift likelihood distribution (Fig. 1), we consider this a marginal claim. Deeper spectroscopic observation, such as those that will become feasible with JWST, are needed to confirm this redshift.

Figure 5.

Top: 2D MOSFIRE spectrum of GN-z9-1. Illustrated features follow those in Fig. 4. Bottom: The 1D extracted spectrum in a 1.0 arcsec aperture is shown in blue. The grey line displays the 1σ error. A tentative (3σ) emission line is highlighted in yellow which may be Lyman-α at z = 9.89.

4.3 UVISTA-1212

UVISTA-1212, one of the brightest z ≥ 8 candidates known to date, was observed with FLAMINGOS-2 on the Gemini-South telescope via the Fast Turnaround programme (ID: GS-2020A-FT-102, PI: Roberts-Borsani). Observations were accomplished in service mode in 2020 February in J-band using the R3K grism and a 4 pixel wide slit. From an allocation of 9.8 h, a 8.8 h on source exposure was taken in good seeing conditions.

The data were reduced using the FLAMINGOS-2 Python cookbook for longslit observations (specifically the reduce_ls script, which primarily uses the PyRaf module as a wrapper around IRAF functions to reduce the data). More specifically, the nsflat, nsreduce, and nscombine were used to reduce the dark, flat, and arc frames, as well as the standard star and science observations. The final 1D spectrum is shown in Fig. 6. After a careful inspection of the reduced spectra, no clear emission line was found over the entire J-band spectrum.

Figure 6.

The portion of the FLAMINGOS2 spectrum where Lyman-α may lie according to the photometric redshift uncertainties for UVISTA-1212. No convincing emission line is found. The grey rectangles show the position of sky lines.

4.4 GS-z9-1

GS-z9-1 has a photometric redshift of 9.26|$^{+0.41}_{-0.42}$| and we have obtained 9.8 h of observing time (8.0 h on source) from an allocation of 13.5 h with X-Shooter/VLT (ID: 106.20ZH.001, PI: R. Ellis). Observations were undertaken in service mode in good seeing conditions (<1.1 arcsec) between 2020 October and 2021 January.² We carefully checked the centring of the source in the X-Shooter slit and applied the same data reduction procedure as described in Section 4.1. Unfortunately no clear emission line is seen after a first inspection of the data (Fig. 7). At the expected position of Ly α at z_phot, we measure a 1σ upper limit of 7 × 10⁻¹⁹ cgs.

Figure 7.

The portion of the X-Shooter spectrum where Lyman-α may lie according to the photometric redshift uncertainties for GS-z9-1. No convincing emission line is found. The blue line shows the 1D extracted spectrum at the position of the candidate, grey regions show the 1σ noise extracted in the same sized aperture, and light blue rectangles show the position of sky lines. The red arrow displays the position of Ly α at z_phot = 9.26.

5 PHYSICAL PROPERTIES

5.1 Updating the properties

We now reconsider the physical properties of our sample, taking into account, where available, the newly acquired spectroscopic redshifts (Table 3). In addition to MACS1149-JD1, we now have additional redshifts for MAC0416-JD and GN-z10-3, representing half of our original sample. Although we have marginal evidence that GN-z9-1 lies beyond z ≃ 9, for it, GS-z9-1 and UVISTA-1212, we will continue to assume the photometric redshifts given in Table 1. A re-analysis of the SED is particularly relevant for GN-z10-3 whose spectroscopic redshift is now confirmed to lie below z = 9, implying a possible significant contribution to the IRAC excess from rest-frame optical emission lines. The updated results based on BAGPIPES fits to the SEDs, now constrained with a fixed spectroscopic redshift, are given in Table 4. In the three cases where the spectroscopic redshift is confirmed to be above z ≃ 9, the derived properties are naturally similar to those given in Table 1. With the exception of UVISTA-1212, which appears to be an outlier in our sample, all galaxies have a stellar mass ranging from ∼ 1 to 4 × 10⁹ M_⊙, a star formation rate ranging from ∼10 to 30 M_⊙ yr⁻¹, and a low dust attenuation (<0.5 mag). However, for GN-z10-3, despite a spectroscopic constraint of z = 8.78, a re-analysis of the SED incorporating this value may still be consistent with a mature stellar population. This may appear surprising given the IRAC excess could be readily explained with strong [O iii] 5007 Å emission.

Table 3.

Spectroscopic observations of z ≥ 9 candidates. For the [O iii] 88 |$\mu$|m line, we report the peak flux.

Target	Redshift	Instruments	Feature	Flux / 1σ
MACS1149-JD1⋆	9.11	X-SHOOTER	Ly α	(4.3 ± 1.1) × 10−18 cgs
MACS0416-JD	9.28	X-SHOOTER	Ly α	<1.8 × 10−18 cgs
GN-z10-3	8.78	MOSFIRE	Ly α	(1.37 ± 0.24) × 10−18 cgs
GN-z9-1	9.89?	MOSFIRE	Ly α	(2.2 ± 0.7) × 10−18 cgs
U-1212	–	FLAMINGOS2	Ly α	<4 × 10−18 cgs
GS-z9-1	–	X-SHOOTER	Ly α	< 7 × 10−19 cgs
MACS0416-JD	9.28	ALMA	[O iii] 88 |$\mu$|m	0.65 ± 0.13 mJy beam−1
MACS1149-JD1	9.11	ALMA	[O iii] 88 |$\mu$|m	0.84 ± 0.11 mJy beam−1

Note.^⋆ from Hashimoto et al. (2018)

Table 4.

Updated properties for those sources in Table 1 with a spectroscopic redshift. The final column displays the percentage of the presently observed stellar mass already formed before z =10.

Target	zspec	|$\rm {M_{\star }}(\times 10^9$|M⊙)	Age (Myr)	fM(z>10) (per cent)
MACS0416-JD	9.28	1.37|$^{+0.75}_{-0.55}$|	351|$^{+115}_{-153}$|	71.8|$^{+10.4}_{-20.3}$|
MACS1149-JD1	9.11	0.66|$^{+0.09}_{-0.04}$|	512|$^{+10}_{-12}$|	93.3|$^{+4.2}_{-5.1}$|
GN-z10-3	8.78	1.27|$^{+0.70}_{-0.61}$|	275|$^{+184}_{-175}$|	39.1|$^{+24.6}_{-38.9}$|

Fig. 8 shows a re-analysis of the SED shown in Fig. 1 with the redshift appropriately constrained at z = 8.78. The posterior distribution in age and stellar masses is now bimodal and thus it is not possible to distinguish between a young population whose [O iii] 5007 Å emission contributes to the IRAC excess, and an older solution with a contribution from a Balmer break. Similar interpretational ambiguities were seen for z < 9 sources with IRAC excesses by Roberts-Borsani et al. (2020) who argued this degeneracy may be broken with independent measures of [O iii] 88 |$\mu$|m emission. In the case of the young solution for GN-z10-3, the required equivalent width (EW) of [O iii] + H β is ∼1000 Å which is consistent with the range observed in some high redshift galaxies (e.g. Bowler et al. 2020; Strait et al. 2021). The flat continuum between F160W and 360 |$\mu$|m would then imply that [O ii] has an EW >50 Å. Recent studies suggest this is close to the maximum possible for a young starburst (e.g. Yang et al. 2017). Although we are unable to distinguish between these two fits with the current data, we adopt the properties of an older stellar population in Table 4.

Figure 8.

Best fit to the spectral energy distribution of GN-z10-3 now constrained at the spectroscopic redshift of z = 8.78. The blue points and limits represent the observed photometry and red crosses show the expected measurements. The posterior distribution of the age and stellar mass is now bimodal, indicating it is not possible to distinguish the relative contributions of [O iii] emission and a Balmer break to the IRAC excess (see the text).

5.2 The z ∼ 9 Lyman-α fraction

We now consider the question of the visibility of Lyman-α deep in the reionization era. This can be done in terms of a Lyman-α emission fraction for our sample. Thus we have two targets with a secure Lyman-α detection and one (GN-z9-1) which is tentative (see Section 4). All candidates have reliable photometric redshifts at z ∼ 9 or a far-infrared emission line, suggesting the likelihood of contamination of the sample by low redshift interlopers is small. Despite obvious uncertainties for such a small sample, the Lyman-α emission fraction calculated for the sample of |$6\, z\sim 9$| galaxies presented in this paper is very high: |$\chi _{\rm {LAE}} = 0.33_{-0.21}^{+0.28}$| if we include GN-z9-1 or |$\chi _{\rm {LAE}} = 0.50^{+0.26}_{-0.26}$| otherwise.

This high fraction of detections continues a trend of recent results (Roberts-Borsani et al. 2016; Stark et al. 2017; Mason et al. 2018) that have weakened the originally claimed rapid decline above z ≃ 7 (Pentericci et al. 2011; Schenker et al. 2012, 2014). Our results are compared with these earlier trends in Fig. 9. As with Roberts-Borsani et al. (2016) and Stark et al. (2017), we explicitly selected bright galaxies amenable to ground-based spectroscopy with red Spitzer/IRAC colours (3.6–4.5 |$\mu$|m > 0.5, Section 2). However, while the 100 per cent visibility of Lyman-α in Roberts-Borsani et al. (2016)’s sample might be attributed to a correspondingly strong rest-frame optical [O iii] emission responsible for the IRAC excess, for the three z > 9 targets in the present sample, only a Balmer break synonymous with a mature stellar population (e.g. Hashimoto et al. 2018) could be responsible. This suggests that the enhanced Lyman-α fraction arises from the fact that some luminous galaxies, regardless of their [O iii] emission, are capable of producing early ionized bubbles. As discussed by Mason et al. (2018), much of the early work probing the redshift-dependent fraction necessarily targeted sub-luminous galaxies in order to efficiently exploit multislit spectrographs with limited fields of view.

Figure 9.

Fraction of Lyman-break galaxies revealing Lyman-α emission with an equivalent width |$\rm {EW}\gt 25$| Å. Full symbols represent the more luminous population with −21.75 < M_UV < −20.25 and the empty symbols indicate galaxies with −20.25 < M_UV < −18.75. The fraction for galaxies in the present sample is indicated by the red symbols.

There is good evidence that intense [O iii] emitters at low and intermediate redshifts have high LyC escape fractions (Izotov et al. 2018; Fletcher et al. 2019; Nakajima et al. 2020) and/or ionizing efficiencies (e.g. Nakajima et al. 2016), consistent with the creation of ionized bubbles at high-redshift. What remains unclear is how z > 9 galaxies with IRAC excesses indicating a mature stellar population can emit enough Lyman continuum photons to ionize large bubbles. The Strömgren radius of the ionized bubble created by a galaxy scales with ∝ (t_emf_{esc, LyC})^1/3 (e.g. Cen & Haiman 2000), where f_{esc, LyC} is the escape fraction of Lyman continuum photons and t_em the time the galaxies has been emitting such photons.

Assuming an average ionizing efficiency and |$f_{\rm {esc,LyC}} \lt 3{{\ \rm per\ cent}}$|⁠, the galaxies in our z ∼ 9 sample, which all sustained star formation over ≳ 100–200 Myr (see Table 1), can create an ionized bubble as large as those of COLA1 (Hu et al. 2016; Matthee & Schaye 2018) or A370p_z1 (Meyer et al. 2021), two actively star forming z > 6 galaxies with very high escape fractions. Such a large ionized bubble may explain why Lyman-α is observed in MACS1149-JD1 blueshifted by ∼−450 km s⁻¹ with respect to [O iii] 88 |$\mu$|m. In summary therefore, the high detected fraction of Lyman-α emission at z ≃ 9 in luminous galaxies may offer further support of their mature stellar populations.

5.3 How representative is our sample?

Early evidence for mature stellar populations at z > 7.5 (Hashimoto et al. 2018; Roberts-Borsani et al. 2020) has been hard to reconcile with numerical simulations of the first galaxies. While Binggeli et al. (2019) conclude that z ∼ 9 Balmer breaks must be very rare, Katz et al. (2019) find in the contrary that they can be formed easily in their simulations, although dust might complicate the interpretation. In light of our results, it is thus important to assess observationally the prevalence of z ∼ 9 galaxies with older stellar populations.

While we cannot compare Balmer-break objects with a larger population of spectroscopically confirmed z ∼ 9 galaxies, we can assess the prevalence of Balmer-break objects in a parent sample of z > 9 photometric candidates. In order to investigate the relative abundances of the two populations, we have carefully selected z > 9 galaxies candidates using the same procedure and magnitude cuts as in Section 2, but removing the IRAC excess criterion. The candidates were visually inspected by two authors (NL, RAM). Due to the absence of filters bluewards of <8500 Å in ULTRAVISTA, we restricted our search to the Frontier Fields and CANDELS. This additional search results in 14 photometric candidates (including the five Balmer-break objects presented above). Balmer-break objects thus appear to represent a significant fraction (⁠|$\simeq 0.36^{+0.17}_{-0.14}$|⁠) of the z > 9 bright galaxy population, in contrast to the results of numerical simulations. Indeed this Balmer-break fraction could represent a lower limit since the denominator comprises galaxies without spectroscopic redshifts: several could be low redshift interlopers or 8 < z < 9 galaxies (given the large photometric redshift errors) whose IRAC excess may not arise from a Balmer break. This is also the case for the nominator (e.g. GS-z9-1), although to a lower extent.

We can compare the rest-frame UV luminosities of this parent sample with those of our Balmer break galaxies. As most of the candidates are drawn from the GOODS-North/South fields (11/14), we restrict the comparison to the CANDELS fields, ensuring we can avoid the need to make corrections for lensing magnification in the cluster fields. We derive the rest-frame UV magnitude of our candidates using the F125W magnitude, the best-fitting photometric redshift and a UV slope α = 2, (f_ν ∝ ν^α) for the K-correction. We compare the rest-frame UV magnitudes for the two samples in Fig. 10. Although the statistics are poor, our objects with red IRAC colours are not notably brighter or fainter than blue IRAC colour objects. A two-sided KS test is inconclusive and the null hypothesis that the two distributions are different cannot be rejected (p = 0.12). Thus, given the modest sample size to date, our z > 9 Balmer break sources remain broadly representative of the larger z > 9 galaxy population. We note that all z > 9 candidates represent the brightest objects existing at that early time due to obvious limited depth of the imaging data. The hypothesis that Balmer-break objects are less frequent in the fainter population is plausible, but can only be tested with JWST.

Figure 10.

A comparison of the rest-frame UV magnitudes for z ∼ 9 candidates in CANDELS selected without regard to their IRAC excesses with those for the three GOODS-S/N sources discussed in this article.

6 DISCUSSION

In this article, we have argued that from a sample of six z ≃ 9 galaxies with Spitzer/IRAC excess for which spectroscopic redshifts have now been secured for three, possibly four, we infer past star formation histories consistent with significant contributions beyond a redshift z = 10.

We can estimate the earlier rest-frame UV luminosity for each of the six galaxies in our sample assuming the best-fitting SFH. As described in Section 3, the best SED fit is obtained for most of the objects with an evolved stellar population. Table 1 shows the best-fitting parameters for all galaxies, updated in Table 4 with the spectroscopic redshift for 3 of them. Based on these SFHs, we compute the rest-frame |$\rm {L_{1500}}$| luminosity as a function of redshift from z ∼15 to the observed redshift taking into account the dust reddening. For those galaxies without spectroscopic redshifts, the photometric value is assumed. Fig. 11 shows the evolution of |$\rm {m_{1500}}$| as a function of redshift and demonstrates that each galaxy is as luminous in its earlier life. While these past star formation histories are clearly illustrative, they serve to demonstrate the feasibility of searches for earlier progenitors with JWST (e.g. Roberts-Borsani et al. 2021).

Figure 11.

Earlier redshift evolution of the rest-frame 1500 Å luminosity computed from our SED fits for the six galaxies in our sample. The filled regions show the 5σ sensitivity of NIRCam/JWST filters (F150W, F200W, and F277W) in 3 h.

From our SFHs, we can also determine the redshift-dependent history of the assembled stellar mass providing independent insight into the evolution of the UV luminosity density within the first 500 Myr (Fig. 12). Two contrasting hypotheses are still debated in the literature: a rapid (Oesch et al. 2018) or a smooth (McLeod, McLure & Dunlop 2016) decline of the UV luminosity density over 8 < z < 11. For the case of a fast decline (∝ (1 + z)⁻¹¹), the luminous galaxies probed would have formed ∼50 per cent of their stellar mass by redshift z =10, whereas for the smooth decline case (∝ (1 + z)⁻⁴) ∼80 per cent of the stellar mass is formed at z = 10. Although our sample is admittedly small, it represents the only relevant data currently available. Averaging the SFHs of our six galaxies, we determine that they formed ∼70 per cent of their mass by z = 10, which seems to favour a smooth decline.

Figure 12.

Redshift evolution of the assembled stellar mass computed from the SFH deduced from our SED fits and averaged over the six galaxies discussed in this paper (red). This is compared with that predicted for two contrasting measures of the rate of decline in the UV luminosity density deduced from large photometric surveys (Oesch et al. 2014; McLeod et al. 2016). We estimate that ≃ 70 per cent of the stellar mass was already in place by z = 10 for the six galaxies in our z ≥ 9 sample.

ACKNOWLEDGEMENTS

NL acknowledges support from the Kavli Foundation. RAM acknowledges support from the ERC Advanced Grant 740246 (Cosmic_Gas). RAM and RSE acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 669253). BER was supported in part by NASA program HST-GO-14747, contract NNG16PJ25C, and grant 80NSSC18K0563, and NSF award 1828315. We thank Adam Carnall for providing an updated version of BAGPIPES allowing high ionization parameter.

Based on observations made with ESO Telescopes at the Paranal Observatory under programme ID 0104.A-0028. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.00061.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. Based on observations obtained at the international Gemini Observatory, a program of NSF’s NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnologia e Innovación (Argentina), Ministério da Ciencia, Tecnologia, Inovaces e Comunicaces (Brazil), and Korea Astronomy, and Space Science Institute (Republic of Korea). The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. ID: 2019

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

https://archive.stsci.edu/prepds/frontier/lensmodels/#magcalc

2

We thank ESO staff for observing during both Christmas and New Year.

Author notes

Hubble Fellow