On the Stellar Populations of Galaxies at z = 9–11: The Growth of Metals and Stellar Mass at Early Times

We use the Flexible Stellar Population Synthesis (FSPS) package (Conroy et al. 2009; Conroy & Gunn 2010) with the Modules for Experiments in Stellar Astrophysics Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016). The MIST isochrones include the effects of rotation that boost the ionizing flux production of massive stars in a manner similar to the effect of binaries (Choi et al. 2017). Furthermore, throughout this work, we assume a Chabrier (2003) initial mass function.

Our fiducial physical model consists of 14 free parameters describing the contribution of stars, gas, and dust (Table 2). While not all 14 parameters are constrained by the photometric data, the use of a highly flexible model together with physically motivated priors prevents the results from being overinterpreted. Therefore, the choice of the priors is important. As we show in Section 5, a key conclusion of the paper is that the inferred early mass growth of galaxies heavily depends on the prior on the SFH.

Table 2. Summary of 14 Parameters and Priors for the Fiducial Physical Model with a Flexible Star Formation History (SFH) within Prospector

Parameter	Description	Prior
zphot	redshift	prior from EAZYor fixed to zspec
$\mathrm{log}(Z/{Z}_{\odot })$	stellar metallicity	uniform: $\min \ =\ -2.0$, $\max \ =\ 0.19$
$\mathrm{log}({M}_{\star }/{M}_{\odot })$	total stellar mass formed	uniform: $\min \ =\ 6$, $\max \ =\ 12$
SFH	flexible SFH: ratio of the SFRs in adjacent time bins of the NSFH-bin SFH (NSFH − 1 parameters total, with default choice NSFH = 6); parametric SFH: delayed-τ model with two free parameters	see Section 3.2 for details
n	power-law modifier to shape of the Calzetti et al. (2000) attenuation curve of the diffuse dust	uniform: $\min \ =\ -1.0$, $\max \ =\ 0.4$
${\hat{\tau }}_{\mathrm{dust},\ 2}$	diffuse dust optical depth	clipped normal: $\min \ =\ 0$, $\max \ =\ 4$, μ = 0.3, σ = 1
${\hat{\tau }}_{\mathrm{dust},\ 1}$	birth-cloud dust optical depth	clipped normal in (τdust, 1/τdust, 2): $\min \ =\ 0$, $\max \ =\ 2$, μ = 1, σ = 0.3
$\mathrm{log}({Z}_{\mathrm{gas}}/{Z}_{\odot })$	gas-phase metallicity	uniform: $\min \ =\ -2.0$, $\max \ =\ 0.5$
$\mathrm{log}(U)$	ionization parameter for the nebular emission	uniform: $\min \ =\ -4.0$, $\max \ =\ -1$
fIGM	scaling of the IGM attenuation curve	clipped normal: $\min \ =\ 0$, $\max \ =\ 2$, μ = 1, σ = 0.3

Download table as:  ASCIITypeset image

In our fiducial runs, we adopt the EAZY posterior (Section 2; F21) as a prior for the photometric redshift or fix the redshift to the spectroscopic redshift z_spec when available. Three galaxies have a spectroscopic redshift: EGS-6811 with z_spec = 8.678 (Zitrin et al. 2015), EGS-44164 with z_spec = 8.665 (R. L. Larson et al. 2022, in preparation), GOODSN-35589 with z_spec = 10.957 (Oesch et al. 2016; Jiang et al. 2020). The main motivation for us to assume the EAZY posterior as redshift prior is that it allows us to focus on posterior sampling on the stellar population part instead of the redshift space (Prospector is rather expensive to run in terms of time) and also to propagate the redshift uncertainty into the Prospector modeling. We model the chemical enrichment histories of the galaxies with a delta function, i.e., assuming that all stars within the galaxy have the same metal content with scaled-Solar abundances. This single metallicity is varied with a prior that is uniform in $\mathrm{log}(Z/{Z}_{\odot })$ between −2.0 and 0.19, where Z_⊙ = 0.0142 (Asplund et al. 2009).

One of the key strengths of the SED fitting code Prospector is the possibility to adopt flexible SFHs. Specifically, we adopt SFHs in our fiducial model, which do not assume a certain shape with time ¹⁷ and are simply partitioned into time bins. The SFHs are characterized by the ratios of the SFRs in adjacent time bins. There are five free parameters for six time bins in addition to the total stellar mass. Furthermore, we also explore a parametric, delayed-τ model (two free parameters). Details about the SFH prior are given in Section 3.2. For the total stellar mass M_⋆, we assume a flat prior in log-space in the range of $6\lt \mathrm{log}({M}_{\star }/{M}_{\odot })\lt 12$. Throughout this work, the stellar mass M_⋆ denotes the integral of the SFH, i.e., it is the mass of all stars ever formed.

We model dust attenuation using a two-component dust attenuation model with a flexible attenuation curve (see Charlot & Fall 2000). The first component is a birth-cloud component in our model that attenuates nebular emission and stellar emission only from stars formed in the last 10 Myr (attenuation law is a power law with a slope of −1). The second component is a diffuse component that has a variable attenuation curve and attenuates all stellar and nebular emission from the galaxy. We use the prescription from Noll et al. (2009) with a Kriek & Conroy (2013) attenuation curve, where the slope n of the curve (dust index) is a free parameter and is directly linked to the strength of the UV bump. The dust index n is modeled as an offset from the slope of the UV attenuation curve from Calzetti et al. (2000). In total, the attenuation prescription has three free parameters: (i) the slope n (flat prior between −1.0 < n < 0.4); (ii) the normalization ${\hat{\tau }}_{\mathrm{dust},\ 2}$ of the diffuse dust component (flat prior between $0\lt {\hat{\tau }}_{\mathrm{dust},\ 2}\lt 4.0$); and (iii) the normalization ${\hat{\tau }}_{\mathrm{dust},\ 1}$ of the birth-cloud component, which we model as a ratio with respect to the diffuse component (prior is a clipped normal centered on 1 with width of 0.3 in the range of $0\lt {\hat{\tau }}_{\mathrm{dust},\ 1}/{\hat{\tau }}_{\mathrm{dust},\ 2}\lt 2.0$, motivated by Calzetti et al. (1994), Price et al. (2014)).

The nebular emission (emission lines and continuum) is self-consistently modeled (Byler et al. 2017). We have two parameters: the gas-phase metallicity (Z_gas) and the ionization parameter (U). We assume a flat prior in log-space for the metallicity ($-2.0\lt \mathrm{log}({Z}_{\mathrm{gas}}/{Z}_{\odot })\lt 0.5$) and ionization ($-4\lt \mathrm{log}(U)\lt -1$) parameters. Importantly, we do not link the gas-phase metallicity Z_gas and the stellar metallicity Z_⋆, i.e., Z_gas and Z_⋆ are decoupled from each other. We choose to decouple the gas-phase from the stellar metallicity because it allows us to cover both cases where the Z_gas is smaller or larger than the Z_⋆. Both cases are expected in the evolution of galaxies. Specifically, in the case of a closed-box chemical model, we expect the stellar metallicity to be always smaller than the gas-phase metallicity. This might be true in certain phases of the galaxy's lifetime. However, galaxies also accrete new gas, which typically has a lower metallicity than gas (and stars) already present in the galaxy, leading to a gas-phase metallicity that can be lower than the stellar metallicity. The main consequence of this assumption is an overall larger flexibility in the SED modeling, in particular in regards to the emission line strengths.

As the photometry probes the rest-frame λ < 1216 Å spectrum at high redshifts, we include a z-dependent IGM attenuation following Madau (1995). This includes a free parameter that scales the total IGM opacity (f_IGM), intended to account for line-of-sight variations in the total opacity. We adopt for f_IGM a clipped Gaussian prior distribution centered on 1, with a dispersion of 0.3 and clipped at 0 and 2.

In order to explore the robustness of our inferred mass assembly histories, we want to explore the dependence of our results on the assumed SFH prior. As noted above, the strength of Prospector is the possibility of adopting a flexible SFH (see also Iyer & Gawiser 2017; Iyer et al. 2019 for another SED fitting code with a flexible SFH approach). We assume four different priors for the SFH: three are flexible SFHs ("continuity prior", "bursty continuity prior", and "Dirichlet prior"), while the fourth is a parametric SFH with the shape of the delayed-τ model ("parametric prior"). The strength of the flexible SFH priors is that they are not parametric functions of time (in contrast to the parametric delayed-τ prior), which allows for a large flexibility regarding the shape of the SFH. Figure 1 illustrates the behavior of these four different priors by plotting the median trend of the SFH and individual draws.

Zoom In Zoom Out Reset image size

For the flexible SFHs, we assume that the SFH can be described by N_SFH time bins, where the SFR within each bin is constant. We fix N_SFH = 6 and specify the time bins in lookback time. The first bin is fixed at 0–10 Myr to capture variation in the recent SFH of galaxies, while the other bins are spaced equally in logarithmic time between 10 Myr and a lookback time that corresponds to z = 20, i.e., we assume the SFR = 0 M_⊙ yr⁻¹ at z > 20 (a reasonable assumption given what observational constraints and theoretical predictions exist for this epoch; Maio et al. 2010; Bowman et al. 2018; Jaacks et al. 2019). These time bins are plotted as vertical dashed yellow lines in Figure 1. We fit for the ratio between the time bins (N_SFH − 1 free parameters) and the total stellar mass formed, which has a flat prior in log-space in the range of $6\lt \mathrm{log}({M}_{\star }/{M}_{\odot })\lt 12$.

Impact of the choice regarding the number of time bins has been extensively discussed in (Leja et al. 2019a; see their Appendix A). They explored varying the number of time bins between N_SFH = 4–14 and show that the results of the mock analysis are largely insensitive to the number of bins as long as N_SFH ≳ 4. Although their mock analysis cannot be translated to our analysis one-to-one, we think that this main conclusion still holds because they investigated lookback times of 10 Gyr, much longer than the age of the universe at the epochs of our objects. Therefore, our log-spaced bins with a width of typically less than 100 Myr should be enough to convey all of the necessary information in the data.

