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

We present a detailed stellar population analysis of 11 bright ($H<26.6$) galaxies at $z=9-11$ (three spectroscopically confirmed) to constrain the chemical enrichment and growth of stellar mass of early galaxies. We use the flexible Bayesian spectral energy distribution (SED) fitting code Prospector with a range of star-formation histories (SFHs), a flexible dust attenuation law and a self-consistent modeling of emission lines. This approach allows us to assess how different priors affect our results, and how well we can break degeneracies between dust attenuation, stellar ages, metallicity and emission lines using data which probe only the rest-frame ultraviolet to optical wavelengths. We measure a median observed ultraviolet spectral slope $\beta=-1.87_{-0.43}^{+0.35}$ for relatively massive star-forming galaxies ($9<\log(M_{\star}/M_{\odot})<10$), consistent with no change from $z=4$ to $z=9-10$ at these stellar masses, implying rapid enrichment. Our SED-fitting results are consistent with a star-forming main sequence with sub-linear slope ($0.7\pm0.2$) and specific star-formation rates of $3-10~\mathrm{Gyr}^{-1}$. However, the stellar ages and SFHs are less well constrained. Using different SFH priors, we cannot distinguish between median mass-weighted ages of $\sim50-150$ Myr, which corresponds to 50\% formation redshifts of $z_{50}\sim10-12$ at $z\sim9$ and is of the order of the dynamical timescales of these systems. Importantly, the models with different SFH priors are able to fit the data equally well. We conclude that the current observational data cannot tightly constrain the mass-buildup timescales of these $z=9-11$ galaxies, with our results consistent with SFHs implying both a shallow and steep increase of the cosmic SFR density with time at $z>10$.


INTRODUCTION
The past decade has seen observational studies leap into the epoch of reionization, the time in the early universe when energetic photons (presumably from early star formation) ionized the gas in the intergalactic medium (IGM). Advances in near-IR imaging both in space with the Hubble Space Telescope (HST) and from the ground (with, e.g., Subaru and VISTA) have allowed the discovery of large samples of dropout galaxy candidates at 6 < z < 11 (e.g., Oesch et al. 2010;Ellis et al. 2013;Oesch et al. 2018;Finkelstein et al. 2015bFinkelstein et al. , 2021Bouwens et al. 2015Bouwens et al. , 2021bBowler et al. 2015Bowler et al. , 2017McLeod et al. 2016;Livermore et al. 2017;Atek et al. 2018;Harikane et al. 2021). Studying the properties of these galaxies, which exist at a time less than 1 Gyr after the Big Bang, can provide key constraints on the buildup of both stellar mass and heavy elements in early galaxies. In particular, the stellar population of galaxies at the earliest probed cosmic times (z > 8) ought to supply crucial information on the formation of the first stars and galaxies.
For instance, the rest-frame ultraviolet (UV) colors of these early galaxies can inform us about the earliest phases of chemical enrichment. The UV color is sensitive to dust attenuation, stellar ages, and stellar metallicities (e.g., Wilkins et al. 2011). At these early times, dust attenuation is believed to dominate, though very low metallicities can result in extremely blue colors. Early results at z ∼ 7 found some evidence that the faintest galaxies in the Hubble Ultra Deep Field (m AB ∼ 29; stellar mass of log(M /M ) ∼ 7 − 8) had rest-frame UV colors consistent with essentially no metals (e.g., Bouwens et al. 2010;Finkelstein et al. 2010), though follow-up studies with larger samples accounting for selection biases (e.g., Dunlop et al. 2012) were able to rule out Population III-dominated galaxies at this epoch (Finkelstein et al. 2012;Dunlop 2013;Bouwens et al. 2014).
Looking at the full dynamic range of the galaxy population, correlations of the rest-UV color have been found with both the UV luminosity Stefanon et al. 2021) and stellar mass (Finkelstein et al. 2012;Bhatawdekar & Conselice 2021), where more luminous/massive systems have redder observed colors. In particular, Finkelstein et al. (2012) found that the most massive galaxies in their sample (log(M /M ) ∼ 9 − 10) had similarly red rest-UV colors from z = 4 − 7, indicating a roughly constant level of dust attenuation in these galaxies. Pushing these measurements to even earlier cosmic epochs can constrain exactly when dust began forming in the early universe, potentially constraining the respective efficiencies of different dust production mechanisms (e.g., Valiante et al. 2011Valiante et al. , 2014Mancini et al. 2015Mancini et al. , 2016Popping et al. 2017;Aoyama et al. 2018;Graziani et al. 2020).
Additionally, the rest-frame UV emission of these high-redshift galaxies, which can be currently probed with HST in the near-IR, contains a wealth of information regarding the ages of the stars: relatively young stars (∼ 10 7 yr) will dominate the observed emission in both the far-and near-UV rest frame, while older stars (∼ 10 9 yr) will contribute more to the near-UV than they do to the far-UV (e.g., Conroy 2013). As at early cosmic times the stellar populations are younger (with an upper limit given by the age of the universe), the rest-UV light can be used to infer stellar ages and stellar masses. A major challenge in using the UV as a tracer for the stellar age and the star-formation history in general is the degeneracy with other galaxy properties, including the wavelength-dependent attenuation of the UV emission by dust and the metallicity content of the stars (e.g., Papovich et al. 2001).
Extending the wavelength coverage further into the rest-frame optical is helpful to constrain the stellar populations and break some of these degeneracies. This can currently be done with Spitzer/IRAC for z > 8 galaxies (e.g., Stefanon et al. 2019), though it remains challenging due to low signal-to-noise and deblending issues. Furthermore, emission lines such as Hβ, [O III] and [O II] can contaminate the IRAC bands and therefore can be confused with a strong Balmer/4000Å break (Labbé et al. 2013;Finkelstein et al. 2013;Smit et al. 2015;Faisst et al. 2016;Roberts-Borsani et al. 2016;De Barros et al. 2019;Endsley et al. 2021).
Several studies have made use of combining HST with Spitzer/IRAC data in order to constrain the stellar populations of z > 8 galaxies. Stefanon et al. (2019) measured for 18 bright z = 8 galaxies an average stellar mass of M = 10 9.1 +0.5 −0.4 M , star-formation rate (SFR) of SFR = 32 +44 −32 M yr −1 , and stellar age of 22 +69 −22 Myr. At higher redshifts, MACS1149-JD1 at z spec = 9.11 gained a lot of attention due to its red IRAC color, which was attributed to old stellar populations (Zheng et al. 2012). In particular, Hashimoto et al. (2018) inferred an age of 290 ± 150 Myr by fitting a young and old stellar population to the data (see also Roberts-Borsani et al. 2020). Laporte et al. (2021) studied six z ∼ 9 galaxies selected to have 4.5µm flux excesses (out of which 3 have a spectroscopic redshift) and found stellar ages (here referring to the time since beginning of star formation) of 200−500 Myr (age of MACS1149-JD1 is consistent with 500 Myr), with the best fit being always obtained for a delayed or constant star-formation history (SFH).
Constraints on these ages provide our crucial first glimpse into the buildup of stellar mass at z > 10. One of the major systematic uncertainties on the ages of these galaxies from previous studies is the choice of the SFH. The derived ages are crucially dependent on this assumption (e.g., Papovich et al. 2011;Curtis-Lake et al. 2013;Schaerer et al. 2013;Buat et al. 2014;Leja et al. 2019b;Lower et al. 2020;Tacchella et al. 2021). This choice is also directly apparent when deriving the evolution of the cosmic SFR density from these ages and inferred SFHs. The cosmic SFR density is usually inferred from the UV luminosity function, with some studies suggesting the cosmic SFR declines with redshift more steeply at z > 8 than at 4 < z < 8 (e.g., Oesch et al. 2018;Bouwens et al. 2021a), while others suggest the evolution continues with a more shallow decline (e.g., McLeod et al. 2016;Finkelstein 2016). Laporte et al. (2021) investigated this by averaging the best-fit SFHs of their six galaxies, determining that these galaxies formed ∼ 70% of their mass by z = 10, which favors a smooth increase in the cosmic SFR density with time. However, the extent to which the priors on the assumed SFH and stellar population parameters impact the results have not yet been deeply investigated.
We present a new analysis of the properties of z ∼ 9 galaxies using a more flexible treatment of the SFH. We use a newly published sample of moderately bright z > 8.5 galaxies (Finkelstein et al. 2021, hereafter F21). These sources are selected in the HST CANDELS fields (Grogin et al. 2011;Koekemoer et al. 2011) and thus have an array of deep HST imaging available; yet they are also selected to be bright (m F160W < 26.6), allowing meaningful Spitzer/IRAC constraints on their restframe optical fluxes, which are crucial to constrain their stellar populations.
We perform a careful inference on the stellar populations by using Prospector, a flexible Bayesian spectral energy distribution (SED) fitting code (Johnson et al. 2021). In particular, we expand upon previous z > 6 SED investigations by adopting a range of simple and flexible models for the SFHs, a flexible dust attenuation law, self-consistent modeling of emission lines, and variable IGM absorption. We explore the dust reddening in these galaxies and thoroughly investigate how our inferred stellar ages depend on the adopted SFH prior. We conclude that our data are unable to meaningfully constrain the SFHs of these high-z galaxies, consistent with findings at lower redshifts (Strait et al. 2021). Specifically, the SFHs can be consistent with either a rapid or slow increase in the cosmic SFR density with time at z > 9.
The paper is structured as follows. Section 2 introduces the galaxy sample and its selection. Section 3 describes in detail the assumptions in our SED modeling. Section 4 discusses our key results concerning the chemical enrichment, while Section 5 focuses on the inferred growth of stellar mass and its implication on early cosmic star formation. We conclude in Section 6.
Throughout this work, all magnitudes are presented in the AB system, and we assume for the cosmological parameters H 0 = 67.74 km s −1 Mpc −1 , Ω M = 0.309 and Ω Λ = 0.691, consistent with the recent Planck Collaboration et al. (2020) measurements.

