Early Results from GLASS-JWST. XVII. Building the First Galaxies—Chapter 1. Star Formation Histories for 5 < z < 7 Galaxies

The buildup of stellar mass in a galaxy is expressed as its star formation history (SFH), that is, the rate of star formation over an epoch of cosmic time. Measuring the SFHs of galaxies over cosmic time has been a long-sought goal of extragalactic astronomy, extending back at least as far as the 1950s (Baade 1951). The best and most direct method, measuring the Hertzsprung–Russell diagrams of populations of stars for our own galaxy and its nearest neighbors—spectacularly advanced by the Hubble Space Telescope (HST; e.g., Tosi et al. 1989; Gallart et al. 1996; Tolstoy & Saha 1996; Aparicio & Hidalgo 2009; Cignoni & Tosi 2010)—is unfortunately limited to a number of present-day galaxies (e.g., Dalcanton et al. 2009, 2012) in the local universe.

To look back in time to galaxies at cosmic distances has been critically challenging because unresolved stellar populations cannot constrain SFHs for ages older than 2 Gyr, beyond which the dependence of an integrated spectrum on age and metallicity thwarts any comprehensive analysis (e.g., Worthey 1994). The ability to resolve SFHs over <2 Gyr of "lookback" is due to the rapid evolution of light-dominating A stars over this time period (Dressler & Gunn 1983), resulting in a tight correlation of SFH and spectral energy distribution (SED) that quickly vanishes for ages τ > 2 Gyr (see also Poggianti et al. 1999). The discovery of "late bloomer galaxies"—i.e., galaxies with rising SFRs (rather than falling as conventionally described; Dressler et al. 2018)—made use of this by observing the growth of galaxies from z ∼ 0.75 to 0.45, a significant (2 Gyr) fraction of their lifetimes. The prevailing model for SFHs has been for homologous lognormal forms (Gladders et al. 2013; rapidly rising and slowly declining) over a Hubble time (Diemer et al. 2017), but Dressler et al. (2018) found that ∼25% of Milky Way–mass galaxies more than doubled their mass during this relatively late period in cosmic history. This study made clear that the SFHs of the first billion years of star formation, 5 < z < 20, would greatly benefit from the strong correlation of SED and SFH. Therefore, the potential of the James Webb Space Telescope (JWST) to produce such data has been eagerly awaited (see, e.g., Tacchella et al. 2022, for state-of-the-art HST-based work).

In this paper, we show very early results of this potential through the analysis of 5 < z < 7 SEDs derived from seven-band NIRCam images obtained as part of the GLASS-JWST Early Release Science program (Treu et al. 2022a). We start from the fluxes for ∼6600 galaxies extracted by Merlin et al. (2022, Paper II) in NIRCam parallel observations of NIRISS spectroscopy of the cluster A2744. We here exploit a program written expressly for the purpose of reconstructing the histories of galaxies at 5 < z < 7 from the rich information encoded in their SEDs for the first time on real data. A recent paper by Whitler et al. (2023) applies similar techniques to JWST data of galaxies at z > 8, when the universe was less than 500 Myr old. By focusing on a sample at slightly lower redshift, we have access to redder rest-frame wavelengths and give sufficient time for older stellar populations to emerge.

These first results demonstrate that this technique will be a powerful tool in studying and understanding the buildup of stellar mass—the growth—of the first galaxies in our universe.

The paper is organized as follows. Section 2 and Appendix A explain and summarize the methodology in advance of a forthcoming paper that will provide more detail. Section 3 presents the observations and the selection of the sample. Section 4 shows 24 examples of the SEDs of 5 < z < 7, with a focus on examples of various SFHs that we find—mainly long and short falling SFHs and prominent bursts. We describe the character of these SFHs and comment on the fidelity and robustness of the results, aided by Appendixes B and C, which also investigate the reliability of our identifications. Section 5 briefly contemplates the bounty expected from deeper, wider surveys already on the JWST observing schedule.

A standard cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹ is adopted.

The method we apply to derive SFHs follows the work of Kelson et al. (2016), who developed a maximum-likelihood Python code to deconstruct observed SEDs from the Carnegie-Spitzer-IMACS Survey (Kelson et al. 2014), in terms of best-fit stellar population templates. The program proved particularly effective in isolating the light from younger (<2 Gyr) populations of A stars, confirming an early photometric study (Oemler et al. 2013) that showed that 25% of Milky Way–mass galaxies at z ∼ 0.5 had rising SFRs (rather than falling at z < 1, as conventionally described), producing at least 50% of their stellar mass in the previous 2 Gyr—so-called "late bloomers."