The priors for flexible SFHs are extensively discussed in (Leja et al. 2019a; see also Tacchella et al. 2022), while parametric SFHs are explored in Carnall et al. (2019). We adopt a continuity prior as well as the Dirichlet prior. For the continuity prior, we directly fit for the ${\rm{\Delta }}\mathrm{log}(\mathrm{SFR})$ between adjacent time bins. We adopt the Student's t-distribution ${\rm{\Delta }}\mathrm{log}(\mathrm{SFR})$. For the "continuity prior", we assume a Student's t-distribution with σ = 0.3 and ν = 2, which weights against sharp transitions and is motivated by simulated SFHs at z ∼ 1 (Leja et al. 2017). This is our fiducial SFH prior. For the "bursty continuity prior", we adopt σ = 1.0 and ν = 2, which leads to a more variable (i.e., bursty) SFH. In the case of the "Dirichlet prior", the fractional sSFR for each time bin follows a Dirichlet distribution (Leja et al. 2017). We assume a concentration parameter of 1, which weights toward smooth SFHs. As shown in Figure 1, both the continuity and the Dirichlet prior include a symmetric prior in age and sSFR and an expectation value of constant SFR(t). The key difference from the Dirichlet prior is that the continuity prior explicitly weights against sharp changes in SFR(t).

Finally, for the parametric SFH, we assume a delayed-τ model of the form:

$$\begin{eqnarray}&&\mathrm{SFR}(t)=(t-{t}_{{\rm{a}}}){e}^{-(t-{t}_{{\rm{a}}})/\tau }.\end{eqnarray} \tag{ 1 }$$

The parameter τ is varied as $\mathrm{log}(\tau )$ within a uniform prior in the range of $-1.0\lt \mathrm{log}(\tau )\lt 10.0$, and the parameter t_a with a uniform prior between 1 Myr and the age of the universe at the galaxies' redshift z_phot (t_H(z_obs)). Despite this large prior space for the parameters τ and t_a, the resulting SFH from the parametric prior follows a specific shape of an increasing SFH with time, as shown in Figure 1, consistent with constraints on SFHs in the epoch of reionization (e.g., Papovich et al. 2011).

After setting up the physical galaxy SED model with 14 free parameters, we fit this model to the photometric data (Section 2) within the Prospector framework using the dynamic nested sampling algorithm dynesty (Speagle 2020), which allows us to perform an efficient sampling of the high-dimensional and complex parameter space. A strength of Prospector together with dynesty is its ability to infer full posterior distributions of the SED parameters and their degeneracies. We discuss these SED parameters and their inferred properties, such as the stellar mass, metallicity, and SFH, in the upcoming sections, but see Table 3 for a summary of the main physical parameters. Here, we focus on the resulting SEDs and compare them with the measured photometry in order to assess the goodness of the fits. Then we briefly discuss the resulting posteriors of the dust attenuation parameters.

3.3.1. SEDs and Goodness of Fit

Figure 2 shows the observed and modeled posterior SEDs for our 11 z = 9–11 galaxy candidates. The blue circles show the detected photometric bands, while the arrows mark the upper limits. The red solid lines and shaded regions indicate the median and the 16–84th percentile of the posterior SED. These are the results for our fiducial SED run, where we adopt the continuity prior for the SFH and the EAZY posterior as the redshift prior (if no spectroscopic redshift is available). Although we include emission lines in all the fits (Section 3.1), they are typically not visible in the posterior SED because they are "smeared" out when marginalizing over the redshift posterior distribution. Nominal exceptions are the galaxies for which we fix the redshift to the spectroscopic redshift (e.g., EGS-6811, EGS-68560, and GOODSN-35589): for those objects, the emission lines are clearly visible.

Zoom In Zoom Out Reset image size

The detections below the Lyα break in COSMOS-20646 and COSMOS-47074 are discussed in F21 (see their Figure 15). Briefly, the detections are only marginally significant (S/N = 2–3) and are also slightly offset from the source position. Consistent with the Prospector analysis here (see Section 3.4), when including these fluxes in the photometric redshift modeling (along with all the nondetections/upper limits), F21 also found the preference is still for a high-redshift solution, although a low-redshift solution is possible for both sources at a low (10%) probability level.

The lower panels of Figure 2 show the χ values for the individual passbands, which is the difference between observed to model fluxes normalized by the observed errors. The individual χ values are typically around 1. We also quote the total ${\chi }_{\mathrm{tot}}^{2}$, which we estimate by summing the individual χ values of the detected bands. We do not quote the reduced χ² values as the number of degrees of freedom is not well defined in a nonlinear model, such as considered here (Andrae et al. 2010). In summary, the model is able to reproduce the observational data well within the observational uncertainties.

These results are for our fiducial, continuity SFH prior, but the other SFH priors are able to reproduce the observational data equally well. Specifically, we find very similar (differences amount to less than 20%) χ values for all the four SFH priors (Section 3.2). Furthermore, none of these priors is preferred by the data: the Bayes factor, i.e., the ratio of the evidences between the different models, is around 1. Specifically, the median and 16–84th percentile of Bayes factor of the bursty continuity prior, of the Dirichlet prior and of the parametric prior (all relative to the fiducial continuity prior) are ${0.7}_{-0.1}^{+0.7}$, ${0.9}_{-0.3}^{+0.1}$, and ${1.1}_{-0.3}^{+0.3}$, respectively. This inability to identify a preferred model can be attributed to both the small sample size and the limited information content of the observational data.

3.3.2. Attenuation Curve

As highlighted in the previous section, the rest-frame wavelength coverage is limited due to the high-redshift nature of these sources. Specifically, we only cover the rest-frame UV and the Balmer/4000 Å break, though the latter is only constrained by the IRAC photometry, which suffers from systematic uncertainties related to deblending (see also F21 and Appendix A) and can also be contaminated by strong emission lines. Therefore, in order to constrain the buildup of stellar mass (i.e., SFHs), we need to properly interpret the rest-UV emission.

Flexibility in the attenuation law is motivated from observations (e.g., Johnson et al. 2007; Kriek & Conroy 2013; Battisti et al. 2016; Salmon et al. 2016; Salim et al. 2018) and theory (e.g., Seon & Draine 2016; Narayanan et al. 2018; Shen et al. 2020). Specifically, Katz et al. (2019b) use a cosmological radiation hydrodynamics simulation to show that dust preferentially resides in the vicinity of the young stars, thereby increasing the strength of the measured Balmer break. Therefore, we adopt a flexible attenuation law (Section 3.1) so that we can marginalize over the uncertainty of an unknown attenuation law when constraining the SFHs and stellar ages among other physical properties.

Figure 3 shows the resulting posterior distributions of the attenuation as a function of the wavelength of all galaxies in our sample. The median and the 16–84th percentiles are shown as solid red lines and red shaded regions, respectively. For comparison, we also plot the Calzetti et al. (2000) attenuation curve and an SMC-like (Pei 1992) attenuation curve. Furthermore, the blue solid line marks the median of the prior, while the blue shaded region indicates the 16–84th percentile of the prior.

Zoom In Zoom Out Reset image size

We find that the attenuation laws of our galaxies are consistent with a Calzetti et al. (2000) law with R_V = 3.1, though significant variations from galaxy to galaxy are present. Interestingly, our obtained attenuation laws are typically shallower than the input prior. Furthermore, the uncertainty for individual galaxies is large, indicating that the attenuation law is not well constrained by our data and that degeneracies with, for example, the stellar age and metallicity exist (see also Figure 22). Nevertheless, the posteriors of the attenuation laws of individual galaxies look physically sensible, which supports the choice of priors (Table 2). Importantly, in the remainder of this paper, we marginalize over the uncertainty in the attenuation parameters (i.e., attenuation law).

As discussed in Section 2 and F21, the photometric-redshift code EAZY has been used to perform our redshift estimation and to select our candidate z = 9–11 galaxies. Although EAZY allows linear combinations of any number of provided templates, the explored parameter space is limited. In this section, we explore the photometric-redshift constraints that we obtain with Prospector and compare them with the EAZY-based photometric redshifts, finding excellent agreement.

In order to fit for the photometric redshift and also allow low-z solutions, we have to extend the model that we introduced in Section 3.1. We call this the "free-z run", where we let the redshift be free and assume a flat prior between 0.1 and 13. Second, we add dust emission and active galactic nucleus (AGN) emission in order to add flexibility to the SED in the infrared in order to reproduce dusty, lower-z galaxies that have similar SEDs as the high-z dropout candidates. In particular, the dust and AGN emission can dominate the near-IR flux, i.e., the emission in the IRAC bands at low redshifts. At higher redshifts (z > 3), this emission contributes to bands at longer wavelengths than IRAC covers; hence, we do not have to consider this emission in our fiducial model.

We follow the description of Leja et al. (2017, 2018) and Tacchella et al. (2022), which adds five new free parameters to our fiducial model, giving us a model with a total of 19 free parameters. Briefly, the three new parameters for the dust emission are γ_e (mass fraction of dust in high radiation intensity; log-uniform prior with minimum and maximum of 10⁻⁴ and 0.1), ${U}_{\min }$ (minimum starlight intensity to which the dust mass is exposed; clipped normal prior with a mean of 2, a standard deviation of 1, minimum and maximum of 0.1 and 15), q_PAH (percent mass fraction of PAHs in dust, uniform prior with minimum and maximum of 0.5 and 7.0). The two new parameters for the AGN emission are f_AGN (AGN luminosity as a fraction of the galaxy bolometric luminosity, log-uniform prior with minimum and maximum of 10⁻⁵ and 3) and τ_AGN (optical depth of AGN torus dust, log-uniform prior with minimum of 5 and maximum of 150).

Figure 4 shows the photometric redshift posteriors obtained from EAZY (black lines) and Prospector (red lines). The results from Prospector assume the aforementioned free-z run with a uniform redshift prior between z = 0.1–13. Three galaxies (EGS-6811, EGS-68560, and GOODSN-35589) have spectroscopic redshifts, which are indicated in blue. We find overall good agreement between EAZY and Prospector. The two approaches return photometric redshifts that are consistent with each other within the uncertainty. This can also be seen directly in Figure 5, which shows a comparison of the Prospector-based and the EAZY-based photometric redshifts. The circles and the error bars indicate the median and 16–84th percentiles of the redshift posterior, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Importantly for this work here, Prospector confirms the high-redshift nature of these galaxies. The probability of lying beyond z = 8, P(z > 8), is larger than 90% for all galaxies except EGS-20381 (P(z > 8) = 0.75) and EGS-40898 (P(z > 8) = 0.71), for which a tail in the z_phot posterior toward z ∼ 6 exists. Furthermore, for the galaxies that show minor peaks at z ∼ 1–3, the posteriors all have P(z < 6) < 0.1, i.e., <10% of the posterior volume is at low redshift.