GALAXY SAMPLE
In this work, we study the sample of 11 bright (H < 26.6) galaxy candidates selected in the CANDELS fields by F21. Many of these sources were also presented in Oesch et al. (2018) and Bouwens et al. (2019). In F21 the authors created new photometric catalogs for each of the five CANDELS fields measuring accurate colors and total fluxes for all available HST imaging bands and obtained deblended photometry in the IRAC/Spitzer bands using TPHOT, following Song et al. (2016) and Merlin et al. (2016). Photometric redshifts were measured using EAZY (Brammer et al. 2008), using a large set of templates including a very blue template to match the expected colors of some high-redshift galaxies. Candidates were initially selected using a combination of criteria designed to select well-detected objects with photometric redshifts robustly constrained to be at z > 8. F21 noted that using Spitzer/IRAC in tandem with HST in the initial selection process more robustly removed potential contaminating systems but also resulted in tighter redshift constraints for likely high-redshift candidates.
This initial sample of galaxies was vetted in a variety of ways, including several screens against non-galactic sources (noise, persistence, stellar sources), the addition of ground-based photometric constraints, and follow-up HST imaging in additional filters. The final sample of 11 galaxies continued to satisfy all stringent criteria for a likely high-redshift nature. We list these 11 sources in Table 1, and make use of the photometry published in the tables in F21. We note that three of these sources are spectroscopically confirmed, as noted in Table 1. We refer the reader to F21 for further details on the photometric measurements, photometric redshifts, and sample validation. Note-The sample of photometric redshift selected z > 8.5 galaxies studied in this work are taken from F21. * These objects have spectroscopic redshifts, as listed in Table 3.

CONSTRAINING STELLAR POPULATION POSTERIORS WITH PROSPECTOR
We constrain the stellar populations by using Prospector (Johnson et al. 2021), a fully Bayesian inference code to derive stellar population properties from photometric and/or spectroscopic data. Prospector has been mainly employed on galaxies at lower redshifts (e.g., Leja et al. 2019b;Webb et al. 2020;Belli et al. 2021;Tacchella et al. 2021). The Prospector fit for one high-z galaxy (GOODSN-35589) has been presented in Johnson et al. (2021) as a demonstration. We adopt a similar physical model for the galaxy SED as in Johnson et al. (2021) with details given in Section 3.1, with Section 3.2 highlighting the prior on the SFH. Section 3.3 assesses the goodness of the SED fits and shows the dust attenuation curve posteriors. Finally, Section 3.4 compares the photometric redshifts from EAZY to the ones obtained here with Prospector.