Dressler has followed this approach by writing SEDz*, a Fortran code for analyzing stellar populations for z > 5 galaxies, all of which have stellar ages of <1 Gyr. It was developed to answer a basic question: what do the SFHs of the earliest galaxies look like? That is, its purpose is to quantify SFHs—their forms—for a sample of very young galaxies. Other codes measure SFHs, of course, for example, PROSPECTOR (Johnson et al. 2021), EAZY (Brammer et al. 2008), and BEAGLE (Chevallard & Charlot 2016). They are essential tools in the study of galaxy evolution; they can provide information on many properties and behaviors for galaxies over a wide range of cosmic histories. The SEDz* program was developed to focus on SFHs over a special epoch. It takes advantage of a unique astrophysical opportunity: the light from galaxies in the first billion years is dominated by A stars, and A stars are the best clock for studying young stellar populations (Dressler & Gunn 1983; Couch & Sharples 1987). Because they evolve rapidly over a gigayear, the SFHs of A star–dominated populations can be derived (rather than inferred) from SEDs, and vice versa. There are no complexities; in addition to providing the proper clock, A stars are among the simplest stars (Morgan & Keenan 1973), main-sequence stars with relatively simple atmospheres. Hydrogen absorption dominates the opacity, not metal lines, thereby much reducing the influence of metal abundance and the environments of star formation. These are significant advantages for SEDz* in measuring SFHs during the first billion years of cosmic history.

Operationally, the data inputs to SEDz* are SEDs—flux measurements in NIRCam's wide bands. ¹⁶ The SEDz* program uses a nonnegative least-squares (NNLS) engine (Lawson & Hanson 1995), which combines custom star formation templates based on work by Robertson et al. (2010) and further developed by Stark et al. (2013 ¹⁷ ; see Section 3) and is used by the NIRCam team to simulate observations of galaxies in the GTO deep fields (Williams et al. 2018). The custom templates made for SEDz* cover the redshift range 3 < z < 12. ¹⁸ For each template, there is a starting redshift with fluxes proportionate to stellar mass along with the evolution of this template, that is, the expected fluxes when observing that stellar population at subsequent epochs. The program divides the 3 < z < 12 epoch into integer steps (i.e., 12, 11, 10, ..., 3) and operates with two sets of SED templates; one set is for 10 Myr bursts starting at epochs 5 < z < 12, and the other is for continuous star formation (CSF) starting at redshifts as early as z = 10 but not after z = 4.

The process is to build up the stellar population by combining templates, epoch by epoch. For each step, NNLS works to improve the fit of combined fluxes to the observed SED, with each step evolving the preceding star formation forward and adding the fluxes of subsequent populations, like building a wave. For example, starting with a single 10 Myr burst at z = 12, SEDz* uses NNLS to find the mass that minimizes the χ² ${\tilde{\chi }}^{2}$ of the fit to the observed SED. Next, using the template of the z = 12 burst evolved to z = 11 and adding a population starting at z = 11, NNLS finds the stellar mass for each epoch (as small as zero) that further lowers the ${\tilde{\chi }}^{2}$. The process continues for bursts at 10, 9, ..., 5 and adds potential CSF populations (in units of 1.0 M_⊙ yr⁻¹ over the z-integer-bound epoch) starting at z = 10 and as late as z = 4, again, as needed to improve the fit of the composite stellar populations to the observed SED. Detecting a CSF population defines the epoch of observation (OE) and ends the possibility of later bursts. To summarize, from z = 12 to 4, a sum of bursts and the possible onset of CSF improves the goodness of fit of the model-to-observed SED, measured as reduced ${\tilde{\chi }}^{2}$.

The SEDz* program does not prejudge the SFH. By using a combination of individual bursts and CSF templates, SEDz* has the flexibility to describe rising or declining, as well as bursty, SFHs. Whatever the elements contributing to the best fit to the observed SED, there is no preference or bias for any particular history. The SEDz* program provides error bars for each mass point in the SFH. If they are large enough, multiple histories are possible, but this is not the case here. This means that, for a sample like this one, each SFH generated by SEDz* is unique; it is not meaningful to ask if another SFH would fit as well or better. With SEDz*, the minimum ${\tilde{\chi }}^{2}$ solution defines the SFH; there are no additional indicators or diagnostics of whether this is the "real" SFH or not. ¹⁹

In this first exploration, we neglect the potential impact of dust. We will revisit the issue in future papers, when we will analyze larger samples in detail. However, we point out that most of the mass we find in running SEDz* on the present data is in comparatively old populations (hundreds of megayears), for which dust is not likely. In support of this, we notice that galaxies in our sample seem to be well described by SEDs with minimal dust, even for the youngest populations, and do not resemble the strong attenuation of ultraviolet light that is characteristic of populations with even moderate amounts of dust. Our choice is consistent with the results of several papers that show these initial JWST selected samples are uniformly blue (Nanayakkara et al. 2022, Paper XVI) and fairly dust-free (see Finlator et al. 2011; Jaacks et al. 2018; Mason et al. 2022, for a theoretical perspective), except for a few rare examples (Barrufet et al. 2022).