As mentioned above, we assume a flat redshift prior. This might actually not be the ideal prior as a flux-limited survey will contain many more low-z than high-z galaxies. We could therefore think of more complicated priors that, for example, weight according to the luminosity function and consider also the selection function of the survey. A detailed investigation of this is beyond the scope of this work. Nevertheless, as we have probably significantly overestimated the high-z prior volume, we have also performed two additional "free-z run", where we split the redshift prior in half by running z_phot in the range of 0.1–7.0 and 7.0–13.0. Although for some galaxies a viable low-z solution is identified, the high-z solution is preferred for all objects considering both the total ${\chi }_{\mathrm{tot}}^{2}$ values as well as the Bayes factor (i.e., ratio of the evidences in the Bayesian analysis). Specifically, we find that the high-z run has a lower ${\chi }_{\mathrm{tot}}^{2}$ by a factor of 2–4 than the low-z run. Furthermore, we obtain for the Bayes factor a median of ${{ \mathcal Z }}_{\mathrm{low}-{\rm{z}}}/{{ \mathcal Z }}_{\mathrm{high}-{\rm{z}}}=3\times {10}^{-2}$ over the whole sample, with all galaxies ${{ \mathcal Z }}_{\mathrm{low}-{\rm{z}}}/{{ \mathcal Z }}_{\mathrm{high}-{\rm{z}}}\lt 0.5$. This shows that the high-z model is preferred for all galaxies, which is consistent with our findings above.

Finally, as mentioned in F21, we have performed the IRAC photometric deblending in two different ways, once with TPHOT (fiducial) and once with GALFIT. We discuss the results of changing from TPHOT to GALFIT IRAC photometry in Appendix A and Figure 18. In summary, we find consistent photometric redshift estimates for all galaxies except COSMOS-20646 (${z}_{\mathrm{phot}}={2.63}_{-0.16}^{+0.16}$) and EGS-20381 (${z}_{\mathrm{phot}}={6.79}_{-0.13}^{+0.57}$). This is consistent with the EAZY-based results with this photometry from F21 for these two objects. We also find a significant difference for EGS-6811 (${z}_{\mathrm{phot}}={7.40}_{-0.04}^{+0.05}$), where this alternative z_phot estimate is inconsistent with the available spectroscopic redshift (z_spec = 8.68).

4.1.1.  β as a Proxy for Dust Attenuation

The UV spectral slope β (defined as f_λ ∝ λ^β ; Calzetti et al. 1994) is often used to quantify stellar populations in the high-redshift universe as it is a straightforward probe of the color of the emergent light from the young, massive stars in these early galaxies. It is also a relatively easy measurement—β is readily measurable if a given galaxy has detections in at least two photometric bands probing the rest-frame UV (free of both the Lyα break introduced by the neutral IGM and Balmer/4000 Å break). While a number of physical factors can affect the rest-UV color, the observed slope β is generally interpreted as a proxy for dust attenuation.

This dust-heavy interpretation of the UV slope is especially true at the highest redshifts we discuss here, as the stellar ages are limited by the very short time since the end of the cosmic dark ages, and metallicities are similarly limited both by the lack of time for significant chemical enrichment and the relative insensitivity of β to changing stellar metallicities. In Figure 6, we use Bruzual & Charlot (2003) models to show how the inferred value of the UV spectral slope β changes with increasing dust attenuation, stellar population age, and metallicity. For a given curve showing the change in one property, we keep the other two properties fixed, with the fixed values being A_V = 0.2 mag, age = 50 Myr, and Z = 0.2 Z_⊙. For dust attenuation, we consider both a starburst (Calzetti et al. 2000) and an SMC-like (Pei 1992) attenuation curve, while for stellar population age we consider a constant (τ = ∞ ), rising (τ = −300 Myr), and extreme burst (τ = 0.1 Myr) SFH.

Zoom In Zoom Out Reset image size

Figure 6 shows that β is typically much more sensitive to dust attenuation than it is to age or metallicity, similar to previous analyses (e.g., Cortese et al. 2008; Bouwens et al. 2009; Wilkins et al. 2016; Jaacks et al. 2018; Tacchella et al. 2018b). The UV slope β does get redder with increasing metallicity or stellar population age (at fixed dust attenuation), but the changes are relatively small. The change from Z = 0.005 to 1.0 Z_⊙ is Δβ = 0.4, while the change from t = 10 Myr to ∼200 Myr (representing a formation redshift of z = 13 for an observation redshift of z = 9) is Δβ = 0.3. In comparison, a change in SMC-law V-band dust attenuation from 0 to 0.7 mag results in Δβ = 1.3. The exception is for the burst SFH, where all of the stellar mass is formed in 0.1 Myr. This population is still fairly blue at t = 10 Myr (β = −2.4) but becomes very red (β = −1.4) by t = 100 Myr. However, as this galaxy has not formed any stars since its initial burst, its luminosity fades rapidly, dropping three magnitudes from t = 10 to t = 100 Myr. Such a galaxy, at $\mathrm{log}({M}_{\star }/{M}_{\odot })=9.5$, would be below the detection limits of our sample. This highlights that precise measures of β can be very sensitive to changes in dust attenuation, especially at these high redshifts where changes in the UV slope due to stellar ages are minimized due to the young age of the universe.

However, β can still inform on chemical enrichment; Figure 6 highlights that, for very young and dust-free galaxies, β reaches a minimum of −2.7 for Z = 0.005 Z_⊙. While this minimum is somewhat model dependent, the search for galaxies with even bluer spectral slopes (β ≲ −3) has been an important part of high-redshift studies since the advent of deep near-IR imaging (e.g., Bouwens et al. 2010; Finkelstein et al. 2010). Such blue values would indicate ultra-poor metallicities, potentially even metal-free Population III galaxies. The likelihood for such a discovery is complex, however, as enrichment from the initial burst of Population III stars alone may significantly redden the observed colors of galaxies (e.g., Jaacks et al. 2018). Nonetheless, our sample of well-observed z ∼ 9–10 galaxies presents an excellent opportunity to measure the UV slope β and constrain chemical enrichment at some of the highest redshifts yet probed.

4.1.2. Measurements of β

While the UV slope β can in principle be measured by a single color, additional colors increase the accuracy of the resulting measurements. Finkelstein et al. (2012) performed simulations to assess best practices for measuring this quantity, comparing a single color, a power-law fit to multiple colors, and measuring β directly from the best-fitting SED model spectrum. They found that, when many colors are available (e.g., at lower redshifts), both the power-law and SED-fitting method outperform the single-color method, while at higher redshifts when information is more limited, the SED-fitting method results in both smaller scatter and smaller bias. We therefore elect to use the SED-fitting method here. We note that as our galaxies are fairly bright, we do not expect photometric scatter to result in a bias toward bluer measured UV slopes as found for fainter galaxies (Dunlop et al. 2012).

This calculation is done by taking the Prospector model spectra (using the fiducial fit with the continuity SFH prior), converting them to f_λ in the rest frame, and fitting a slope to the spectrum in wavelength windows from Calzetti et al. (1994) designed to omit spectral emission and absorption features. The Calzetti et al. (1994) windows span 1268–2580 Å; however here we omit the three bluest windows to avoid potential contamination from the Lyα break due to the photometric redshift uncertainties, so our bluest window begins at 1407 Å. We apply this measurement to the spectra from 100 random draws of the Prospector posteriors such that the uncertainty on β includes all model uncertainties (including uncertainties on the redshift when relevant). From these 100 draws, we calculate the median value and 68% confidence range on β.

The results for each galaxy are listed in Table 3. While our measured values span a wide range, interestingly all galaxies have β > −2.5, implying measurable dust attenuation in every galaxy in our sample. While this may not be unexpected in such relatively massive systems, it implies that significant dust production must be taking place at z > 10 to be observable in this epoch.

Table 3. Results for the Main Physical Parameters from our Fiducial Run Assuming the Continuity Prior

ID	Redshift	MUV, obs	UV slope β	$\mathrm{log}\,{M}_{\star }$	$\mathrm{log}\,{\mathrm{SFR}}_{50}$	$\mathrm{log}\,{\mathrm{sSFR}}_{50}$	AV	$\mathrm{log}\,Z$
		(mag)		(M⊙)	(M⊙ yr−1)	(Gyr−1)	(mag)	(Z⊙)