Physical model for the galaxy SED
We use the Flexible Stellar Population Synthesis (FSPS) package (Conroy et al. 2009;Conroy & Gunn 2010) with the MIST stellar evolutionary tracks and isochrones (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). Fur-thermore, 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.
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 (Larson et al., submitted), GOODSN-35589 with z spec = 10.957 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 a 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 log(Z/Z ) 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 1 and are simply partitioned into time bins. The SFHs are characterized by the ratios of the SFRs in adjacent time bins. There are 5 free parameters for 6 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 < log(M /M ) < 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 birth-cloud dust optical depth clipped normal in (τ dust,1 /τ dust,2 ): min = 0, max = 2, µ = 1, σ = 0.3 log(Zgas/Z ) gas-phase metallicity uniform: min = −2.0, max = 0.5 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 (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τ dust,2 of the diffuse dust component (flat prior between 0 <τ dust,2 < 4.0); and (iii) the normalizationτ 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 <τ dust,1 /τ dust,2 < 2.0, motivated by Calzetti et al. 1994;Price et al. 2014). The nebular emission (emission lines and continuum) is self-consistently modeled . We have two parameters: the gas-phase metallicity (Z gas ) and the ionization parameter (U ). We assume a flat prior in logspace for the metallicity (−2.0 < log(Z gas /Z ) < 0.5) and ionization (−4 < log(U ) < −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 de-coupled 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.

Priors for the SFH
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;  Redshift Figure 1. Different choices for the SFH prior produce different behavior in SFR(t) (see also Leja et al. 2019a andTacchella et al. 2021). The panels from left to right show 100,000 random draws from the continuity prior (our fiducial SFH prior), the bursty continuity prior, the Dirichlet prior and the parametric prior (Section 3.2). The solid black lines mark the median in time for these draws, while the gray shaded regions indicate the 16-84th percentiles. The vertical yellow lines show the six time bins for the flexible SFHs. In each panel, seven individual draws (i.e., SFHs) are shown as blue lines to illustrate the behavior of different priors. The bursty continuity prior is weighted in order to produce multiple bursts and quenching episodes in comparison with the smoother continuity prior and Dirichlet prior. The parametric prior introduces a specific shape with an increasing SFH with time. All of these priors are able to fit the data equally well.
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. 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 < log(M /M ) < 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 oneto-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. 2021, 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 ∆ log(SFR) between adjacent time bins. We adopt the Student's t-distribution ∆ log(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: The parameter τ is varied as log(τ ) within a uniform prior in the range of −1.0 < log(τ ) < 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).

Resulting posteriors
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 Note-We list the galaxy identifier (ID), the redshift (photometric if not zspec specified), the absolute UV magnitude at restframe 1500Å (M UV,obs ), the UV spectral slope (β), the stellar mass (M ), the SFR averaged over past 50 Myr (SFR 50 ), the specific SFR averaged over past 50 Myr (sSFR 50 ), the dust attenuation at 5500Å (A V ), and the stellar metallicity (Z).
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. The detections below the Lyman-α break in COSMOS-20646 and COSMOS-47074 are discussed in F21 (see their Figure 15). Briefly, the detections are only marginally significant (SNR = 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 non-detections/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 χ 2 tot , which we estimate by summing the individual χ values of the detected bands. We do not quote the reduced χ 2 values as the number of degrees of freedom is not well defined in a non-linear 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.7 −0.1 , 0.9 +0.1 −0.3 , 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.

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. (2019a) 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 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 SMClike (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.
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).

Confirmation of the photometric redshifts
As discussed in Section 2 and F21, the photometricredshift 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), Leja et al. (2018) and Tacchella et al. (2021), which adds 5 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 −4 and 0.1), U min (minimum starlight intensity to which the dust  Prospector assume the free-z setup (see Section 3.4) with a uniform redshift prior between z = 0.1 − 13. Each panel shows the EAZY and Prospector z phot posteriors in black and red, respectively. Three galaxies (EGS-6811, EGS-68560 and GOODSN-35589) have spectroscopic redshifts, which are indicated as blue vertical lines. We find good agreement between EAZY and Prospector; see also Figure 5 for a direct comparison. 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 −5 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 errorbars in-dicate the median and 16-84th percentiles of the redshift posterior, respectively.
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 towards 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 fluxlimited survey will contain many more low-z than highz 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, since we have probably significantly overestimated the highz 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 χ 2 tot 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 χ 2 tot by a factor of 2-4 than the low-z run. Furthermore, we obtain for the Bayes factor a median of Z low−z /Z high−z = 3 × 10 −2 over the whole sample, with all galaxies Z low−z /Z high−z < 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 phot = 2.63 +0.16 −0.16 ) and EGS-20381 (z phot = 6.79 +0.57 −0.13 ). 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 phot = 7.40 +0.05 −0.04 ), where this alternative z phot estimate is inconsistent with the available spectroscopic redshift (z spec = 8.68).  Figure 6. Dependence of the UV spectral slope β on individual galaxy physical properties. The horizontal axis shows the fractional change (e.g., ∆) in a given parameter relative to the maximal value considered (with the actual values given by the upper horizontal axis). The red, green, and blue curves show the change in β for changing dust attenuation, stellar population age, and stellar metallicity, respectively (with the other two parameters held fixed, using fiducial fixed values of AV = 0.2 mag, age = 50 Myr, and Z = 0.2 Z , a constant SFH, and a Calzetti attenuation curve). While rising metallicity and age can affect β, changes in dust attenuation are much more significant, thus using β to study dust attenuation is warranted. This is especially true at these early epochs where the maximal stellar age is limited by the short period since the onset of star formation (∼ 200 − 300 Myr).