As discussed below, the main outputs of SEDz* are redshifts and SFHs. ²⁰ To evaluate the performance in determining redshifts, a simulated catalog of galaxies following the distribution of Williams et al. (2018) was prepared in the context of the JWST Data Challenge 2 (DC2) to exercise the GTO pipeline. Simulated NIRCam images of many deep fields were created to test a wide range of analysis tools, with SEDz* among them. Although redshift recovery is not its purpose, SEDz* is competitive with the commonly used photoz codes (based on DC2 results), even to the extent of identifying z < 3 galaxies, something SEDz* was not designed to do and is actually outside its native ability (because stellar populations exceed an age of 1 Gyr).

To test the performance of SEDz* on SFHs, a code that accepts parametric SFHs or generates stochastic SFHs was written and applied to the DC2 sample; SEDz* was found to recover declining SFHs, as well as single and multiple bursts of star formation, with good fidelity but exhibiting a limited dynamic range for SFHs rising from higher redshift, as would be expected (older populations are fading while younger populations ascend). Examples of these simulations are given in Appendix A, along with a description of how SEDz* is used to make these tests that provides further insight into how the program works.

A more detailed description of the features of the code and experience in its development, along with the results of DC2 and further SFH tests, is planned for a forthcoming paper (A. Dressler et al. 2022, in preparation). The SEDz* program has passed its early tests of measuring redshifts and SFHs; the purpose of this paper is to demonstrate its potential by applying it to some of the first real data from GLASS-JWST.

We use the NIRCam data obtained on 2022 June 28–29 in the context of the GLASS-JWST-ERS survey (Treu et al. 2022a). They consist of images centered at R.A. = 35017025 and decl. = −303375436 and taken in seven wide filters: F090W, F115W, F150W, F200W, F277W, F356W, and F444W. The 5σ depths in a 030 diameter aperture are in the range 29–29.5 AB mag. The final images are PSF-matched to the F444W band.

Details of the NIRCam data quality, reduction, and photometric catalog creation can be found in Paper II. As in other papers of this focus issue, we neglect the potential impact of lensing magnification, which is expected to be modest at the location of the NIRCam parallel field. We will revisit lensing magnification in future work. However, we note that the effect is small, and that lensing cannot affect SFH shapes, only their magnitude.

Out of the 6590 detections, we select sources with a signal-to-noise ratio (S/N) >5. This S/N is obtained by averaging the S/N of three bands: F200W, F277W, and F356W. We run SEDz* on this sample and obtain a redshift for all of the sources. From now on, we consider only galaxies with 5 < z < 7, for a total of 123 objects. ²¹ We visually inspect the photometry of these galaxies, removing those with artifacts or chip boundaries affecting them in at least one band. We also select only galaxies with an SEDz* fit (${\tilde{\chi }}^{2}\lt 2$), for a total of 43 galaxies. We further reject 10 galaxies with clear (single-band) photometry errors, poor fits of the SED, or ambiguous ${\tilde{\chi }}^{2}$ (redshift) measurements. From the remaining 33, we select 24 that best illustrate the distribution and diversity of SFHs we find. ²²

In this section, we investigate a sample of SEDs from the JWST-GLASS NIRCam observations. By way of example, Figure 1 shows the SED of galaxy 1423 from the GLASS catalog of Paper II. This galaxy has a redshift of z = 7.5, so it does not enter the sample, but it is otherwise indistinguishable from the 33 galaxies we have chosen and a good example to show the potential of SEDz* on a wide range of observations.

Zoom In Zoom Out Reset image size

The data for galaxy 1423 are shown in three panels, the format for SEDz* results. The middle panel shows the SED, with crosses marking the fluxes with seven NIRCam "wide" filters in nJy, with errors. The SEDz* program has been run as described in the previous section on this SED, yielding the reduced ${\tilde{\chi }}^{2}$ distribution shown in the left panel, with a well-defined minimum at z ≈ 7.5, between the two lowest points. Although SEDz* calculates the SED fit in integer redshifts, there is sufficient precision in the fits to interpolate between adjacent values; the first decimal in the redshifts is significant. Returning to the middle panel, the magenta band that "connects" the points is generated by randomly perturbing the input SED with 1σ flux errors (typically a few nJy) and repeating the NNLS fitting 21 times. The width of the band represents the quartile range of NNLS fits with this relatively low-S/N example. The run of ${\tilde{\chi }}^{2}$ presented in the left panel shows a significant drop between z = 6 and 7, interpolated by the program to derive the z = 7.5 value recorded in the central panel, in reasonable agreement with the estimate of the redshift as 6.9 obtained with the code EAzY (Brammer et al. 2008) with the default V1.3 spectral template and a flat prior in Paper X.