COSMOS-20646	${9.77}_{-0.19}^{+0.19}$	$-{22.10}_{-0.09}^{+0.09}$	$-{0.62}_{-0.12}^{+0.11}$	${10.9}_{-0.2}^{+0.2}$	${2.5}_{-0.3}^{+0.3}$	${0.7}_{-0.4}^{+0.3}$	${0.9}_{-0.5}^{+0.6}$	$-{0.4}_{-0.7}^{+0.3}$
COSMOS-47074	${9.83}_{-0.39}^{+0.28}$	$-{21.01}_{-0.12}^{+0.12}$	$-{2.11}_{-0.18}^{+0.20}$	${9.3}_{-0.4}^{+0.3}$	${1.0}_{-0.2}^{+0.3}$	${0.8}_{-0.4}^{+0.3}$	${0.1}_{-0.1}^{+0.3}$	$-{1.2}_{-0.6}^{+0.7}$
EGS-6811	zspec = 8.68	$-{22.10}_{-0.05}^{+0.05}$	$-{1.61}_{-0.12}^{+0.18}$	${10.6}_{-0.3}^{+0.2}$	${2.0}_{-0.3}^{+0.4}$	${0.5}_{-0.4}^{+0.3}$	${0.7}_{-0.4}^{+0.6}$	$-{0.5}_{-0.7}^{+0.4}$
EGS-44164	zspec = 8.66	$-{21.87}_{-0.05}^{+0.05}$	$-{1.87}_{-0.11}^{+0.11}$	${10.2}_{-0.2}^{+0.2}$	${1.6}_{-0.3}^{+0.4}$	${0.5}_{-0.4}^{+0.4}$	${0.4}_{-0.2}^{+0.5}$	$-{1.3}_{-0.4}^{+0.6}$
EGS-68560	${9.16}_{-0.24}^{+0.28}$	$-{21.47}_{-0.09}^{+0.09}$	$-{2.37}_{-0.11}^{+0.15}$	${9.1}_{-0.3}^{+0.2}$	${1.0}_{-0.2}^{+0.2}$	${1.0}_{-0.2}^{+0.2}$	${0.1}_{-0.0}^{+0.1}$	$-{1.5}_{-0.3}^{+0.5}$
EGS-20381	${8.51}_{-0.49}^{+0.33}$	$-{21.12}_{-0.14}^{+0.18}$	$-{1.60}_{-0.20}^{+0.22}$	${10.0}_{-0.4}^{+0.3}$	${1.5}_{-0.3}^{+0.4}$	${0.6}_{-0.4}^{+0.4}$	${0.6}_{-0.3}^{+0.5}$	$-{1.0}_{-0.6}^{+0.7}$
EGS-26890	${9.06}_{-0.28}^{+0.29}$	$-{21.25}_{-0.10}^{+0.10}$	$-{1.94}_{-0.21}^{+0.14}$	${9.6}_{-0.3}^{+0.3}$	${1.3}_{-0.3}^{+0.3}$	${0.7}_{-0.4}^{+0.3}$	${0.2}_{-0.2}^{+0.4}$	$-{1.1}_{-0.6}^{+0.7}$
EGS-26816	${9.40}_{-0.33}^{+0.36}$	$-{21.28}_{-0.13}^{+0.13}$	$-{1.80}_{-0.24}^{+0.33}$	${9.7}_{-0.4}^{+0.3}$	${1.4}_{-0.3}^{+0.3}$	${0.8}_{-0.4}^{+0.3}$	${0.3}_{-0.2}^{+0.4}$	$-{1.1}_{-0.6}^{+0.8}$
EGS-40898	${8.48}_{-0.54}^{+0.39}$	$-{20.70}_{-0.16}^{+0.23}$	$-{1.45}_{-0.28}^{+0.28}$	${9.9}_{-0.4}^{+0.4}$	${1.4}_{-0.4}^{+0.4}$	${0.6}_{-0.4}^{+0.4}$	${0.6}_{-0.3}^{+0.6}$	$-{1.0}_{-0.7}^{+0.6}$
GOODSN-35589	zspec = 10.96	$-{21.71}_{-0.09}^{+0.08}$	$-{2.34}_{-0.14}^{+0.13}$	${9.1}_{-0.2}^{+0.3}$	${1.1}_{-0.2}^{+0.2}$	${1.0}_{-0.2}^{+0.2}$	${0.1}_{-0.0}^{+0.1}$	$-{1.1}_{-0.6}^{+0.8}$
UDS-18697	${9.54}_{-0.10}^{+0.09}$	$-{22.14}_{-0.07}^{+0.07}$	$-{1.41}_{-0.06}^{+0.08}$	${11.0}_{-0.2}^{+0.4}$	${1.7}_{-0.6}^{+0.6}$	$-{0.4}_{-0.5}^{+0.4}$	${0.5}_{-0.4}^{+0.9}$	$-{0.0}_{-0.5}^{+0.1}$

Note. The values are the median of the posterior, while the errors indicate the 16–84th percentiles. We list the galaxy identifier (ID), the redshift (photometric if not z_spec specified), the absolute UV magnitude at rest-frame 1500 Å (M_{UV, obs}), the UV spectral slope (β), the stellar mass (M_⋆), the SFR averaged over past 50 Myr (SFR₅₀), the specific SFR averaged over past 50 Myr (sSFR₅₀), the dust attenuation at 5500 Å (A_V), and the stellar metallicity (Z).

Download table as:  ASCIITypeset image

The left panel of Figure 7 compares these β measurements to each galaxy's rest-UV absolute magnitude (taken from F21), compared to previous results at z ∼ 9–10 from Wilkins et al. (2016) and Bhatawdekar & Conselice (2021), as well as to the derived trends between β and M_UV from Bouwens et al. (2014). As our galaxies span a relatively small dynamic range in M_UV, there is no correlation visible within our small sample (Pearson r = −0.32), though the bulk of our galaxies have measured β consistent with similarly bright galaxies at lower redshifts (z ≈ 4–8). Our faintest galaxies, at M_UV ∼ −21 have colors that are also consistent with the rest-UV colors measured for a stack of bright z ∼ 8 galaxies from Stefanon et al. (2021), who measured J − H ∼0 (for β ∼ −2) at M_UV = −21.

Zoom In Zoom Out Reset image size

In the right panel of Figure 7, we plot β versus our Prospector-derived stellar mass, also including points from Bhatawdekar & Conselice (2021) and the derived correlations between β and stellar mass at lower redshifts from Finkelstein et al. (2012). We see that our sample of z = 9–10 galaxies appears to exhibit a strong correlation between the stellar mass and β (Pearson r = 0.85), where more massive galaxies have redder UV spectral slopes, similar to the correlations found by Finkelstein et al. (2012). Furthermore, the measured values of β for our sample of galaxies at $\mathrm{log}({M}_{\star }/{M}_{\odot })\gtrsim 9.5$ are consistent with those measured for similarly massive galaxies at z = 4–8. We simulated whether photometric scatter could cause this trend and found that, at the brightness range of our sample, scatter does not appear to input any bias.

We explore this further in Figure 8, where we plot the median values of β from z = 4–8 from Finkelstein et al. (2012) in mass bins versus redshift alongside the results from our z = 9–11 galaxy sample. We calculate the median β of our sample by combining the measured values of β from the 100 posterior draws for each object into a single array and then calculating the median and 68% confidence range. For our full sample of 11 galaxies, this calculation yields a median of β = −1.76${}_{-0.49}^{+0.42}$. However, we acknowledge that two of our sources appear fairly red and also have the highest derived stellar masses (UDS-18697 and COSMOS-20646). As discussed further in F21, the proximity to bright neighbors makes the IRAC photometry for these sources less reliable; thus it is possible that residual light from the neighbors is contributing to the high stellar mass measurement. If we exclude these two galaxies from this median measurement, we find $\beta \,=-{1.87}_{-0.43}^{+0.35}$. While the median is not highly dependent on the inclusion of these two galaxies, we consider this latter value our fiducial value.

Zoom In Zoom Out Reset image size

Comparing our z = 9–11 galaxies to the results at lower redshift for similarly massive galaxies in Figure 8, we find the surprising result that even though we are probing a few 100 Myr closer to the Big Bang, these relatively massive galaxies appear similarly red to comparable mass galaxies at lower redshifts. This implies that not only does dust build up to significant values very rapidly in modestly massive galaxies, but that this level of attenuation is relatively invariant with redshift at 4 < z < 10 at these fixed high stellar masses.

In Figure 9, we compare our measured values of β to the Prospector-derived values of the V-band dust attenuation (A_V), stellar metallicity, and stellar age t₅₀. ¹⁸ The stellar age t₅₀ is the lookback time at which 50% of the stellar mass has formed, and it is very similar to the mass-weighted age of the SFH. Starting with dust attenuation, our sample of 11 galaxies exhibits a strong, and nearly monotonic, positive correlation between dust attenuation and β. This is consistent with what we expected from Figure 6, which implied that β should inform most strongly about dust attenuation. With the exception of the two bluest galaxies, all galaxies are constrained (at the 1σ level) to have nonzero levels of dust attenuation.

Zoom In Zoom Out Reset image size

While the lack of a strong observed correlation between β and the stellar metallicity in the middle panel is not surprising, as Figure 6 shows that β is not very sensitive to changes in stellar metallicity, we do find that our bluest galaxy (EGS-68560) has the tightest constraints on a low metallicity with Z < 0.1 Z_⊙ (1σ). This is consistent with the idea that very blue values of β (<−3) will imply very low metallicities. While we do not yet see such blue galaxies in our sample, as noted above these are fairly massive systems thus one might expect to see such blue colors (and thus relatively unenriched galaxies) at lower masses at this same epoch. In the right-hand panel, we see a similar result, where the bluest galaxies have the tightest constraints on a young average age, while at β > −2, there is little correlation. However, as discussed further in Section 5, these ages are highly dependent on the SFH prior.

The average inferred dust attenuation in our sample of 11 galaxies is A_V = 0.4 ± 0.3 mag. This is larger (though only at ∼1σ significance) than the average attenuation found in a sample of four fainter galaxies by Wilkins et al. (2016), who found an average A₁₅₀₀ = 0.5 ± 0.3, which corresponds to A_V = 0.12 ± 0.07. This is consistent with the expectation observed at lower redshifts that fainter galaxies have less dust attenuation, though we note that Wilkins et al. (2016) used the locally derived relation by Meurer et al. (1999) to convert between β and the dust attenuation (while our attenuation is derived from SED fitting using a flexible attenuation curve). However, this is not surprising as these fainter galaxies are presumably less massive, and simulations (e.g., Graziani et al. 2020) predict that the dust mass grows rapidly at higher stellar masses.

One of the main conclusions we can make in this section is that the rest-UV colors of z = 9–10 galaxies at M_UV < −21 and $\mathrm{log}({M}_{\star }/{M}_{\odot })=9\mbox{--}10$ are similar to those at the same UV luminosities and masses at z = 4–8. This has implications for the interpretation of the evolution of the UV luminosity function. Evidence has been growing that the bright end of the rest-UV luminosity function changes little, if at all, from z = 7 to 10. This idea was introduced by Bowler et al. (2020), and even more recent luminosity function measurements (including using this same sample here) are continuing to find a higher-than-expected number density of these bright systems (e.g., Rojas-Ruiz et al. 2020; Finkelstein et al. 2021). As the bright end of the luminosity function does appear to decline in abundance from z = 4 to 7 (e.g., Bouwens et al. 2015; Bowler et al. 2015; Finkelstein et al. 2015a), this apparent lack in evolution at higher redshift points to a physical change in the galaxies themselves.

The most obvious potential physical change would be in dust attenuation: if more distant galaxies at fixed UV magnitude are less attenuated, then the bright end of the UV luminosity function would evolve more slowly than the faint end, which is exactly what observations suggest. However, our results here cast doubt on this physical interpretation, as we find that the bulk of these bright massive galaxies have similar UV spectral slopes as their z ∼ 7–8 counterparts, and thus by extrapolation likely have similar levels of dust attenuation. While our sample is small, if this result holds with larger samples from robust observations with the James Webb Space Telescope (JWST) it implies another physical explanation will be needed, such as changes in the star formation efficiency (e.g., Finkelstein et al. 2015b; Stefanon et al. 2019), or time/mass scales for the onset of quenching (e.g., Tacchella et al. 2018a; Bowler et al. 2020).