CHEMICAL ENRICHMENT IN EARLY BRIGHT GALAXIES
4.1. The UV Spectral Slope

β 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 Lyman-α 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 signif-icant 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. 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 log(M /M ) = 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.

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 powerlaw 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 a smaller scatter and a 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 towards 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.
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)  . The measured UV spectral slope β for our sample of 11 z = 9 − 10 galaxies (Section 4.1.2) versus their derived UV absolute magnitude (left) and stellar mass (right). We show previously published results for predominantly fainter galaxies as small symbols. In the left panel, the colored lines show the measured correlations between β and MUV at z = 4 − 8 from Bouwens et al. (2014), while with similar lines in the right panel we show the measured correlations between β and the stellar mass at z = 4 − 8 from Finkelstein et al. (2012), converting from Salpeter to a Chabrier IMF. Our sample of z = 9 − 10 galaxies appears to exhibit a strong correlation with the stellar mass (Pearson r = 0.85) and little-to-no correlation with MUV (Pearson r = −0.32).
β 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) 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 the 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 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 log(M /M ) 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.42 −0.49 . 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 β = −1.87 +0.35 −0.43 . While the median is not highly dependent on the inclusion of these two galaxies, we consider this latter value our fiducial value.
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. Simulations FLARES (Vijayan+21) SAM (Yung+19) THESAN (Kannan+21) Figure 8. The UV spectral slope β versus redshift. We show each of our individual galaxies as small purple circles. The large dark purple square shows the median measured β from our sample (calculated by stacking the posterior distributions on β), excluding the two galaxies with log(M /M ) > 10.7 as their mass may be biased high due to residuals from neighbors in the IRAC imaging (white filled). We compare to results at z = 4-8 from Finkelstein et al. (2012), shown by the lighter-colored squares, using color to denote the stellar mass. We find that our sample of observed modestly massive galaxies (log M/M = 9.1-10.6) have measured values of β comparable to similarly massive galaxies at z = 4-8. This implies that galaxies of these masses can grow their dust reservoirs in a relatively short period of time, as we are observing many of these galaxies <500 Myr from the Big Bang. This is consistent with predictions from multiple simulations (a semi-analytic model (Yung et al. 2019b), the FLARES simulations (Lovell et al. 2021;Vijayan et al. 2021a), and the THESAN radiation-magneto-hydrodynamic simulation (Kannan et al. 2021a;Garaldi et al. 2021;Smith et al. 2021), which predict significant dust reservoirs in these early massive galaxies.

Comparison of β to SED-fitting Results
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 50 2 . The stellar age t 50 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 in-form most strongly about dust attenuation. With the exception of the two bluest galaxies, all galaxies are constrained (at the 1σ level) to have non-zero levels of dust attenuation.
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 un-enriched galaxies) at lower masses at this same epoch. In . The dependence of the Prospector-derived physical quantities on the UV-slope β, showing V-band dust attenuation (left), log stellar metallicity (middle), and stellar population age t50 (right), the lookback time when 50% stellar mass has been formed. The color-coding denotes increasing log stellar mass as indicated with the colorbar on the right. The squares mark the galaxies with spectroscopic redshifts, while the circles show the galaxies with photometric redshifts. Consistent with our expectation from Figure 6, β appears to correlate most strongly with dust attenuation. We do not see a strong correlation with stellar metallicity, though the uncertainties on the Prospector-derived metallicities are high, so conclusions from this middle panel are not significant. We do see that the bluest galaxies appear to have the youngest derived stellar population ages, though at β > −2 there is no visible trend between β and age.
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 1500 = 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.

Implications on Evolution of UV LF
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 log(M /M ) = 9 − 10 are similar to those at the same UV luminosities and masses at z = 4 − 8 redshift. 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;Finkelstein et al. 2015b;Bowler et al. 2015), 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 starformation efficiency (e.g, Finkelstein et al. 2015a;Stefanon et al. 2019), or time/mass scales for the onset of quenching (e.g., Tacchella et al. 2018a;Bowler et al. 2020).