The real payoff for these data and our methodology is the SFH shown in the right panel; the green boxes record the stellar mass added to the galaxy at succeeding epochs. Error bars for these boxes do not include probable systematic errors (for example, the errors associated with photometry at these faint levels), but the relative uniformity of this clearly declining rate of star formation suggests that they are not severely underestimated. As we will see, this slowly falling SFH, starting at z = 12 (and likely earlier), shows a combination of features we find in this early study, both parametric and stochastic. The total mass of 10⁹ M_⊙ is shown at z = 7.5, with the OE marked by the red dot/blue circle, while the squares showing the mass buildup are at their integer redshifts by necessity. The galaxy shows a factor of ∼3 buildup in stellar mass added after its formation at z ∼ 12. This relatively low mass is consistent with the 10–20 nJy fluxes measured for standard mass-to-light ratios (Santini et al. 2022, Paper XI).

We continue by showing the SEDs of 24 galaxies chosen to illustrate the full sample. The details of such galaxies are given in Table 1, while color composite cutouts are presented in Figure 2.

Zoom In Zoom Out Reset image size

Table 1. Summary of the Main Galaxy Properties of the Sample Discussed in the Paper

ID	R.A.	Decl.	z	F444w	$\mathrm{log}({M}_{* })$	Notes
	(J2000)	(J2000)		(μJy)	(M⊙)
1423	00:14:03.29	−30:21:46.7	7.5	0.019 ± 0.003	9.1
733	00:14:01.51	−30:22:12.6	6.5	0.031 ± 0.004	9.5	L
2199	00:14:04.77	−30:21:14.8	6.9	0.046 ± 0.003	9.5	L
3459	00:14:01.05	−30:17:41.7	5.7	0.032 ± 0.004	9.4	L
3665	00:14:00.96	−30:18:00.8	5.2	0.019 ± 0.003	9.0	L
4761	00:13:58.16	−30:19:23.5	6.4	0.094 ± 0.005	9.7	L, (2)
5429	00:13:58.11	−30:19:03.9	7.0	0.023 ± 0.004	9.2	L
5879	00:13:54.26	−30:18:51.8	5.1	0.018 ± 0.003	8.9	L
6155	00:14:02.99	−30:18:47.8	6.8	0.012 ± 0.003	9.2	L
4144	00:14:02.29	−30:18:31.9	5.9	0.024 ± 0.004	9.2	L, *
483	00:14:05.72	−30:22:22.5	6.5	0.036 ± 0.004	9.2	B, *
1663	00:14:04.39	−30:21:38.5	5.6	0.007 ± 0.002	8.3	B, *
2104	00:14:02.85	−30:21:19.3	5.6	0.037 ± 0.003	9.1	B, *
2321	00:14:06.58	−30:21:18.2	5.8	0.014 ± 0.003	9.0	B, *
2626	00:14:04.50	−30:20:58.6	6.2	0.056 ± 0.004	9.1	B, *
2827	00:14:03.49	−30:20:51.0	5.6	0.013 ± 0.003	8.7	B, *
2850	00:14:03.02	−30:20:49.0	5.9	0.026 ± 0.004	9.3	B
2993	00:14:00.58	−30:20:40.8	6.7	0.021 ± 0.003	9.1	B
3824	00:13:53.63	−30:18:09.7	6.6	0.019 ± 0.003	9.0	B
3918	00:14:02.51	−30:18:15.2	5.5	0.016 ± 0.003	8.8	B
4026	00:13:52.80	−30:18:21.8	5.7	0.198 ± 0.008	9.8	B, (2)
4587	00:14:02.65	−30:19:16.2	5.7	0.028 ± 0.004	9.1	B
4736	00:14:02.01	−30:19:22.6	6.3	0.012 ± 0.002	8.6	B
5394	00:14:00.89	−30:19:04.7	6.3	0.076 ± 0.005	9.3	B, *
6324	00:14:01.78	−30:18:13.3	6.6	0.054 ± 0.004	9.4	B
680	00:14:05.63	−30:22:14.9	5.0	0.02 ± 0.003	8.9	S
1105	00:14:01.61	−30:21:58.6	5.4	0.077 ± 0.002	8.6	S
1127	00:14:02.43	−30:21:57.6	5.6	0.009 ± 0.002	8.8	S
1356	00:14:04.14	−30:21:48.7	5.0	0.017 ± 0.003	8.6	S
2006	00:13:59.17	−30:21:23.8	5.7	0.013 ± 0.003	8.8	S
2838	00:14:02.83	−30:20:48.6	6.7	0.13 ± 0.005	9.6	S, (1)
2864	00:14:05.89	−30:20:47.9	6.4	0.062 ± 0.005	9.4	S, (2)
2872	00:14:03.49	−30:20:51.0	5.6	0.013 ± 0.002	8.7	S
4455	00:14:00.00	−30:19:10.3	5.1	0.036 ± 0.006	9.2	S, *