The most surprising result from this analysis is that the UV spectral slope β for relatively massive UV-selected systems ($9\lt \mathrm{log}({M}_{\star }/{M}_{\odot })\lt 10$) changes very little from z ∼ 4 to z ∼ 9–10. This implies that the dust attenuation at this fixed stellar mass is roughly constant with redshift (though, as the effects of reddening due to age and stellar metallicity should be less at higher redshift, it is possible the actual attenuation at the highest redshifts is even higher at fixed β). Although the most recent constraints from the Atacama Large Millimeter Array imply that dust-obscured star formation is not dominant in the epoch of reionization (e.g., Zavala et al. 2021), finding evidence that relatively massive galaxies at these early epochs have significant levels of dust attenuation is not surprising in and of itself, as there are a growing number of direct individual detections of dust emission at z ≳ 8 (e.g., Watson et al. 2015; Laporte et al. 2017; Tamura et al. 2019; Bakx et al. 2020; Fudamoto et al. 2021; Schouws et al. 2021). A number of theoretical models have explored these results, finding that with a variety of assumptions the implied dust masses at early times could be formed with our current understanding of dust formation physics (e.g., Mancini et al. 2015, 2016; Popping et al. 2017; Behrens et al. 2018; Graziani et al. 2020; Sommovigo et al. 2020).

An important caveat to highlight is that the connection between dust attenuation and the physical properties of the dust in galaxies (such as the dust mass and the grain properties) is nontrivial. Neglecting the effect of geometry and orientation on attenuation can severely bias the interpretation (e.g., Padilla & Strauss 2008). For example, Chevallard et al. (2013) show that geometry and orientation effects have a stronger influence on the shape of the attenuation curve than changes in the optical properties of dust grains. Similarly, several studies show that galaxy shape and inclination are the major factors in determining the observed amount of dust attenuation, and not the galaxy dust mass (Maller et al. 2009; Kreckel et al. 2013; Zuckerman et al. 2021). Although these studies focus on lower-redshift systems (z < 3), similar effects might drive some of the observed effects we see at z > 4 regarding β, A_V, and the attenuation curve (Section 3.3.2). Although parts of this caveat can be alleviated by including far-IR constraints (modulo the assumption regarding energy conservation), this should be kept in mind in the following paragraphs when connecting the attenuation to the physical properties of dust.

The young age of the universe at these observed epochs could in principle constrain the efficiencies of different dust production mechanisms. The formation of dust grains can happen via multiple sources, which have their own timescales, and uncertainties due to various physical assumptions (see, e.g., Dayal & Ferrara 2018 and references therein). For example, while dust formation in the ejecta of supernovae (SNe) could lead to the formation of the first dust grains at extremely early times (e.g., Todini & Ferrara 2001; Schneider et al. 2004; Bianchi & Schneider 2007; Sarangi & Cherchneff 2015; Sluder et al. 2016; Marassi et al. 2019), the dust destruction timescales are not well constrained (e.g,. Bianchi & Schneider 2007; Silvia et al. 2010; Bocchio et al. 2016; Micelotta et al. 2016; Martínez-González et al. 2019; Slavin et al. 2020), especially in the early universe (e.g., Hu et al. 2019).