Implications on Dust Formation in Early Galaxies
The most surprising result from this analysis is that the UV spectral slope β for relatively massive UVselected systems (9 < log(M /M ) < 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;Schouws et al. 2021;Fudamoto 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. 2015Mancini et al. , 2016Popping et al. 2017;Behrens et al. 2018;Sommovigo et al. 2020;Graziani 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 non-trivial. 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 highredshift 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 lowski et al. 2010;Valiante et al. 2011;Asano et al. 2013;Mancini et al. 2015;Leśniewska & Micha lowski 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. 2015Marassi et al. , 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. 2018Jaacks et al. , 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 log(M /M ) > 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 semi-numerical (e.g., Mancini et al. 2015Mancini et al. , 2016 models. We compare our observations to predictions from a semi-analytic model (SAM; Yung et al. 2019b), the FLARES simulations (Lovell et al. 2021;Vijayan et al. 2021a), and the THESAN radiation-magnetohydrodynamic simulation (Kannan et al. 2021a;Garaldi et al. 2021;Smith et al. 2021) to our observations in Figure 8. The β values for the THESAN simulation are presented in Kannan et al. (2021b). While THESAN does track dust formation and destruction, it does not have any galaxies as massive as those we observe at z > 9. FLARES and the SAM 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. 2019b;Shen et al. 2020Shen et al. , 2021Vijayan et al. 2021b), 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.

GROWTH OF STELLAR MASS AND IMPLICATIONS ON EARLY COSMIC STAR FORMATION
This section presents the key results concerning the early mass growth histories inferred from our SEDmodeling analysis. In particular, we present the stellar mass and SFR measurements in Section 5.1, followed by an exploration of the inferred SFHs (Section 5.2) and the stellar ages and star-formation timescales (Section 5.3). Finally, we look into the fraction of mass formed beyond redshift 10 (Section 5.4) and the implications for the early evolution of the cosmic SFR density (Section 5.5). An important conclusion is that our SFHs depend on the assumed SFH prior, which we introduced in Section 3.2. We then discuss whether our galaxies are overly massive for the ΛCDM universe (Section 5.6) and how we can make progress in the future with JWST (Section 5.7).

Stellar masses and star-formation rates (SFRs)
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 50 . 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): where t SF = 50 Myr is the timescale of the SFR indicator. As an example, for a M = 10 10 M galaxy, the maximum SFR is SFR max = 200 M yr −1 . Similarly, the maximum specific SFR (sSFR) is given by:  Figure 10. Star-formation rate (SFR) and stellar mass (M ) properties of our z = 9 − 11 galaxies. The left and right panels show our galaxies in the SFR50 − M plane and the sSFR50 − z plane, respectively. The SFR50 and sSFR50 are averaged over the past 50 Myr. The red symbols show the fiducial continuity prior, while the smaller blue, purple and green symbols indicate the results for the bursty continuity prior, the Dirichlet prior, and the parametric prior. The squares mark the galaxies with spectroscopic redshifts, while the circles show the galaxies with photometric redshifts. The gray shaded regions indicate the forbidden parameter spaces, where SFR50 would be too high given 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 sSFR max = 20 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 50 − M and the sSFR 50 − z planes, respectively. The black lines and the gray shaded regions indicate the maximum SFR and sSFR mentioned above. The red datapoints and errorbars 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 datapoints 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.
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 9 −10 11 M at z > 4, where the "bending" of the star-forming main sequence plays presumably a minor role. Salmon et al. continuity prior bursty continuity prior Dirichlet prior parametric prior Figure 11. Posteriors of the star-formation rate (SFR50) and stellar mass (M ). SFR50 is averaged over the past 50 Myr. The red, blue, purple and green colors indicate the results from assuming the continuity prior, the bursty continuity prior, the Dirichlet prior, and the parametric prior for the SFH. These different priors give rise to similar stellar masses and SFRs (within a factor of 3). An exception is UDS-18697, where SFR50 varies by more than 1 order of magnitude, which stems from the degeneracy with the amount of attenuation. . 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 9 − 10 10 M . We measure sSFR values in the range of 3 − 10 Gyr −1 , 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.6 −0.3 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 −1 .
We have also added to the right panel of Figure  steeper increase with 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 50 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. 2021). The bursty continuity SFH prior shows the largest spread in the posterior space, in particular along the SFR axis.

Star-formation histories (SFHs)
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. Figure 12 emphasizes when most of the mass has formed, and it stresses that the SFH 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 para-metric 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 50 ) 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 50 on top and τ SF on the bottom. The colors and symbols are the same as those in Figure 10.

Inferred ages and star-formation timescales
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 three times younger ages (and shorter star- Redshift Figure 12. Star-formation histories (SFHs) obtained from Prospector assuming different priors. The SFHs are plotted as a fraction of the stellar mass formed. The adopted priors include the continuity prior (red), bursty continuity prior (blue), Dirichlet prior (purple) and parametric prior. These priors are discussed in Section 3.2. The lines and shaded regions show the median and 16-84th percentile of the SFH posterior, respectively. In most cases (exception is the most massive galaxy UDS-18697), the bursty continuity prior leads to more recent star formation than the continuity and Dirichlet priors, which both roughly follow a constant SFH. The parametric prior lies in between those extrema. The data do not prefer any of those priors: the Bayes factor for all models with respect to each other is roughly 1, highlighting that the adopted prior heavily effects the resulting posterior SFH.
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 50 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 50 ∼ 10 − 12. These ages are consistent with expectations from the semi-analytical model by Yung et al. (2019b) and the empirical halo model by Tacchella et al. (2018a), which find stellar ages t 50 for galaxies at z ≈ 10 of 60 Myr and 40 Myr, respectively.
In Appendix C, Figures 21 and 22 show the posterior distribution of M versus t 50 and UV attenuation A UV versus t 50 , respectively. We show that there is a degeneracy between M and t 50 : 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 anti-correlation 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 . Stellar age t50 (top panels) and star-formation timescale τSF (bottom panels) as a function of stellar mass M (left panels) and redshift z (right panels) of our z = 9 − 11 galaxies. The star-formation timescale τSF is defined as the time it took to increase the stellar mass from 20% to 80%, i.e., it is a measure of how quickly a galaxy formed its stellar mass. The pink shaded region in the bottom right panel indicates 30 − 40% of the age of the universe at a given redshift z. The red symbols show the fiducial continuity prior, while the smaller blue, purple, and green symbols indicate the results for the bursty continuity prior, the Dirichlet prior, and the parametric prior. The squares mark the galaxies with spectroscopic redshifts, while the circles show the galaxies with photometric redshifts. The exact values of t50 and τSF depend heavily on the assumed SFH prior: the bursty continuity prior implies younger ages and shorter star-formation timescales, while for the fiducial continuity and Dirichlet prior τSF older ages and longer τSF. Independent of the SFH prior, we find older ages for more massive galaxies and longer star-formation timescales for the lower-redshift galaxies.
used the FAST code (Kriek et al. 2009) and assumed a fixed sub-solar 0.2 Z metallicity, a constant SFH, and a fixed Calzetti et al. (2000) attenuation law. Our 10 9 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 V = 0.15 +0.30 −0.15 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 3 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 50 (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 20 (lookback age at which 20% of the Table 4. Stellar age (t50), star-formation timescale (τSF), and fraction of mass formed before z = 10 (fM) for our galaxy candidates at z = 9 − 11. We list the results for the four different SFH priors (Section 3.2). 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 +59 −105 Myr and 265 +58 −90 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 20 when adopting the bursty continuity (91 +183 −72 Myr) or parametric (163 +117 −99 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.

Fraction of mass formed at z > 10
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 50 . 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). 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 > 10) = 53 +45 −51 %, 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 +18 −28 %, 39 +61 −39 %, 82 +14 −21 %, and 65 +31 −45 % 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.

Implications for the cosmic SFR density
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) −11 (e.g., Oesch et al. 2014, continuity prior bursty continuity prior Dirichlet prior parametric prior Figure 14. Fraction of mass formed previous to redshift z = 10 (fM(z > 10)). For each galaxy, we plot the posteriors of fM(z > 10) and stellar age (t50). As expected, we find a strong correlation between fM(z > 10) and age: older galaxies have a higher fM(z > 10). More importantly, fM(z > 10) (as the age) depends heavily on the assumed prior (shown with the different colors). Since GOODSN-35589 lies at zspec = 10.96, fM(z > 10) is equal to 1.
However, there are also other models that prefer a rather slow increase of 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, dash-dotted, and dashed black lines indicate f M (z) for a constant SFR, the rapidly increase cosmic SFR density (∝ (1 + z) −11 ), and a slowly increase cosmic SFR density (∝ (1 + z) −4 ), respectively.
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 yet 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 (Fig. 12), which is of interest for the interpretation of the 21cm  Figure 15. Cosmic SFR density implied by our inferred SFHs. We compute the fraction of stellar mass formed from our posterior SFHs averaged over all 10 galaxies that lie at z < 10 (excluding GOODSN-35589). The red, blue, purple and green lines show the continuity prior, the bursty continuity prior, the Dirichlet prior and the parametric prior for the SFH, respectively. The dotted black line shows the evolution of fM for a constant SFR, while the dashed and dash-dotted lines indicate the rate of decline in the UV luminosity density deduced from large photometric surveys (Oesch et al. , 2018Bouwens et al. 2015;McLeod et al. 2016). The SFH prior plays a crucial role: assuming the continuity prior or the Dirichlet prior leads to a rather constant SFH, consistent with the rather slow increase with time (∝ (1 + z) −4 ; dash-dotted line) as inferred in some UV luminosity studies. On the other hand, assuming the bursty continuity prior or the parametric prior leads to a more rapid and recent increase in fM, consistent with the steep increase with time (∝ (1 + z) −11 ; dashed line) as inferred in other UV luminosity studies.
signal and the formation and evolution of black holes in the early universe (e.g., Dayal & Ferrara 2018).

Overly massive galaxies
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 6 Mpc 3 , 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 less likely in Planck CDM > 10 6 Mpc 3 continuity prior bursty continuity prior Figure 16. Are the z = 9 − 11 galaxy candidates overly massive? We plot the stellar mass M as a function of redshift z for the galaxies in our sample. The red and blue symbols mark the measurement adopting the continuity and the bursty continuity SFH prior, respectively. The red and blue lines show the inferred mass growth histories from the SFH measurements. The solid black line indicates the threshold stellar mass for a cumulative number density of Φ = 10 −6 Mpc −3 (Behroozi & Silk 2018, converted to Chabrier IMF), which is roughly our survey volume. The two galaxies with stellar masses larger than the black line (i.e., in the gray shaded region) are in tension with ΛCDM, though the significance of this tension is low as we discuss in the text.
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 −6 Mpc −3 , 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 since 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). 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 (log(M /M ) = 10.6 +0.2 −0.3 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 11 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 log M , 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.

Towards JWST
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 mediumband 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 errorbars at the bottom right illustrate the 5σ uncertainty for ∼ 15 ksec exposures, estimated from the point source limit with a 0.1 aperture and medium background (Williams et al. 2021).
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β+ [OIII] ) 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 EW Hβ+[OIII] = 675 +2213 −578Å ,   −2388Å , respectively. These different EW measurements are then directly reflected in the color-color diagram since 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 NIR-Spec/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).