Note. L: galaxies shown in Figure 3. B: galaxies shown in Figure 4. S: galaxies shown in Figure 5. (1) Galaxy also in the sample selected in Paper X; estimated redshifts are consistent within the uncertainties. (2) Galaxies also in the sample of Paper XVI; estimated redshifts are consistent within the uncertainties. Asterisk: galaxies not shown in the figures.

Download table as:  ASCIITypeset image

Figure 3 shows a selection of eight galaxies at 5 < z < 7 that display long, often continuous star formation back to z ∼ 11–12—early times in galaxy building. Unsurprisingly, then, these SFHs are either level or declining, remarkably steady for such an early birth and growth. Their total mass is typically a few billion M_⊙, with only one (4761) approaching 10¹⁰ M_⊙. Though they resemble the smooth, settled SFHs of much older galaxies, these SFHs are not the most common in our sample, accounting for only 27% of our sample.

Zoom In Zoom Out Reset image size

These galaxies with extended SFHs straddle the end of reionization (Fan et al. 2006; Planck Collaboration et al. 2016), a full gigayear since the probable beginning of galaxy formation and growth, and appear to have been forming stars continuously from z ∼ 12. There are indications of abrupt order-of-magnitude changes in stellar mass added (more and less), but mostly they are smooth SFHs.

In the sample of 33 galaxies, there are no extended, rising SFHs. However, this is consistent with the results of our simulations of stochastic SFHs (Appendix A), which show the difficulty of detecting an old, fading population behind the light of an ascending, younger population (e.g., Iyer et al. 2020; Whitler et al. 2023). On the other hand, many of the falling SFHs with shorter runs of contiguous star formation (see Figure 5) could be masking rapidly rising star formation prior to their detection. Larger and deeper samples will be needed to characterize such behavior and, in general, the distribution of SFH types through what will essentially be continuity constraints.

A more thorough inspection of the SEDs themselves in Figure 3 reveals that all have a substantial older stellar population, in agreement with the derived SFHs. Also consistent is the absence of strong emission lines (raising the flux in a wide band) that has been seen in some of the stochastic or burst-dominated cases we look at next.

Figure 4 shows examples of galaxies with what look like more stochastic SFHs—bursty, repeating, and noncontinuous episodes of <200 Myr. The top two cases are a single and perhaps an extended single burst. The third and fifth examples are galaxies with widely separated bursts of star formation. Each formed roughly half of their stellar mass at z = 11–12, then experienced a comparatively quiescent period with a comparable burst of star formation added near OE; possibly these are interrupted versions of the long SFHs of Figure 3. The blue upturn of each SED relative to longer wavelengths is representative of the younger, ongoing star formation added to a much older population, with little in between.

Zoom In Zoom Out Reset image size

The fourth and sixth examples are two nearly identical triple bursts (with a mass difference of a factor of 3) that look like (at least for the two pairs of early bursts) independent events, possibly followed by decaying star formation for each at the following (adjacent) epoch. The strong resemblance of the SEDs and SFHs for these two cases demonstrates the strong SED 〈 − 〉 SFH relation for the 5 < z < 12 galaxies that SEDz* exploits. The last two panels show two SFHs that suggest extended star formation episodes—three and four starbursts spread over a few hundred megayears, continuous or episodic. These might be cases of extended burst histories, but, according to the tests done for Figures 8 and 9 of Appendix B, they could possibly be single early single bursts, z = 9 and 10, that are badly resolved by the SEDz* fitting of the SED. Arguing against this interpretation is star formation at OE in both cases, unlike the typical simulated cases in Figure 8.