Dust can also form in the atmospheres of asymptotic giant branch (AGB) stars (e.g., Gehrz 1989; Ferrarotti & Gail 2006; Zhukovska et al. 2008; Ventura et al. 2012; Nanni et al. 2013; Ventura et al. 2018; Dell'Agli et al. 2019), with yields being sensitive to the mass and metallicity of their progenitor stars. Depending on the SFH and on the metallicity of the stellar population, Valiante et al. (2009) found that AGB stars could plausibly contribute to 30%–50% of the total dust budget in high-redshift galaxies in ≈300 Myr. Finally, dust can grow in the cold/warm phase of the ISM on seed grains formed by early SNe (e.g., Draine & Salpeter 1979; Draine 2009; Hirashita & Kamaya 2000; Michałowski et al. 2010; Valiante et al. 2011; Asano et al. 2013; Mancini et al. 2015; Leśniewska & Michałowski 2019). The timescale for this process may be quite short if dust is formed via the first SNe, thus this formation pathway may be significant at early times. However, we still lack a full understanding of this process at the atomic level, and we equally do not know the phase of the ISM where the process may occur, e.g., molecular (Ferrara et al. 2016; Ceccarelli et al. 2018) versus warm atomic (Zhukovska et al.2018).

As the grain growth timescale is thought to be density dependent (Asano et al. 2013; Schneider et al. 2014; Mancini et al. 2015; Popping et al. 2017), the expected higher density of star-forming clouds in these early galaxies could lead to this mode of dust production being more efficient. As ISM grain growth also requires initial seed grains, more efficient grain growth at earlier times could point to more efficient dust production by core-collapse SN explosions from low-metallicity massive stars (e.g., Marassi et al. 2015, 2019) or pair-instability SN explosions from Population III stars (Nozawa et al. 2003; Schneider et al. 2004), due to either higher yields (e.g., Schneider et al. 2004) or earlier Population III star formation times (e.g., Jaacks et al. 2018, 2019). However, these seed grains require some chemical enrichment to form, so this entire process is also dependent on the metallicity of the gas.

Graziani et al. (2020) explored dust formation at high redshift by including dust formation and evolution in a hydrodynamic simulation, accounting for dust formation via both stellar sources (e.g., SNe and AGB stars) and grain growth in the ISM, and several sources of dust destruction. They found that, at z > 10, dust produced via stellar sources dominates the total cosmic dust mass, with grain growth not playing a significant role until z < 9. This in principle might predict that massive galaxies at z > 9 should begin to appear significantly less dusty as they are not yet enriched via grain growth, seemingly at odds with our observation. However, they also found that grain growth becomes dominant for systems with stellar masses of $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 8.5$, consistent with the mass range of our observed reddened galaxies. It is thus plausible that even at these early epochs, these massive systems have their total dust content enriched via stellar dust production and grain growth, maybe aided by more favorable conditions in their ISM (e.g., higher densities), which is also consistent with the predictions from semianalytical (Popping et al. 2017; Vijayan et al. 2019) and seminumerical (e.g., Mancini et al. 2015, 2016) models.

We compare our observations to predictions from a semianalytic model (SAM; Yung et al. 2019b), the FLARES simulations (Lovell et al. 2021; Vijayan et al. 2021), and the THESAN radiation-magneto-hydrodynamic simulation (Garaldi et al. 2021; Smith et al. 2021; Kannan et al. 2022) to our observations in Figure 8. The β values for the THESAN simulation are presented in Kannan et al. (2021). While the SAM does track dust formation and destruction, it does not have any galaxies as massive as those we observe at z > 9. FLARES and THESAN do not directly track dust, rather both models use the ISM metal abundance to derive a dust attenuation, which should be kept in mind when comparing these models to our observations. It is interesting however that the FLARES predictions agree well with our observations. While the SAM seems to overpredict our observations, we also need to account for our observational bias. If we apply our observational cut of H < 26.6 to the SAM galaxies, we find predicted β values well in agreement with our observations (darker red dashed line). This model thus would predict a population of more dusty high-redshift massive systems missed by our UV-selection, similar to those recently discovered by Fudamoto et al. (2021). Galaxy selection at redder wavelengths, as will soon be possible with JWST, will alleviate this potential selection bias.

While our observations cannot alone distinguish between the various competing physical processes, they do point to fairly efficient dust growth in massive galaxies at early times, which can in turn be used to further constrain detailed simulations (e.g., McKinnon et al. 2017; Aoyama et al. 2018; Graziani et al. 2020; Vogelsberger et al. 2020). The abundance of JWST Cycle 1 programs targeting the early universe should both allow measures of the rest-UV colors of larger samples of massive galaxies at z ∼ 9–10, as well as push to lower-mass systems at z > 10 for the first time. Together with radiative transfer simulations (e.g., Behrens et al. 2018; Katz et al. 2019a; Shen et al. 2020, 2022; Vijayan et al. 2022), this will allow a detailed, more direct comparison between theoretical dust models and observations over a wide range of different galaxies, and thereby constrain the physical processes related to dust growth and chemical enrichment in early galaxies.

We present in this section the stellar mass (M_⋆) and SFR measurements of our z = 9–11 galaxy candidates. If the SFR varies with time, it is important to specify the timescale over which the SFR is measured (e.g., Caplar & Tacchella 2019). We choose a timescale of 50 Myr, which is roughly the timescale that the UV light at 1500–3000 Å probes. We label this SFR as SFR₅₀. As the galaxies at these early cosmic times are young, i.e., it is plausible that all the stellar mass of a galaxy has formed within the past 50 Myr, it is useful to consider what this maximal SFR is given the stellar mass M_⋆ (Tacchella et al. 2018a):

$$\begin{eqnarray}&&{\mathrm{SFR}}_{\max }=\displaystyle \frac{{M}_{\star }}{{t}_{\mathrm{SF}}},\end{eqnarray} \tag{ 2 }$$

where t_SF = 50 Myr is the timescale of the SFR indicator. As an example, for a M_⋆ = 10¹⁰ M_⊙ galaxy, the maximum SFR is ${\mathrm{SFR}}_{\max }=200\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Similarly, the maximum specific SFR (sSFR) is given by:

$$\begin{eqnarray}&&{\mathrm{sSFR}}_{\max }=\displaystyle \frac{{\mathrm{SFR}}_{\max }}{{M}_{\star }}=\displaystyle \frac{1}{{t}_{\mathrm{SF}}}.\end{eqnarray} \tag{ 3 }$$

This implies that the maximum sSFR is independent of mass (and cosmic epoch) and only depends on the SFR timescales. In our case, the maximum sSFR is ${\mathrm{sSFR}}_{\max }=20\,{\mathrm{Gyr}}^{-1}$. A corollary is that when considering long SFR timescales (relative to the ages of the galaxies), a perfect correlation between the SFR and M_⋆ is introduced by construction—important to consider when studying the star-forming main sequence.

After these general considerations, we plot the inferred SFRs and M_⋆ in Figure 10. The left and right panels of Figure 10 show the SFR₅₀–M_⋆ and the sSFR₅₀–z planes, respectively. The black lines and the gray shaded regions indicate the maximum SFR and sSFR, respectively, mentioned above. The red data points and error bars show the median and 16–84th percentiles of our fiducial run with the continuity prior for the SFH. The exact values are also given in Table 3. The smaller blue, purple, and green data points indicate the results of the bursty continuity prior, the Dirichlet prior and the parametric prior, respectively. The circle symbols show the galaxies with photometric redshift estimates, while the squares mark the objects with a spectroscopic redshift.

Zoom In Zoom Out Reset image size

Despite the large uncertainty in the individual measurements, we find that more massive galaxies have a higher SFR (left panel of Figure 10). The slope of the M_⋆–SFR relation—estimated with the orthogonal distance regression (ODR) taking into account the uncertainties in M_⋆ and SFR of each galaxy—is 0.70 ± 0.17 for our fiducial continuity prior, i.e., the higher mass galaxies have typically a lower sSFR than their lower mass counterparts. Although our slope estimates include the propagation of the errors of the inferred M_⋆ and SFR via the ODR, we do not perform a fully hierarchical Bayesian approach to measure the slope, intercept, and scatter of the main sequence (Curtis-Lake et al. 2021). Furthermore, the exact value of this slope depends on the SFH prior: the bursty continuity prior typically leads to a decrease in M_⋆ and an increase in the SFR measurements for the low-mass galaxies, which flattens the M_⋆–SFR relation. Specifically, the bursty continuity prior results in a slope of 0.45 ± 0.27, while the Dirichlet prior and the parametric prior lead to a slope of 0.79 ± 0.16 and 0.80 ± 0.12, respectively.

Measurements of the star-forming main sequence slope have been mainly published at slightly lower redshifts. We focus here on the stellar mass range of 10⁹–10¹¹ M_⊙ at z > 4, where the "bending" of the star-forming main sequence plays presumably a minor role. Salmon et al. (2015) found a rather shallow slope of 0.54 ± 0.16 at z ∼ 6, while Pearson et al. (2018) and Khusanova et al. (2021) inferred a slope of 1.00 ± 0.22 and 0.66 ± 0.21 at z ∼ 5.5, respectively. Based on a large literature compilation of z < 7 studies, Speagle et al. (2014) inferred a steeper slope of 0.84 ± 0.02 − (0.026 ± 0.003) × t, where t is the age of the universe in Gyr (at z = 10, this inferred slope is 0.83 ± 0.02). All of these estimates are consistent with our estimate when considering the uncertainty. Theoretical models typically produce steeper slopes, closer to 1 (Somerville et al. 2015; Tacchella et al. 2018a; Behroozi et al. 2019; Yung et al. 2019b). A more careful comparison between observations and theory of the star-forming main sequence slope (and in particular its scatter) will be useful to shed more light onto the star formation efficiency in low-mass halos (e.g., Tacchella et al. 2020) and the underlying assembly of dark matter halos (e.g., Dayal et al. 2015; Khimey et al. 2021).

The right panel of Figure 10 shows the redshift evolution of the sSFR of galaxies with masses of M_⋆ ≈ 10⁹–10¹⁰ M_⊙. We measure sSFR values in the range of 3–10 Gyr⁻¹, which indicate a mass-doubling timescale of ∼100–300 Myr under the assumption of a constant sSFR, roughly consistent with our age estimate (Section 5.3). An exception is UDS-18697 for which we measure a low sSFR of ${0.4}_{-0.3}^{+0.6}\,{\mathrm{Gyr}}^{-1}$—an interesting galaxy that seems to have gone through an intense episode of star formation early on and now is showing a declining SFH and a high stellar mass (Section 5.2). An important caveat to this object is that the IRAC deblending uncertainty is large (Appendix A). Despite this uncertainty and even if this object does not have such a low sSFR, it viscerally shows that your prior does not exclude such "old" solutions if the data warrant it. For the bursty continuity SFH prior, the sSFR values are typically larger and in some cases reach the maximum sSFR of 20 Gyr⁻¹.

We have also added to the right panel of Figure 10 the lower-redshift measurements from González et al. (2014), Stark et al. (2013), Salmon et al. (2015), and Stefanon et al. (2021), and the predictions from the models by Yung et al. (2019b), Behroozi et al. (2019), and Tacchella et al. (2018a). The Behroozi et al. (2019) and Tacchella et al. (2018a) models both track the evolution of the dynamical timescale of dark matter halos, while the model by Yung et al. (2019b) predicts a steeper increase with the redshift. Our fiducial sSFR values are consistent with the lower redshift estimates and lie slightly below the expected, increasing evolution of the theoretical models. This, however, depends on the assumed prior: the burstier SFH prior leads to higher sSFR, slightly higher—but still consistent—with theoretical expectations.

Figure 11 shows the detailed posterior distribution for SFR₅₀ and M_⋆ for all the SFH priors. Each panel shows an individual galaxy. It is difficult to draw a single conclusion from this figure, as each galaxy seems to show a different dependence on varying the prior. A common and important feature is that the posteriors of the different priors overlap, i.e., the resulting posteriors are consistent with each other. Furthermore, the typical differences in the median values of the SFR and M_⋆ measurements are of the order of a factor of 2–3 (except UDS-18697, for which the SFR varies by over 3 orders of magnitude). The parametric SFH prior typically leads to lower masses (as the ages are younger), consistent with findings at lower redshifts (Leja et al. 2019b; Tacchella et al. 2022). The bursty continuity SFH prior shows the largest spread in the posterior space, in particular along the SFR axis.

Zoom In Zoom Out Reset image size

As highlighted in the Introduction and Section 4.1, the UV contains plenty of information regarding the age of the stellar population. The key challenge is to differentiate age-related effects from other effects, such as dust attenuation or metallicity variations. Therefore, we have chosen to use a rather flexible SED model (Section 3.1), which includes a variable dust attenuation law and a range of different SFH priors (Section 3.2). Here we focus on the inferred SFHs and their dependence on the prior, while the next section focuses on the degeneracy between the stellar age and the attenuation.

Figure 12 shows the inferred SFHs by plotting the fraction of the mass formed, f_M , as a function of lookback time. Each panel shows an individual galaxy, with the red, blue, purple, and green lines showing the median SFHs obtained from the continuity, the bursty continuity, the Dirichlet, and the parametric SFH priors, respectively. The shaded regions show the 16–84th percentiles. The different SFH priors result in different SFH posteriors, i.e., it is important to fully understand how the inferred SFHs are affected by the choice of the prior.

Zoom In Zoom Out Reset image size

Figure 12 emphasizes when most of the mass has formed, and the SFH stresses that the mass heavily depends on the assumed SFH prior. The linear increase in f_M found when adopting the continuity and the Dirichlet prior is consistent with a constant SFH. Both the parametric and the bursty continuity SFH priors imply a stronger increase in f_M in more recent times relative to the fiducial continuity prior. A nominal exception is UDS-18697, where all the priors consistently find an early burst of star formation and little mass growth in the past ∼100–200 Myr. In Appendix B, we also show the SFHs plotted as SFR as a function of time (Figure 20). There, it is more clearly visible that a few galaxies (i.e., COSMOS-20646, EGS-68560, GOODSN-35589, and UDS-18697) show significant variation in the past ∼100 Myr.

Importantly, the otherwise rather constant behavior for the fiducial continuity and Dirichlet priors is expected, as the expectation value of these two priors are a constant SFR(t) (Section 3.2 and Figure 1). For the parametric SFH prior, we find for all except one galaxy (UDS-18697) an increasing SFH, something we expect from theoretical models (e.g., Tacchella et al. 2018a). However, again, this is the expected behavior of the prior. This underscores the worry that the current data provide little constraining power when it comes to the SFH. The different priors, all of which provide equally good fits to the data, are producing rather different SFH posteriors.

Figure 13 plots the mass- and redshift-dependence of the stellar ages (t₅₀) and star formation timescales (τ_SF). The stellar age is the lookback time at which 50% of the final mass is formed. The star formation timescale is defined as the time it takes to increase the stellar mass from 20% to 80% of the final mass, i.e., it is a measure of how quickly a galaxy formed the bulk of its stellar mass. Figure 13 shows t₅₀ on top and τ_SF on the bottom. The colors and symbols are the same as those in Figure 10.

Zoom In Zoom Out Reset image size

Independent of the SFH prior, we find that more massive galaxies have typically older ages. There is also a hint that galaxies at higher redshifts have younger ages, though the scatter and uncertainties are large. Parts of this trend can probably be explained by an "envelope effect", where the galaxy age and the onset of star formation are required to be less than the age of the universe. The star formation timescales do not show a convincing trend with the stellar mass. However, the trend with the cosmic epoch is more pronounced and roughly follows the age of the universe: these galaxies form their stars on timescales of roughly 30% to 40% of the age of the universe, as indicated by the shaded band in the lower-right panel of Figure 13. The exact value of the ages and the star formation timescales depend on the SFH prior: the bursty continuity prior results in roughly 3 times younger ages (and shorter star formation timescales) than the continuity prior. The exact values of the ages and star formation timescales are given in Table 4. We find a median age t₅₀ over the whole sample of ∼150 Myr for the continuity and Dirichlet priors, while the median age is only ∼50 Myr for the bursty continuity prior. The parametric SFH prior produces ages that lie in-between with a median age of ∼90 Myr. This means we cannot distinguish between median ages of ∼50–150 Myr, which corresponds to 50% formation redshifts of z₅₀ ∼ 10–12. These ages are consistent with expectations from the semianalytical model by Yung et al. (2019b) and the empirical halo model by Tacchella et al. (2018a), which find stellar ages t₅₀ for galaxies at z ≈ 10 of 60 Myr and 40 Myr, respectively.

Table 4. Stellar Age (t₅₀), Star Formation Timescale (τ_SF), and Fraction of Mass Formed before z = 10 (f_M ) for our Galaxy Candidates at z = 9–11

	Continuity Prior				Bursty Continuity Prior				Dirichlet Prior				Parametric Prior
ID	t50	τSF	fM		t50	τSF	fM		t50	τSF	fM		t50	τSF	fM
	(Myr)	(Myr)	(%)		(Myr)	(Myr)	(%)		(Myr)	(Myr)	(%)		(Myr)	(Myr)	(%)
COSMOS-20646	${137}_{-75}^{+51}$	${170}_{-45}^{+28}$	${94}_{-16}^{+5}$		${111}_{-49}^{+98}$	${70}_{-47}^{+88}$	${99}_{-10}^{+1}$		${165}_{-53}^{+45}$	${160}_{-47}^{+41}$	${96}_{-4}^{+2}$		${107}_{-47}^{+44}$	${148}_{-63}^{+54}$	${92}_{-6}^{+3}$
COSMOS-47074	${121}_{-73}^{+59}$	${168}_{-60}^{+32}$	${95}_{-11}^{+4}$		${27}_{-15}^{+57}$	${18}_{-12}^{+100}$	${53}_{-51}^{+45}$		${146}_{-58}^{+55}$	${167}_{-58}^{+39}$	${94}_{-6}^{+3}$		${79}_{-49}^{+56}$	${114}_{-69}^{+74}$	${97}_{-6}^{+1}$
EGS-6811	${186}_{-77}^{+57}$	${215}_{-50}^{+34}$	${74}_{-21}^{+13}$		${104}_{-77}^{+136}$	${59}_{-45}^{+111}$	${51}_{-51}^{+39}$		${193}_{-74}^{+61}$	${203}_{-64}^{+50}$	${75}_{-19}^{+12}$		${81}_{-52}^{+72}$	${119}_{-77}^{+90}$	${40}_{-40}^{+26}$
EGS-44164	${184}_{-84}^{+57}$	${216}_{-57}^{+28}$	${73}_{-23}^{+13}$		${61}_{-42}^{+88}$	${35}_{-25}^{+108}$	${11}_{-10}^{+59}$		${184}_{-64}^{+69}$	${204}_{-67}^{+48}$	${74}_{-18}^{+11}$		${100}_{-49}^{+61}$	${141}_{-66}^{+79}$	${49}_{-34}^{+18}$
EGS-68560	${80}_{-48}^{+72}$	${188}_{-124}^{+44}$	${60}_{-26}^{+16}$		${23}_{-10}^{+15}$	${10}_{-5}^{+28}$	${1}_{-1}^{+12}$		${127}_{-53}^{+69}$	${195}_{-76}^{+44}$	${71}_{-16}^{+11}$		${54}_{-31}^{+48}$	${81}_{-46}^{+66}$	${49}_{-41}^{+22}$
EGS-20381	${183}_{-95}^{+71}$	${225}_{-68}^{+44}$	${68}_{-25}^{+15}$		${63}_{-46}^{+119}$	${43}_{-34}^{+128}$	${11}_{-11}^{+63}$		${179}_{-59}^{+76}$	${209}_{-70}^{+49}$	${69}_{-18}^{+13}$		${97}_{-59}^{+61}$	${138}_{-82}^{+82}$	${42}_{-42}^{+21}$
EGS-26890	${141}_{-75}^{+65}$	${193}_{-67}^{+34}$	${74}_{-24}^{+15}$		${42}_{-26}^{+63}$	${30}_{-23}^{+91}$	${14}_{-14}^{+53}$		${163}_{-57}^{+69}$	${189}_{-62}^{+49}$	${78}_{-17}^{+11}$		${97}_{-51}^{+58}$	${139}_{-72}^{+76}$	${64}_{-32}^{+13}$
EGS-26816	${128}_{-82}^{+70}$	${178}_{-77}^{+38}$	${81}_{-28}^{+14}$		${29}_{-17}^{+80}$	${22}_{-15}^{+85}$	${24}_{-23}^{+73}$		${145}_{-56}^{+68}$	${176}_{-62}^{+43}$	${84}_{-11}^{+8}$		${84}_{-51}^{+56}$	${120}_{-70}^{+74}$	${75}_{-33}^{+10}$
EGS-40898	${178}_{-107}^{+76}$	${220}_{-75}^{+55}$	${65}_{-29}^{+17}$		${74}_{-59}^{+183}$	${56}_{-46}^{+103}$	${26}_{-26}^{+67}$		${183}_{-68}^{+72}$	${209}_{-71}^{+53}$	${68}_{-19}^{+14}$		${93}_{-57}^{+70}$	${132}_{-79}^{+90}$	${39}_{-39}^{+24}$
GOODSN-35589	${73}_{-46}^{+48}$	${130}_{-56}^{+27}$	${100}_{-0}^{+0}$		${11}_{-6}^{+97}$	${14}_{-8}^{+74}$	${100}_{-0}^{+0}$		${99}_{-40}^{+49}$	${131}_{-48}^{+30}$	${100}_{-0}^{+0}$		${49}_{-32}^{+49}$	${73}_{-48}^{+65}$	${100}_{-0}^{+0}$
UDS-18697	${225}_{-22}^{+13}$	${121}_{-16}^{+32}$	${99}_{-1}^{+0}$		${262}_{-60}^{+19}$	${64}_{-11}^{+40}$	${100}_{-0}^{+0}$		${207}_{-31}^{+11}$	${116}_{-16}^{+38}$	${97}_{-2}^{+1}$		${218}_{-54}^{+41}$	${205}_{-23}^{+18}$	${96}_{-3}^{+2}$

Note. We list the results for the four different SFH priors (Section 3.2).

Download table as:  ASCIITypeset image

In Appendix C, Figures 21 and 22 show the posterior distribution of M_⋆ versus t₅₀ and UV attenuation A_UV versus t₅₀, respectively. We show that there is a degeneracy between M_⋆ and t₅₀: a younger age implies a lower stellar mass. This is expected as younger stellar populations are typically brighter at fixed stellar mass. Additionally, we also demonstrate that it is challenging to break the dust-age degeneracy with our current observational data. Specifically, we find an anticorrelation between older ages and more UV attenuation. The UV attenuation is overall not well constrained, i.e., we find rather wide posteriors with uncertainties of more than 1 mag, which can at least in part be explained by the degeneracy with the attenuation law.

Previously, Stefanon et al. (2019) measured stellar ages of ${22}_{-22}^{+69}\,\mathrm{Myr}$ (median and 68% confidence interval over their sample) for 18 bright z = 8 galaxies (H_AB < 25) with stellar masses of M_⋆ ∼ 10⁹ M_⊙. They used the FAST code (Kriek et al. 2009) and assumed a fixed subsolar 0.2 Z_⊙ metallicity, a constant SFH, and a fixed Calzetti et al. (2000) attenuation law. Our 10⁹ M_⊙ galaxies are slightly older (50–100 Myr) when considering the fiducial continuity prior (and therefore an SFH that is similar to a constant as Stefanon et al. 2019 assumed), which might imply that fixing the attenuation law leads to this difference as their typical attenuation values are also lower than ours (${A}_{{\rm{V}}}={0.15}_{-0.15}^{+0.30}\,\mathrm{mag}$).

As described in the Introduction, Laporte et al. (2021) (see also Hashimoto et al. 2018) study six z ∼ 9–10 galaxies (out of which three have a spectroscopic redshift) and perform the SED modeling with BAGPIPES (Carnall et al. 2018), assuming a fixed Calzetti et al. (2000) attenuation law and model emission lines self-consistently. They investigate four different SFH prescriptions (delayed, exponential, constant, or burst-like), where the best fit is always obtained for the delayed or constant SFH. They quote ages of 200–500 Myr, but those ages do not correspond to our t₅₀ (or mass-weighted ages), but to the time since star formation has started. This could partially explain why these ages are significantly (factor of 10) older than what Stefanon et al. (2019) found and also a factor of 2–3 older than what we find with our fiducial continuity prior (∼150 Myr). Indeed, when computing t₂₀ (lookback age at which 20% of the stellar mass was formed), which is a better tracer of the first star formation episode, we find a median over the whole sample of ${260}_{-105}^{+59}\,\mathrm{Myr}$ and ${265}_{-90}^{+58}\,\mathrm{Myr}$ for the continuity and Dirichlet priors, respectively. These estimates are consistent with the ones inferred by Laporte et al. (2021). However, we still infer younger t₂₀ when adopting the bursty continuity (${91}_{-72}^{+183}\,\mathrm{Myr}$) or parametric (${163}_{-99}^{+117}\,\mathrm{Myr}$) priors. An exception is UDS-18697, where the results from all priors are consistent with a very early buildup of stellar mass: most of the mass formed around z ≈ 15 and the mass fraction formed before z > 12 is >60% (Figure 23). The reason for this is the strong Balmer break (see Figure 2, but note that we question the reliability of the IRAC photometry of this object), which can only be explained by older stellar populations. This indicates that these age differences between our work and the work by Laporte et al. (2021) might not only be caused by the different methodologies and definitions, but also because of sample selection. Indeed, the selection of the sample in different works differs significantly: Laporte et al. (2021) only selected objects with red IRAC [3.6]−[4.5] colors at z > 9 (and detection in both bands) to specifically search for older stellar populations, while we do not employ such a cut. Importantly, we do not claim that we can rule out old ages; we claim that our current data cannot unambiguously confirm old ages in all galaxies in our sample.

We now focus on the fraction of stellar mass formed at early times. Figure 14 shows the posterior distribution of the fraction of mass formed before redshift 10 (f_M (z > 10)) and stellar age t₅₀. The inferred values of f_M (z > 10) and their uncertainties are also given in Table 4. By definition, f_M (z > 10) is only a useful number if the redshift of the galaxy is below z = 10, i.e., f_M (z > 10) = 100% for galaxy GOODSN-35589 with z_spec = 10.96. Therefore, Figure 23 in Appendix D shows f_M (z > 12).

Zoom In Zoom Out Reset image size

Figure 14 makes the point that the uncertainty in stellar age directly translates into an uncertainty in f_M (z > 10). Therefore, different SFH priors can lead to substantially different estimates of f_M (z > 10). A good example of this is COSMOS-47074, which has f_M (z > 10) > 90% when considering the continuity prior, the Dirichlet prior, or the parametric prior, while the fraction is only ${f}_{M}(z\gt 10)={53}_{-51}^{+45} \%$, albeit a large uncertainty when adopting the bursty continuity prior. Across our full sample, we find a median (and 68% percentile) for f_M (z > 10) of ${79}_{-28}^{+18} \%$, ${39}_{-39}^{+61} \%$, ${82}_{-21}^{+14} \%$, and ${65}_{-45}^{+31} \%$ when adopting the continuity prior, the bursty continuity prior, the Dirichlet prior, and the parametric prior, respectively. In summary, we conclude that the fraction of stellar mass formed at z > 10 (and also z > 12; see Appendix D) is not well constrained by our data.

An interesting application of the derived SFHs is to study the implications for the early mass assembly of stellar mass as this provides an independent insight into the evolution of the UV luminosity density within the first ∼500 Myr. In the literature, it is debated whether the luminosity density over the redshift range 8 < z < 11 declines rapidly with ∝ (1 + z)⁻¹¹ (e.g., Oesch et al. 2014, 2018; Bouwens et al. 2015) or slowly with ∝ (1 + z)⁻⁴ (e.g., Finkelstein 2016; McLeod et al. 2016). From a theory perspective, a constant star formation efficiency model together with the buildup of dark matter halos can reproduce the suggested rapid increase in the cosmic SFR density with time (Tacchella et al. 2013, 2018a; Mason et al. 2015; Mashian et al. 2016; Yung et al. 2019a). However, there are also other models that prefer a rather slow increase in the cosmic SFR density at early times (Moster et al. 2018; Behroozi et al. 2020), consistent with SFHs for the six galaxies studied by Laporte et al. (2021).

Following the same approach, we average all our 10 f_M (t) SFHs (excluding GOODSN-35589 which lies at z_spec = 10.96) and plot their stacked mass assembly history f_M (z) in Figure 15. The important assumption when doing this is that the galaxy sample is representative of the overall galaxy population at this epoch. The red, blue, purple, and green lines show the resulting f_M (z) assuming the continuity, the bursty continuity, the Dirichlet, and the parametric SFH priors, respectively. For reference, the dotted, dashed–dotted, and dashed black lines indicate f_M (z) for a constant SFR, the rapidly increase cosmic SFR density (∝(1 + z)⁻¹¹), and a slowly increase cosmic SFR density (∝(1 + z)⁻⁴), respectively.

Zoom In Zoom Out Reset image size

As expected from the previous sections, our conclusion depends on the adopted SFH prior. If we assume the continuity or the Dirichlet prior, which are overall consistent with a constant SFH, we find they match well with the slowly declining SFR density. This is also compatible with the conclusions of Laporte et al. (2021), which is not surprising, as their best fit is a constant or delayed SFH model. On the other hand, if we adopt the bursty continuity prior or the parametric prior, we find f_M (z) increases more rapidly at more recent times, which is more consistent with a rapid decline in the cosmic SFR density at z > 10. We find qualitatively the same results if we include GOODSN-35589 and perform the analysis at z = 12. In summary, the presently available observational data cannot constrain the SFHs of our ensemble of galaxies well enough to distinguish between a rapid or more smooth decline in the cosmic SFR density at z > 10; therefore the epoch of first galaxy formation remains to be identified. Interestingly, there are indications (even when varying the prior; see in particular UDS-18697) that at least some galaxies have formed significant amounts of stellar mass by z ≈ 15 (Figure 12), which is of interest for the interpretation of the 21 cm signal and the formation and evolution of black holes in the early universe (e.g., Dayal & Ferrara 2018).

As discussed above, our stellar mass estimates are overall robust, including the variation of the SFH prior, but some sources might suffer from IRAC systematics as highlighted below. We now investigate whether these stellar masses are above the expectation of the Planck ΛCDM universe. In particular, given our survey volume of roughly 10⁶ Mpc³, high-redshift galaxy stellar masses can place interesting limits on number densities of massive halos, which itself constrains cosmology (Steinhardt et al. 2016; Behroozi & Silk 2018).

Figure 16 shows our stellar mass estimates as a function of redshift. The red and blue symbols show the measurements adopting the continuity and the bursty continuity SFH priors, respectively. The red and blue solid lines show the mass growth tracks of individual galaxies as inferred from the SFHs presented in Section 5.2. The black solid line marks the threshold stellar mass for a cumulative number density of Φ = 10⁻⁶ Mpc⁻³, which we adopt from Behroozi & Silk (2018). The assumption for this stellar mass threshold is a 100% star formation efficiency, i.e., the halo mass is related to the stellar mass via M_h = M_⋆/f_b, where f_b = 0.16 is the cosmic baryon fraction. This can be regarded as the maximal stellar mass as the average SFHs of galaxies inferred from the present-day stellar-to-halo mass relation are much less than the cosmic baryon fraction (e.g., Moster et al. 2018; Behroozi et al. 2019).

Zoom In Zoom Out Reset image size

We find that three galaxies (EGS-6811, COSMOS-20646, and UDS-18697) lie above the black line and in the gray shaded region that is less likely in the Planck ΛCDM universe. Galaxy EGS-6811 ($\mathrm{log}({M}_{\star }/{M}_{\odot })={10.6}_{-0.3}^{+0.2}$ with a spectroscopic redshift) actually lies on the threshold when considering the uncertainty. Therefore, this galaxy is not challenging ΛCDM. Interestingly, when studying its trajectory in the M_⋆–z plane, we find that the bursty continuity prior leads to a steep mass growth history, which means that it actually falls below the mass threshold by z ∼ 10–11, while it remains in the less likely ΛCDM region when adopting the continuity prior. An important note with respect to this statement is that these growth curves assume that there are no previous galaxy mergers involved in the mass growth. Or, equivalently, they represent the summed M_⋆ associated with all galaxies that merge to form the observed objects.

We have two galaxies (COSMOS-20646 and UDS-18697) that both lie at z ∼ 9–10 and have stellar masses of M_⋆ ≈10¹¹ M_⊙ and therefore lie solidly within the less likely ΛCDM region. However, we acknowledge that these two galaxies (along with EGS-6811) have bright neighbors in their proximity, making the IRAC photometry of these two sources the least reliable of our sample (see F21 and Section A). Thus it is possible that residual light from the neighbors is contributing to the high stellar mass measurement. Nevertheless, taking our fiducial stellar mass and redshift at face value (i.e., ignoring the systematic uncertainty in the IRAC photometry for the moment), we find that COSMOS-20646 and UDS-18697 lie within the less likely ΛCDM region at 3.0σ and 4.6σ significance.

However, as discussed in detail in (Behroozi & Silk 2018; see their Appendix A), the significance of this tension is not as sound as apparent on first sight. Attempts to rule out ΛCDM are limited by both cosmic variance and observational errors. Considering cosmic variance, as we have selected the galaxies from five different survey fields (Section 2), the chance actually significantly increases that one of the fields will have an "outlier" even in a standard ΛCDM universe (e.g., Trenti & Stiavelli 2008). The observational errors (considering an uncertainty of 0.2–0.3 dex in $\mathrm{log}\,{M}_{\star }$, but ignoring the larger systematic uncertainty stemming from the IRAC photometry) will inflate the number density of massive halos compared to the underlying true number density (see also the Eddington bias; Eddington 1913) because objects are preferentially scattered to higher masses when drawn from a steep mass function. Taking these effects into account, we calculate that the probability is ∼ 20% and ∼ 0.4% to find in our survey a galaxy as massive as COSMOS-20646 and UDS-18697, respectively.

Finally, another interesting insight from Figure 16 is about the SFH prior, as hinted at above. In the case of the continuity SFH prior, which typically leads to extended SFHs and older ages (Figure 13), the number of galaxies crossing the threshold increases toward higher redshifts. Contrarily, when considering the bursty continuity SFH prior, we find the opposite behavior, and galaxies depart from the less likely ΛCDM region. This implies that the continuity prior has a larger tendency to violate ΛCDM than the bursty continuity prior. As the uncertainties in the derived SFHs are large, we however cannot rule out any of the priors at the moment.

JWST will transform high-z galaxy evolution studies by providing near-IR (i.e., rest-frame optical) data of unprecedented depth, spatial, and spectral resolution. This will help to better constrain the rest-frame Balmer/4000 Å break and therefore get tighter constraints on the SFHs of z > 8 galaxies. A detailed discussion on the implications for these kinds of measurements is out of the scope of this paper, but see Roberts-Borsani et al. (2021) for an exploration of the improvement of z ∼ 7–11 galaxy property estimates with JWST/NIRCam medium-band photometry.

Here we focus on the implications of the different SFH priors regarding the JWST wavelength coverage. Specifically, we plot in Figure 17 the SED posterior for the four different SFH priors (left panel) and their log differences (middle panel). The orange and red lines on the bottom show the JWST/NIRCam broad- and medium-band filter curves. We focus here on EGS-44164 as this galaxy has a spectroscopic redshift and is detected with Spitzer/IRAC. As we can see from the left panel of Figure 17, the SEDs from the four different SFH priors are similar in the rest-frame UV wavelength range, while they start diverging in the rest-frame optical. We find strong emission lines but a weak 4000 Å continuum break for the parametric prior, while the break is stronger but the emission lines are nearly absent for the bursty continuity prior. The continuity and the Dirichlet prior both have a rather strong Balmer/4000 Å break and emission lines. These features can be directly understood by looking at the SFHs, which is shown as an inset in the figure. The error bars at the bottom right illustrate the 5σ uncertainty for ∼15 ksec exposures, estimated from the point source limit with a 01 aperture and medium background (Williams et al. 2021).

Zoom In Zoom Out Reset image size

The Spitzer/IRAC data cannot currently differentiate between these SEDs and SFHs. For JWST, thanks to its unprecedented sensitivity and its higher spectral resolution, progress can be made, in particular through the inclusion of emission lines. We show the F277W−F460M versus F460M−F444W color–color diagram in the right panel of Figure 17. We give the Hβ+[O iii] equivalent width (EW_{Hβ+[O III]}) distribution in the inset of the panel. The significant differences of the very recent (<20 Myr) SFHs between the four priors leads to contrasting EW distributions. For the continuity prior, the bursty continuity prior, the Dirichlet prior, and the parametric prior, we find ${\mathrm{EW}}_{{\rm{H}}\beta +[{\rm{O}}\,{\rm\small{III}}]}\,={675}_{-578}^{+2213}\,{\rm{\mathring{\rm A} }}$, ${\mathrm{EW}}_{{\rm{H}}\beta +[{\rm{O}}\,{\rm\small{III}}]}={33}_{-14}^{+512}\,{\rm{\mathring{\rm A} }}$, ${\mathrm{EW}}_{{\rm{H}}\beta +[{\rm{O}}\,{\rm\small{III}}]}={673}_{-463}^{+2361}\,{\rm{\mathring{\rm A} }}$, and ${\mathrm{EW}}_{{\rm{H}}\beta +[{\rm{O}}\,{\rm\small{III}}]}={3507}_{-2388}^{+4314}\,{\rm{\mathring{\rm A} }}$, respectively.

These different EW measurements are then directly reflected in the color–color diagram as the Hβ line straddles the medium-band F460M for this specific case. Hence, these red medium bands will help constrain the recent SFH and will generally provide more stringent constraints on the stellar populations (see also Roberts-Borsani et al. 2021). Obviously, having higher spectral resolution information (via for example the NIRSpec/Prism or NIRCam/Grism) will further constrain these emission lines. One caveat to this is the rather large uncertainty on the number of ionizing photons that power those emission lines, in particular related to the escape and dust absorption of those photons (e.g., Kimm et al. 2017; Smith et al. 2017; Glatzle et al. 2019) as well as the production efficiency (stellar binarity and rotation; e.g., Choi et al. 2017; Eldridge et al. 2017).