SUMMARY & CONCLUSIONS
In this work we carefully assess the current constraints on the stellar populations of galaxies at z = 9 − 11. The sample consists of 11 bright (H < 26.6) galaxies, out of which 3 have spectroscopic redshifts (F21). Given the high-redshift nature of these sources, HST and Spitzer /IRAC data only trace the rest-frame UV and the Balmer/4000Å-break of these galaxies.
We perform a careful inference of the stellar populations by using Prospector, a flexible Bayesian SED fitting code (Johnson et al. 2021). In particular, we ex-pand upon previous z > 6 SED investigations by adopting a range of flexible and parametric SFHs, a flexible dust attenuation law, self-consistent modeling for emission lines, and variable IGM absorption. A flexible attenuation law is important because it allows us to assess how degenerate metallicity, age, and SFH constraints behave with the attenuation when extracting this information from a rather limited wavelength coverage at low spectral resolution (e.g., Kriek & Conroy 2013;Battisti et al. 2016;Salim et al. 2018;Tacchella et al. 2021). The self-consistent modeling of emission lines is essential because it allows us to estimate their contribution to the IRAC fluxes and to address the degeneracy with a possible Balmer/4000Å break also present in those bands (e.g., Stark et al. 2013;Hashimoto et al. 2018). Finally, the different SFH priors are crucial as we have little knowledge about the "burstiness" (i.e., starformation variability) of these systems (e.g., Smit et al. 2016;Faucher-Giguère 2018;Faisst et al. 2019;Iyer et al. 2020;Tacchella et al. 2020).
Our SED-modeling approach (Section 3) takes into account all of the aforementioned points. The strength of Prospector is its fully Bayesian nature and its flexibility, which allows us to investigate how different assumptions (i.e., priors) affect our results and conclusions. In particular, the choice of the prior for the SFH (Section 3.2) has important consequences on conclusions regarding the early growth of these z = 9 − 11 galaxies. We investigate the impact of four priors: the continuity prior, the bursty continuity prior, the Dirichlet prior, and the parametric prior. The behaviors of these different priors are shown in Figure 1. The continuity prior and Dirichlet prior are both weighted toward a smooth behavior with a constant SFR as the expectation value. The bursty continuity prior allows for bursty star formation, where most of the mass is formed in one or a few time bins. Finally, the parametric prior assumes a delayed-τ model and -at these early times -typically leads to an increasing SFH. We find that all these SFH priors can equally well describe the observational data, and the data do not prefer any of the priors (Bayes factors are about 1).
We measure the rest-frame UV spectral slope (β) from the Prospector posterior model spectra for our sample of galaxies and find a significant correlation between β and the stellar mass, where more massive galaxies are redder. We do not measure a significant correlation between β and the UV luminosity in our sample. The measured UV slopes of our massive (log(M /M ) = 9 − 10) z ∼ 9-10 galaxies are similar to measured values at the same stellar masses at z = 4-8, indicating that galaxies at these stellar masses rapidly develop a dust reservoir, which then grows more slowly. This is consistent with galaxy formation models that include a significant contribution of dust grain growth. The roughly constant attenuation at these masses across a wide range of redshift implies that the apparent lack of evolution in the number density of UV-bright galaxies at z > 7 is not due to changes in the dust attenuation.
Despite the flexibility in the SED modeling and the investigation of different SFH priors, we find that the stellar masses and the SFRs (averaged over 50 Myr) are rather well constrained with uncertainties of a factor of 2. We find a hint of a star-forming main sequence with a sub-linear slope (0.7 ± 0.2; Figure 10), i.e., more massive galaxies have a higher SFR. The sSFRs are in the range of 3 − 10 Gyr −1 , which indicates a mass-doubling timescale of ∼ 100 − 300 Myr under the assumption of a constant sSFR.
Stellar population parameters such as the stellar metallicity, the stellar age, and the SFH itself are less well constrained. For example, changing the SFH prior leads to changes in the age of a factor of ∼ 3 (Figure 13): the median measured age over our sample is 147 +78 −89 Myr and 53 +153 −40 Myr for the continuity prior and the busty continuity prior, respectively. Even when just selecting one prior, the uncertainty in age for an individual galaxy remains rather large, which can be largely attributed due to the flexibility in the dust attenuation law and the emission line modeling (Figure 22). More generally, we find that the SFH priors have an important impact on the measured SFHs (Figures 12 and 20).
From this, we draw the important conclusion that the current observational data cannot give us tight constraints on how quickly these z = 9 − 11 galaxies are building up their mass (Figure 15). In the case of the bursty continuity prior and the parametric prior, which both prefer younger ages, the inferred stellar mass buildup is consistent with a rapidly increasing cosmic SFR density at z > 8 with time (∝ (1 + z) −11 ). In contrast, the continuity prior and the Dirichlet prior, both preferring older ages, are consistent with a slow increase with time in the cosmic SFR density in this epoch (∝ (1 + z) −4 ). Therefore, the epoch of first galaxy formation remains to yet be identified. JWST with its high sensitivity, larger wavelength coverage, and higher spectral resolution (including medium band imaging and spectroscopy) will help solve this mystery. This in turn will help to constrain galaxy formation models and might even shed light onto the nature of dark matter (Dayal et al. 2015;Khimey et al. 2021).