Lastly, Figure 5 shows further diversity in high-redshift SFHs. All eight are examples of a "run of star formation" over three adjacent epochs, typically at z = 8 − 7 − 6. As just discussed, a detection by SEDz* of adjacent episodes of star formation at higher redshift observed in a galaxy at much lower redshift (e.g., z = 5 or 6) may only be a failure to resolve a single, larger burst. Unlike those, however, for six of the eight cases we show in Figure 5, star formation is spread over lower redshifts. The epochs of 8, 7, and 6 have durations of 125, 171, and 244 Myr, over which time the SED changes markedly, which means that their SFHs are not easily confused with bursts, for example. Rather, these half-gigayear runs of star formation producing between 5 × 10⁸ and 5 × 10⁹ M_⊙ by redshifts 5–6 could be the beginnings of typical M* galaxies. The lack of rising star formation before these usually declining, short SFHs is interesting, since at these epochs, lower rates at z ∼ 9, for example, should be detectable. A possibility is that the first episode of these short SFHs is a burst lasting only tens of megayears, unresolved in the ∼100 Myr of the first epoch.

Zoom In Zoom Out Reset image size

The presence of both short and long continuous runs of star formation, along with a diversity of starbursts also occurring over short or long timescales, suggests that the starburst examples are undergoing significant merger events, while the continuous SFHs, whether short or long, suggest galaxies left largely undisturbed through their early growth.

We note that, given the heterogeneous nature and varying quality of this early sample, we have not tried to estimate the frequencies of the different types of SFHs identified. However, for what it is worth, the three types of SFHs in Figures 3, 4, and 5, including the cases in the 33 galaxy sample that are not shown, amount to 27%, 45%, and 27%, a roughly uniform distribution, with bursty behavior leading the way.

The principal goal of this paper has been to show how broadband imaging with NIRCam on JWST can provide an excellent measurement of both redshift and SFH; the latter is critical to understanding the assembly of the earliest galaxies. High-redshift SFHs will also provide an additional probe of galaxy populations at z = 11–12 and possibly beyond; the SFH of each lower-redshift galaxy probes earlier times, adding to what we will learn from actually observing galaxies at that time.

Future data will include deeper NIRCam imaging with more multiple, highly dithered images that will minimize cosmetic defects and improve photometric measurements; uniformity in all bands is critically important for SED analysis. Improvements in S/N and defect-free imaging will be crucial for optimizing the ability of SEDz* to delineate and distinguish SFHs for the first galaxies.

This initial small sample appears to exhibit the types of SFHs that astronomers would expect, including falling SFHs ("tau models") leading to the near-constant star formation rate that dominates the epoch of greatest growth, z ∼ 3 → 1. But even this small sample suggests that starbursts (stochastic histories) and as-yet-undetected rising SFHs—expected to dominate as we look further back to "the beginning"—are showing up. If the results of this study are indicative of what is to come, astronomers can look forward to the JWST fulfilling its "prime" mission better and sooner than we could have hoped.

A.D. gratefully acknowledges the support of the NIRCam team for the opportunity to contribute to the NIRCam science program on the earliest galaxies. The DC2 "deep field" simulation was a Herculean task that provided the essential guidance for the program described here and the NIRCam GTO program to come. A.D. also thanks the GLASS team for allowing him in, to lift the curtain on the growth of the first galaxies—what the "HST & Beyond" Committee (Dressler et al. 1996) dared to dream.

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program JWST-ERS-1324. The JWST data used in this paper can be found on MAST at DOI: 10.17909/fqaq-p393. We acknowledge financial support from NASA through grant JWST-ERS-1324. D.S., K.M., M.R., and A.D. are supported by JWST/NIRCam contract to the University of Arizona NAS 5-02015. A.M. acknowledges financial support through grant PRIN-MIUR 2017WSCC32.

In this Appendix, we show examples of simulations of SFHs that have been used to test the robustness of the resulting SFH. We first show how the basic types of SFHs we find can be validated by the simulation process by examining six "generated" examples. In Appendix B, we address the uniqueness of the SEDz* solutions, specifically, the robustness of the long SFHs of Figure 3 and the possibility that some or all of these could actually be bursty SFHs misidentified by SEDz*.

The SEDz* program solves for the buildup in stellar mass per epoch by finding a maximum-likelihood solution of the observed SED from combinations of stellar population templates. The process is reversible because of the strong covariance of SFH and SED over the first gigayear of cosmic history. Our tests of SEDz* consist of "composing" SFHs, like the green squares in Figures 3, 4, and 5, over 5 < z < 12 and then using SEDz* "backward" to construct the SED that generates the invented mass profile. This has been done for parametric SFHs, such as smooth rising or falling SFHs over many epochs as well as shorter connected epochs of star formation or bursts. The three steps are (1) invent SFHs like those observed in the figures, (2) use the SEDz* code in reverse to produce the SED that would generate the simulated SFH, and (3) use SEDz* to generate 21 SFHs by perturbing the SED by expected random errors, including error bars on the green boxes to represent the quartile ranges of those 21 solutions. Comparing the simulated and recovered SFH tests the fidelity of the code.

Figure 6 shows six examples. The dark blue dots are the "invented" run of stellar mass, and the small red dots are the integrated mass; the green dot with a red ring indicates the total accumulated mass. In all cases, the simulated mass buildup is well traced by the recovered mass growth, even as the epoch-to-epoch agreement bounces around, as would be seen in the evolution of this diagram stepped through earlier epochs than those shown in this Appendix and Figure 6 (we show examples of this epoch-by-epoch buildup in Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The top three histories were constructed to resemble those of Figure 3: long, falling, mostly continuous runs of star formation since z = 11 or 12. We show, top left to right, CSF until z = 6 or 5. In these cases, the star formation is added in 10 Myr bursts only; there is no constant star formation at OE. ²³ The example on the right was designed to be smooth and declining with a stochastic "noise" of a factor of ∼2, similar to the steeper examples in Figure 3. The middle panel of Figure 3 shows a slightly rising SFH that is recovered as slightly falling, and for the right example, the modestly rising SFH that is simulated is recovered as flat. In summary, the basic character of a long stretch of star formation is recovered, but the slopes are not accurate. This is unimportant for our purposes here, but it suggests that it might be necessary to see what systematic error produces the effect and possibly adjust for it.

The qualitative agreement is good for these top-row examples; the smoothness and falling character are recovered, although in each case, the SFH is higher earlier and lower later than what was simulated. The typical scatter of a factor of 2 comparing model (purple points) and recovered masses (green squares) is remarkably good and certainly sufficient for a basic description of the SFH—the goal of this program.

The bottom panels are simulations of different degrees of bursty behavior. The bottom left case simulates three independent bursts, 10⁹ M_⊙ at z = 12, 3 × 10⁹ M_⊙ at z = 9, and 6 × 10⁸ M_⊙ at z = 7. The agreement is good, except for the alias mass added at z = 10. These simulations suggest that SFHs dominated by bursts, like those in Figure 4, are well recovered by SEDz*. Again, the displacement of the points from the model by factors of 2 is insignificant our purpose.

The example in the middle panel is for the kinds of shorter but continuous runs of star formation shown in Figure 5. Again, the agreement is qualitatively very good, with a slight systematic error in slope. Another alias is seen at z = 11, well before the simulated SFH.

Finally, the bottom right example is for a simulation of four episodes of star formation that are disjoint but likely associated: 2 × 10⁹ M_⊙ at z = 10, 8 × 10⁸ M_⊙ at z = 8, and two at ∼10⁸ M_⊙ at z = 7 and 6. Again, there is a slope error, scatter, and an alias at z = 9, but the character of the SFH is recovered successfully.

Our more extensive simulations, like those in Figure 6, show that SEDz* is successful at recovering the basic form of a wide variety of SFHs. With more experience and understanding of how to optimize this kind of testing and a richer, deeper sample than is now available, SEDz* should be able to characterize the distribution of SFHs and perhaps even how the distribution itself evolves over the epoch 5 < z < 12.

The purpose of this section is to use simulations to further illustrate how SEDz* computes an SFH by combining fluxes from stellar population templates to find the best fit to an observed SED. Specifically, we show that SEDz* finds a unique solution that is the best measurement of the SFH and related uncertainties, free of preconceptions or bias of what that SFH might be.

We begin this section by exploring, epoch by epoch, the SFHs of galaxies shown in Figure 3 that we described as having long and more or less continuous star formation for the ∼1 Gyr z = 12 to OE (5 < z < 7). In Figure 7, we show in three columns three examples from Figure 3: GLASS-5429, GLASS-3665, and GLASS-3459. For each, star formation begins at or before z = 12 and accumulates by OE to ∼1–2 × 10⁹ M_⊙. As we described in Section 2, SEDz* uses two sets of stellar population templates, for 10 Myr bursts and CSF, with fluxes recorded at integer redshifts z = 12–3. For each epoch, these "template tables" give the fluxes expected from 1 M_⊙ yr⁻¹ at that redshift and the subsequent fluxes for that evolving population.

Each of the three tests begins with the bottom panel, with OE at z = 11, and works up to an OE of z = 6 (left case) or 5 (middle and right cases). Using GLASS-5429 (left) as an example and starting at the bottom, SEDz* has attempted to fit the observed SED (OE, z = 6) only with stars born at z = 12 and 11. The best fit is for a mass of 9 × 10⁸ (right panel) at z = 12 and no z = 11 contribution, producing an SED (magenta line, middle panel) that is a very poor fit to the observed NIRCam fluxes (black points with errors), as indicated by inspection and the high ${\tilde{\chi }}^{2}$ (left panel). ²⁴ The next epoch up shows the best fit at z = 10 and demonstrates a key feature of SEDz*: the previous z = 12 and 11 contributions are not simply "carried over" to z = 10. Instead, these are recomputed with new mass estimates for the evolved fluxes from z = 12 and 11 star formation and adding more at z = 10. In other words, the amplitude of star formation at each previous epoch and the "present" epoch are vectors that the NNLS weighs and combines to find the lowest ${\tilde{\chi }}^{2}$. The fit at z = 10 to the observed SED is a bit better, and this improvement continues up the sequence of epochs as star formation is added at z = 9, 8, and 7. The best fit—the lowest ${\tilde{\chi }}^{2}$—that sets the redshift of the SED derives from contributions from most epochs. This contrasts with the situation for bursty SFHs, discussed in the next section.

The other two cases shown here, GLASS-3665 and GLASS-3459, show the same behavior; the good fit of the SED at the minimum ${\tilde{\chi }}^{2}$ develops from z = 11 as additional epochs of star formation are added. These three examples are typical of the full sample of eight galaxies shown in Figure 3; this "continuing growth" behavior played no role in their inclusion in our sample and indeed was not recognized before this "evolution of the SED" experiment. This is our primary evidence that the long, continuous SFHs in Figure 3 are genuine; achieving the best fit requires the addition of stellar populations over a wide redshift range.

In the previous section, we used SFH simulations to show that the long SFHs we have found are, in fact, reproduced by star formation over most of the epochs z = 12–5. Of course, some of these might be better described as a series of four to five bursts with significant drops in star formation in between. Since SEDz* cannot resolve such star formation episodes, a better way to answer this question could be to use a large sample to study the distribution of bursty histories to see if there is a continuum of SFHs that varies between one or two bursts from z = 12–5 and the long histories just discussed. Lacking such a sample at this time, it is useful to ask if well-separated bursts can mimic long histories, at least to the extent that most or all of the examples in Figure 3 could be "noisy" versions of bursts. We investigate this question in this section.

We rigged the SEDz* simulation program to produce random combinations of bursts with masses between 10^8.5 and 10^9.5 M_⊙; redshifts z = 8, 9, 10, 11, or 12; and observed at z = 5, 6, or 7. Fifty such combinations were made; the first 28, which are representative, are shown in Figure 8. For each combination, there is a pair of SEDz* plots: the burst parameters at the chosen epoch (top) and a companion SFH at the chosen OE (bottom). ²⁵ Scanning over these examples, it is readily seen that, at the redshift of the burst, the SFH is usually a single epoch or two adjacent epochs of star formation, and the observed starburst from a much later time is unresolved and typically splits into several epochs straddling the input epoch of the burst. At OE, the SED is always steep for these; there is no star formation that extends to OE at a significant level, so the relatively old population dominates.

Zoom In Zoom Out Reset image size

However, there are six cases out of 28 where we see something like the long SFHs of Figure 3, simulated in the previous section (examples enclosed in blue boxes). We extract these for comparison with six of the observed SFHs in Figure 9. We see that, although superficially like the long SFHs on the right, there are differences. Simulations 6, 14, 16, and 26 all have flat SFHs around the burst epoch (the "resolution" problem), with two rapidly falling values—with large error bars—to lower redshift. GLASS 2199 on the right is the closest to these in steepness, but its longer history is distinguishable. Simulations 8 and 15 are most like the righthand examples, which would make this possibility a <10% occurrence. In point of fact, there is another difference that separates even these mistaken burst SFHs from the examples in Figure 3: their SEDs are all considerably redder (with the exception of GLASS 2199). As expected, the SED of an old burst is not likely to mimic an SFH that has a significant proportion of relatively younger stars.

Zoom In Zoom Out Reset image size

We conclude from this test that a small fraction of strong bursts at higher redshift can appear to stretch to longer histories when viewed at z = 5–7, but they could only account for at most one or two of the long SFHs we have identified; because of the difference in the SED color, even these should be distinguishable.

Reviewing these tests, Figure 6 shows that we can simulate SFHs of the kinds seen in this paper and recover them with sufficient fidelity to distinguish one type from another. Figure 7 shows that the long SFHs we observe extending to lower z are in fact resolved by SEDz* into many epochs of star formation from z = 12 to OE. Figures 8 and 9 show through simulations that 8 < z < 12 bursts are unlikely to be mistaken for long SFHs. Even if a burst cannot be pinned down to a certain redshift and duration when viewed from 5 < z < 7, its derived SFH does not contaminate the sample of the other histories we have shown; it is identifiable as a likely burst in most cases. The ability of SEDz* to leverage stellar populations whose light is dominated by A stars to investigate the variety of SFHs at z > 5 is confirmed.