ACKNOWLEDGMENTS
We thank the referee for a thorough report that improved and strengthened this work. We are grateful to Ben Johnson and Joel Leja for helpful discussions related to Prospector. We thank Nicolas Laporte for clarifications and useful comments. S.T. is supported by the Smithsonian Astrophysical Observatory through the CfA Fellowship and by the 2021 Research Fund 1.210134.01 of UNIST (Ulsan National Institute of Science & Technology). S.L.F. acknowledges support from NASA through ADAP award 80NSSC18K0954. L.G. and R.S. acknowledge support from the Amaldi Research Center funded by the MIUR program "Dipartimento di Eccellenza" (CUP:B81I18001170001). I.J. acknowledges support from NASA under award number 80GSFC21M0002. Support for Hubble Space Telescope program #15862 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Associations of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.   Figure 18. Impact of IRAC photometry on resulting photometric redshift (z phot ) posteriors. Following Figure 4, each panel shows the z phot posteriors obtained by EAZY, by Prospector with TPHOT (fiducial) IRAC photometry, and by Prospector with GALFIT IRAC photometry in black, red, and orange, respectively. Three galaxies (EGS-6811, EGS-68560 and GOODSN-35589) have spectroscopic redshifts, which are indicated in blue. We find a significantly different z phot in 3 out of the 11 galaxies (COSMOS-206464, EGS-6811, and EGS-20381), while the other galaxies have a consistent z phot within the uncertainty. Importantly, the z phot of EGS-6811 with GALFIT-based photometry is inconsistent with its zspec.

A. GALFIT-BASED IRAC PHOTOMETRY
We discuss extensively the uncertainty of the IRAC photometry in F21. In that work, we have performed the deblending and flux measurements of the IRAC photometry with both TPHOT (our fiducial approach) and GALFIT (see Table 6 in F21). We discuss in this section how the systematic uncertainty of the IRAC photometry impacts our results. Figure 18 shows the photometric redshift posterior when adopting TPHOT (solid red lines) or GALFIT (dashed orange lines) IRAC photometry. Similar to Figure 4, we also plot the redshift posterior of EAZY for comparison, and the vertical blue lines indicate the spectroscopic redshift when available. This figure assumes the "z-free model" introduced in Section 3.4.
We find that COSMOS-20646 with GALFIT-based IRAC photometry strongly prefers a z ∼ 2.5 solution, consistent with the finding by F21. Furthermore, EGS-6811 and EGS-20381 shift mildly but significantly toward lower redshifts (z ∼ 7). However, in the case of EGS-6811, the GALFIT-based photometric redshift is inconsistent with the spectroscopic redshift. For all other galaxies, the GALFIT-based photometric redshifts are consistent with the TPHOT-based ones within the uncertainties. Figure 19 investigates the SFRs and stellar masses obtained when swapping between TPHOT-and GALFIT-based IRAC photometry. Here we assume our standard model, i.e., the continuity SFH prior and the EAZY redshift posterior as the redshift prior. For most galaxies, the SFR and M posteriors are consistent with each other, but exceptions are EGS-6811 and UDS-18697. EGS-68560 and UDS-18697 are more massive when using the GALFIT IRAC photometry, continuity prior (TPHOT) continuity prior (Galfit) Figure 19. Posteriors of the star-formation rate (SFR50) and stellar mass (M ) as in Figure 11, but for the TPHOT (fiducial) and GALFIT IRAC photometry in red and orange, respectively. We find that changes in the IRAC photometry lead to differences in the stellar mass and SFR that can be larger than the uncertainties.
while COSMOS-20646, EGS-44164 and EGS-26890 have a tendency to be less massive. The SFRs are overall similar, but we find a significant increase for EGS-6811 and UDS-18697. This also has important implications for the discussion whether our galaxies are too massive at these early cosmic times (Section 5.6). Both COSMOS-20646 and UDS-18697 have stellar masses of M ≈ 10 11 M and therefore lie solidly within the less likely ΛCDM region in Figure 16. Figure 19 shows that the stellar mass of COSMOS-20646 reduces by roughly a factor of 3 when adopting the GALFIT photometry, making it consistent with the ΛCDM boundary when considering the uncertainty. On the other hand, UDS-18697 seems to get even more massive (and star forming) when adopting the GALFIT photometry.
In summary, these differences depending on which photometric method is used highlight a systematic uncertainty when making use of deblended photometry from low-resolution imaging, something which will be alleviated soon with JWST. We adopt throughout this work the TPHOT-based IRAC photometry as it does overall a better job of removing the bright neighboring sources (F21). The object COSMOS-20646 has the largest uncertainty regarding which method is more accurate; we therefore caution the reader that the nature of this galaxy is still somewhat uncertain.
B. STAR-FORMATION HISTORIES: SFR VERSUS TIME For completeness, Figure 20 plots the inferred SFHs as the SFR as a function of time, while Figure 12 in the main text plots the SFHs as fraction of mass formed as a function of 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.
For the fiducial continuity and Dirichlet priors, the resulting SFHs are similar and roughly constant with time. For a few galaxies (i.e., COSMOS-20646, EGS-68560, GOODSN-35589, and UDS-18697), there is significant variation in the past ∼ 100 Myr. The otherwise rather constant behavior is expected, as the expectation value of this prior is 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 Redshift Figure 20. Star-formation histories (SFHs) obtained from Prospector assuming different priors. The adopted priors include the continuity prior (red), bursty continuity prior (blue), Dirichlet prior (purple), and parametric prior. These priors are discussed in Section 3.2. The lines and shaded regions show the median and 16-84th percentile of the SFH posterior, respectively. The continuity and Dirichlet prior typically lead to a rather constant SFH, the parametric prior (delayed-tau parameterization) to an increasing SFH, and the bursty continuity prior to a bursty SFH, which biases the median SFH low (see text). The key point of this figure is that the adopted prior heavily effects the resulting posterior SFH.
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 only provide little constraining power when it comes to the SFH. The median SFH of the bursty continuity prior seems to lie significantly below the other SFHs, which -at first glimpse -implies a lower stellar mass. However, this is not the case as we have showed in Figure 11. The explanation is that this prior tends to form the stars in one or a few time bins, which leads to a peaked SFH. When computing the median (as a function of time), this leads to a bias low, but peaky behavior can still be seen in the 16-84th percentiles (shaded region). In order to circumvent this problem, we can normalize first by mass before computing the median. This is done in Figure 12.  Figure 11 in the main text. Figure 21 shows that there is a degeneracy between M and t 50 : a younger age implies a lower stellar mass. This is expected as younger stellar populations are typically brighter at fixed stellar mass. There is also a clear rank ordering of the SFH prior: the bursty continuity prior produces younger ages, followed by the parametric prior, while the Dirichlet and the continuity priors produce the oldest galaxies. As these galaxies lie within a redshift range of z ≈ 9 − 11, we also give the formation redshift z 50 , i.e., the redshift by which 50% of the mass of the galaxy has formed. We find that these galaxies typically form around z f ≈ 11 − 13, with UDS-18697 forming the earliest at z f ∼ 15. Importantly, z f should not be confused with the epoch COSMOS-20646 z50 = 12.4 +1.5 1.6 z50 = 11.8 +2.8 1.0 z50 = 13.1 +1.5 1.4 z50 = 11.7 +1.1 0.9 zphot = 9.83 COSMOS-47074 z50 = 12.1 +1.6 1.5 z50 = 10.0 +1.0 0.2 z50 = 12.5 +1.7 1.3 z50 = 11.3 +1.3 0.9 zspec = 8.68 EGS-6811 z50 = 11.6 +1.4 1.4 z50 = 10.1 +2.8 1.1 z50 = 11.7 +1.6 1.4 z50 = 9.7 +1.2 0.7 zspec = 8.66 EGS-44164 z50 = 11.5 +1.4 1.5 z50 = 9.4 +1.4 0.5 z50 = 11.5 +1.7 1. EGS-20381 z50 = 11.2 +1.7 1.6 z50 = 9.1 +1.8 0.5 z50 = 11.1 +1.8 1.0 z50 = 9.7 +1.0 0.8 zphot = 9.06 EGS-26890 z50 = 11.3 +1.6 1.3 z50 = 9.4 +0.9 0.3 z50 = 11.6 +1.8 1.1 z50 = 10.5 +1. continuity prior bursty continuity prior Dirichlet prior parametric prior Figure 21. Posteriors of stellar mass (M ) and stellar age (t50). The red, blue, purple, and green colors indicate the continuity prior, the bursty continuity prior, the Dirichlet prior, and the parametric prior for the SFH. We also indicate the formation redshift z50 at which 50% of the stellar mass has been formed. Although both the stellar mass and the age are consistent within the uncertainty with each other, assuming different priors, the absolute values for the ages vary significantly. Consistent with Figure 12, the bursty continuity prior and the parametric prior lead to younger ages than the continuity prior and the Dirichlet prior.
when the star-formation initially started. As we can see in Figure 20, star formation in most galaxies start around z ≈ 20, which is basically by construction as constraining this is out of reach for the current data (Section 3.2). Figure 22 shows that it is challenging to break the dust-age degeneracy with our current observational data. Specifically, we find an anti-correlation 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. This can at least in part be explained by the degeneracy with the attenuation law, over which we marginalize here. Furthermore, we do not find any trend with the SFH prior.

D. FRACTION OF MASS FORMED AT Z > 12
We present in Section 5.4 the results on the fraction of mass formed before z = 10 (f M (z > 10); see also Table 4 and Figure 14). We find that all of our z = 9 − 11 galaxies have formed at least 50% of their mass before z = 10, though the exact number depends on the assumed prior. In particular, for the bursty continuity prior, for some galaxies f M (z > 10) drops to less than 10%.