Outshining by Recent Star Formation Prevents the Accurate Measurement of High-z Galaxy Stellar Masses

We demonstrate that the inference of galaxy stellar masses via spectral energy distribution (SED) fitting techniques for galaxies formed in the first billion years after the Big Bang carries fundamental uncertainties owing to the loss of star formation history (SFH) information from the very first episodes of star formation in the integrated spectra of galaxies. While this early star formation can contribute substantially to the total stellar mass of high-redshift systems, ongoing star formation at the time of detection outshines the residual light from earlier bursts, hampering the determination of accurate stellar masses. As a result, order-of-magnitude uncertainties in stellar masses can be expected. We demonstrate this potential problem via direct numerical simulation of galaxy formation in a cosmological context. In detail, we carry out two cosmological simulations with significantly different stellar feedback models, which span a significant range in SFH burstiness. We compute the mock SEDs for these model galaxies at z = 7 via calculations of 3D dust radiative transfer, and then backward fit these SEDs with prospector SED fitting software. The uncertainties in derived stellar masses that we find for z > 7 galaxies motivate the development of new techniques and/or priors for SFH to model star formation in the early Universe.


INTRODUCTION
The most common method for determining the stellar mass of a galaxy is through ultraviolet-near infrared spectral energy distribution (SED) modeling.This technique, first developed for SED fitting have been explored in recent years (Brammer et al. 2008;Kriek et al. 2009;da Cunha et al. 2008;Chevallard & Charlot 2016;Iyer & Gawiser 2017;Carnall et al. 2018;Johnson et al. 2021).
The assumed form for the model star formation history in SED fitting software is an essential element in deriving galaxy stellar masses.Traditional functional forms for the SFH include constant, exponential declining, burst models, and combinations of these amongst others (Conroy 2013).These parameterized forms for SFH (hereafter, "parametric" SFHs) have parameters describing the (for example) normalization, -folding time and amplitude of bursts that are varied until a match is found between the synthetic SED produced by the SPS model and the observed data.When a solution is found, the model SFH is then used to infer the stellar mass of the observed galaxy (typically assuming a fixed metallicity).Of course, the assumed form of this SFH can severely impact the modeled stellar mass (Michalowski et al. 2012;Simha et al. 2014;Acquaviva et al. 2015;Iyer & Gawiser 2017;Carnall et al. 2018Carnall et al. , 2019;;Dudzevičiūtė et al. 2020).
More recently, a number of codes have explored the impact of more flexible so-called "non-parametric" forms for the model SFH (e.g.Heavens et al. 2000;Tojeiro et al. 2007;Iyer & Gawiser 2017;Johnson et al. 2021).Non-parametric SFH models do not have an explicit functional form as the parametric models, but instead can vary the amplitude of the SFH over a number of redshift or time bins in the modeled history of the galaxy.Iyer et al. (2018) validated the usage of non-parametric SFHs constructed via Gaussian processes by ground-truthing these methods against the the Santa Cruz Semi-Analytic Model (SAM) and the mufasa cosmological simulation (Davé et al. 2016(Davé et al. , 2017)).Similarly, Lower et al. (2020) demonstrated the efficacy of non parametric SFH techniques by ground-truthing modeled mock SEDs from galaxy simulations against their true stellar masses.Lower et al. (2020) found that while traditional parametric forms for the SFH in SED fitting software had uncertainties at the level ∼ 0.4 dex, non-parametric SFH models reduced these uncertainties to a level of ∼ 0.1 dex.Leja et al. (2019b) found that observed galaxies from the 3D-HST catalog (Brammer et al. 2012;Skelton et al. 2014) are systematically more massive and older when using non-parametric SFHs as compared to parametric methods, which bring the observed main sequence in line with theoretical predictions, potentially alleviating the long-standing tension between theory and observations in the SFR of galaxies at  ≈ 2 at a fixed stellar mass (Davé 2008;van Dokkum 2008).
With the successful launch of the JWST in 2021, observations are characterizing the physical properties of galaxies at unprecedented redshifts (e.g.Adams et al. 2023;Atek et al. 2023;Castellano et al. 2022;Donnan et al. 2023;Finkelstein et al. 2022b,a;Harikane et al. 2023;Labbé et al. 2023;Naidu et al. 2022;Robertson et al. 2023).Beyond providing constraints on the stellar mass buildup of some of the earliest galaxies in the Universe, JWST observations of high-redshift galaxies have been used to demonstrate potential tensions with the standard ΛCDM model (e.g.Boylan-Kolchin 2023;Haslbauer et al. 2022;Labbé et al. 2023;Lovell et al. 2019;Mason et al. 2023).What has yet to be explored, however, is the ability of traditional SED fitting methods to accurately derive stellar masses in these earliest galaxies, and in particular, the problem of stellar 'outshining'.
The problem of stellar outshining refers to the light from recent bursts of star formation overwhelming that of older stellar populations, making the inference of their stellar masses difficult.The potential impact of outshining has a rich history in the observational literature.Papovich et al. (2001) demonstrated an uncertainty of ∼ 3 − 8 in the derived stellar masses of  ∼ 2 − 3 Lyman Break Galaxies due to uncertainties in the star formation history modeling.This result has been further studied by Shapley et al. (2001); Daddi et al. (2004); Shapley et al. (2005);Trager et al. (2008); Graves & Faber (2010); Sorba & Sawicki (2018); Tacchella et al. (2022); Topping et al. (2022); Whitler et al. (2023) and Giménez-Arteaga et al. (2023), amongst many other studies.From a theoretical standpoint, Maraston et al. (2010), Pforr et al. (2013), andSuess et al. (2022) constructed mock star formation histories combined with stellar population modeling to study the potential impact of outshining.
The potential impact of stellar outshining by recent star formation in stellar mass inference has yet to be studied in bona fide cosmological hydrodynamic galaxy formation simulations.In this paper, we use numerical simulations of galaxy formation, combined with post-processed radiative transfer and SED fitting software, to study the impact of outshining on stellar mass estimates in galaxies, with a particular emphasis on high- galaxy detections by JWST.We demonstrate that standard SED fitting techniques have a difficult time deriving correct stellar masses in high- galaxies owing to outshining by recent bursts of star formation.We find that early bursts of star formation can contribute significantly to the total stellar mass of a galaxy.By the time the galaxy is detected at relatively late times ( ≈ 7 − 10), the current ongoing star formation outshines the light from evolved stars, making the integrated buildup of stellar mass difficult to measure.We demonstrate this effect using similar techniques to Lower et al. (2020)'s study of the efficacy of non-parametric SFH modeling: we simulate the formation of galaxies in cosmological simulations, forward model their mock SEDs using dust radiative transfer, and then fit these mock SEDs using standard techniques.This allows us to compare to ground truth from the simulations, and assess the efficacy of the fitting procedure.
In the remainder of this paper, we expand on these points.In § 2, we describe our numerical methods; in § 3, we describe the physical and luminous properties of the galaxies that we model; in § 4, we demonstrate the main issues when SED fitting very young galaxies; in § 5, we provide discussion, and in § 6, we summarize our main results.

Summary of Methods
Our main goal is to create mock SEDs from galaxies formed in a cosmological simulation, and then fit those SEDs as an observer would in order to ground-truth SED fitting techniques for  > 7 galaxies.To do this, we first simulate the highredshift evolution of galaxies by conducting cosmological hydrodynamic galaxy evolution simulations.We then "forward model" the emission from these galaxies by coupling them with dust radiative transfer in order to generate their mock SEDs.With these SEDs in hand, we then "backward model" them in order to derive the inferred galaxy physical properties.Through this methodology, we determine the relationship between the inferred stellar masses of our modeled galaxies from their SEDs, and the true stellar mass.

Galaxy Formation Simulations
We employ two rather different cosmological galaxy formation models in order to ensure the robustness of our results.The primary differences in these models as far as this study is concerned are: 1.The stellar feedback model, which impacts the burstiness of the modeled star formation histories.
2. The dust model, which impacts the dust content, extinction, and attenuation of the emergent light.
In detail, we employ the simba cosmological simulation, as already run by Davé et al. (2019), as well as a newly run simulation employing the smuggle explicit stellar feedback model within the arepo code base.We briefly describe these models in detail in turn below.The simba simulation is based on the gizmo cosmological gravity plus hydrodynamic solver (Hopkins 2015), and evolves dark matter and gas elements together, including gravity and pressure forces.Gas cools radiatively using the grackle library (Smith et al. 2017), including both metal line cooling, as well as non-equilibrium evolution of primordial elements.Stars form in molecular H 2 gas following the Krumholz et al. (2009) subresolution prescription for determining the HI and H 2 content in a gas particle, as well as the Kennicutt (1998) prescription for the star formation rate.Here, a star formation efficiency is manually set to  * = 0.02.The gas itself is artificially pressurized in order to resolve the Jeans mass as described in Davé et al. (2016).The upshot of this is that the star formation history for a given galaxy tends to occur more smoothly than in models with an explicit feedback model (and no artificial pressurization of the ISM; Iyer et al. 2020).Dust is modeled within simba following the algorithms outlined in Li et al. (2019).Specifically, dust is included as a single sized particle, though can comprise of multiple species (graphites and silicates).Dust is formed in evolved stars (Dwek 1998), can grow via metal accretion (Asano et al. 2013), and can be destroyed via thermal sputtering in hot gas (Tsai & Mathews 1995;McKinnon et al. 2017;Popping et al. 2017), supernova blastwaves, and astration in star-forming regions.
In addition to the simba simulation, we have run simulation with the smuggle galaxy formation model enabled in the arepo code (Springel 2010;Weinberger et al. 2020).We refer the reader to Marinacci et al. (2019) for a full description of this model, and highlight only the key differences from the simba simulation as they pertain to our study.Star formation occurs only in gravitationally bound molecular gas (Hopkins et al. 2013), and follows a volumetric Kennicutt (1998) law, though with an efficiency  * = 1.This is in contrast to the forced inefficiency of star formation in simba because on long timescales, stellar feedback in the smuggle model selfregulates the star formation rates to result in the relatively low efficiencies observed in molecular clouds (Marinacci et al. 2019).Stellar feedback models include supernovae, radiative feedback, stellar winds, and thermal feedback from HII regions.Like simba, dust is also included in the smuggle model.In contrast, however, dust is modeled to include a spectrum of grain sizes that evolve as the dust grains evolve in the simulation.Beyond the aforementioned dust processes that are included in simba, the dust in smuggle is allowed to grow in size via coagulation (sticking together), as well as fragment into smaller grains via grain-grain shattering collisional processes.The upshot from this dust modeling is that the local extinction law is explicitly computed in a spatially resolved sense in galaxies (Li et al. 2021), and represents therefore a fundamental difference in the forward modeling of radiative transfer between the smuggle and simba models.We have simulated a 25/ℎ Mpc side-length box with periodic boundary conditions, starting from  = 99 with initial conditions generated with music (Hahn & Abel 2011).We have run the model at the same mass resolution as the simba simulation (2 × 512 3 particles), though the initial conditions are generated with different random seeds, so the galaxies are not directly mappable from one simulation to another.We allow the simulation to evolve to redshift  = 6, and restrict our analysis to this redshift range.

Dust Radiative Transfer
In order to generate the mock SEDs (that we will then fit using prospector), we employ the public powderday dust radiative transfer package (Narayanan et al. 2021), which employs yt, fsps and hyperion for grid generation, stellar population synthesis, and Monte Carlo radiative transfer respectively (Turk et al. 2011;Conroy et al. 2009;Robitaille 2011).
Here, stars that form in the galaxy formation simulation emit a stellar spectrum based on their ages and metallicities.This spectrum is computed using fsps (Conroy et al. 2009(Conroy et al. , 2010;;Conroy & Gunn 2010).We assume the mist stellar isochrones (Choi et al. 2016), and a Kroupa (2002) stellar initial mass function (consistent with both sets of hydrodynamic simulations).We note that the choices of these parameters can subtly impact our results.For example, different isochrone models can change the assumed lifetimes of massive stars, which will impact the degree of outshining that we model here.
The light from these stars1 is emitted in an isotropic manner, and can be absorbed, scattered, and re-emitted by dust in the individual cells in the galaxy.For the simba simulations, which are particle based, the dust information is smoothed from the particles onto an adaptive mesh with an octree memory structure, and the radiative transfer occurs on this grid.For simba we assume Weingartner & Draine (2001) dust extinction laws locally.For smuggle, we perform the radiative transfer on a Voronoi mesh built around the dust particles simulated with the active dust model, and compute the extinction laws explicitly in each cell following Narayanan et al. (2023).Here, we assume extinction efficiencies from Draine & Lee (1984) and Laor & Draine (1993) for silicates and carbonaceous grains respectively.
Beyond attenuation by the diffuse dust in the galaxy (which is explicitly modeled in both simba and smuggle), powderday includes the possibility of obscuration by subresolution birth clouds following the Charlot & Fall (2000) formalism as built into fsps.Here, the attenuation only occurs for star particles below a threshold age (we set this to  = 10 Myr when included), and has a user-defined normalization to the attenuation curve, that we set as  V = 1.We discuss this in further detail in § 3.2.

THE PHYSICAL AND LUMINOUS PROPERTIES OF 𝑍 > 7 GALAXIES 3.1. Stellar Masses and Star Formation Histories
We begin our analysis by describing the physical and luminous properties of our model galaxies.In Figure 1, we plot a distribution of the stellar masses at  = 7 of the 20 most massive galaxies in the simba and arepo simulations.In a (25/ℎ) 3 volume, the most massive galaxies at  ≈ 7.5 are  * ∼ 10 8 − 10 9  ⊙ .These stellar masses are, by and large, built up via a series of individual star formation episodes, though the importance of these bursts to the total stellar mass buildup (and, in general, the shape of the SFH) is dependent on the assumed stellar feedback model (e.g.Sparre et al. 2017;Iyer et al. 2020).In Figure 2, we show the SFHs for 5 arbitrarily chosen galaxies from each simulation.The SFH is constructed from the  ≈ 7 snapshot, using the  = 7 stellar ages, correcting for mass loss processes.
The simba model with a manifestly pressurized ISM model broadly has a smoother rising star formation history, with punctuated bursts superposed.In comparison, the smuggle stellar feedback model is dominated by individual bursts at these high redshifts, owing to the explicit nature of the feedback model.The two models bracket the smoother SFHs typically seen in simulations with a pressurized ISM (e.g.eagle, illustris, and simba), and the burstier SFHs seen in explicit feedback models (e.g.fire and smuggle).It is the ability to reconstruct these SFHs with reasonably high fidelity (or lack thereof) that will prove essential for the SED fitting software to accurately derive the stellar masses of galaxies in this epoch.

Model SEDs
In the top row of Figure 3, we show the model SEDs for the same arbitrarily chosen galaxies in Figure 2 at  = 7.5.As The former has an artificially pressurized ISM (as do most cosmological galaxy formation models), while the right represents a class of explicit feedback models.The two models highlight the diversity of modeled star formation history burstiness possible in galaxy formation models.We use both of these models in our analysis to model a range of possibilities when investigating the efficacy of SED fitting in young galaxies.Note that the sudden drop off in the SFH at  < 7 is an artifact of our constructing the SFHs at  = 7, and not physical.a reminder (c.f.§ 2.2), the simba model actively models the dust content of the galaxies on-the-fly in the simulation but does not model the grain size distribution.As a result, the interstellar extinction curves are assumed to be those of Weingartner & Draine (2001).In contrast, the smuggle model includes as grain size distribution as well, and therefore the radiative transfer used to generate the mock SEDs in Two immediate points are clear from Figure 3. First, the model SEDs that we present here provide a reasonable match to observations3, including SEDs that have relatively blue and relatively red optical spectra.This allows us to proceed in our analysis with reasonable confidence in our methods.Second, the observations exhibit a wide range of rest-frame optical colors, with (as a general statement) the Labbé et al. (2023) sources being redder than the Endsley et al. (2022) galaxies.When comparing against our models, it is immediately evident that the majority of reddening in  > 7 galaxies is due to local obscuration at the sites of very young stars; diffuse dust in the ISM is insufficient to provide the required reddening to match the observed NIRCAM photometry for the reddest sources.This latter point is a net win for galaxy SED fitting: many modern SED fitting codes have the ability to include such clouds in their backwards modeling.This reduces the uncertainties incurred by the diverse shapes of diffuse ISM attenuation curves (Narayanan et al. 2018;Trayford et al. 2020;Lower et al. 2022).

RECOVERING THE STELLAR MASSES OF REDSHIFT > 7 GALAXIES VIA SED FITTING
Having established the nature of early universe star formation histories (at least within the context of two reasonably plausible galaxy formation models), as well as the forward modeled SEDs from these galaxies, we now ask how accurately we can recover the stellar masses of these model galax-  ies via SED fitting the mock SEDs.We fit these SEDs with prospector (Johnson et al. 2019), assuming full JWST NIR-CAM and MIRI coverage of the rest-frame UV-NIR SEDs (noting this is assuming a relatively optimistic wavelength coverage).
By using prospector for the backwards modeling, we are able to minimize the uncertainty incurred by many of the stellar population parameter choices as fundamentally prospector and powderday both use fsps under the hood to model stellar populations.We therefore assume the exact same model IMF, spectral libraries, and stellar isochrone models in order to obviate these potential uncertainties in our SED modeling.We additionally fix the redshift to the true redshift of the galaxy to avoid uncertainties in the redshift fit.We assume a uniform (and arbitrarily chosen) signal to noise ratio of 10 across all bands.We allow the birthcloud model in the SED fits to be flexible.
We assume the non-parametric form for star formation history modeling where the star formation histories are constrained with a set of 3 model priors: 1.The Dirichlet prior (Betancourt & Girolami 2013) is parameterized by concentration index  that sets the preference for all stellar mass to be formed in a single Owing to outshining by recently formed stellar populations, early stellar mass buildup can be missed by SED modeling.This will have a significant impact in the derived stellar masses (c.f. Figure 6).
bin vs a smoother distribution of stellar mass formed over the modeled time period.We have run tests ranging  = [0.3,1.0] and found minimal impacts on our results.
2. The Standard Continuity prior fits for Δlog(SFR) between SFH time bins, thus weighting against dramatic variations in the SFH between adjacent time bins.This weighting is parameterized by a Student's t-distribution (Leja et al. 2019a).This is given by  = log (SFR n /SFR n+1 ), where: (1) Here, Γ is the Gamma function,  controls the width of the distribution, and  controls the tail of the distribution function.Following Leja et al. (2019a), we set  = 2,  = 3.These SFH priors are not meant to be comprehensive, but rather to span a range of reasonable priors commonly used in the literature, and to demonstrate the potential impact of these priors on the derivation of galaxy stellar masses in the early Universe.
In Figure 4, we show the example fit for one of our model SEDs, with the observational filters overlaid (we zoom in to the NIRCAM/MIRI wavelengths).The blue shaded region denotes the 16 − 84% percentile confidence intervals in the posterior, while the black line shows the input powderday mock SED.The fit was performed at  ≈ 7. The quality of fit presented in Figure 4 is comparable to all of the fits performed for this study.5) results in up to ∼order of magnitude uncertainties in the derived stellar masses of high- galaxies with systematic trends in over(under) estimating the  * depending on the nature of the exact prior used.The color coding shows  SFR which is the standard deviation in the star formation rates over the lifetime of the galaxy, and is a measure of the 'burstiness' of the system.The uncertainty is measured as the standard deviation of  * (fit)/ * (true), while the bias is measured as the median of  * (fit)/ * (true).Generally, the simba simulation exhibits less burstiness (c.f. Figure 2).The uncertainties and biases vary dramatically, though the priors that favor rising SFHs (e.g.Wang et al. 2023), as well as those that allow for more dramatic SFH variations (e.g.Suess et al. 2022) tend to perform the best.
In Figure 5, we show the model SFH for an arbitrarily chosen galaxy from both the simba and smuggle simulations with the best fit star formation history as derived from the prospector fits with each of these priors imposed (we show this same plot for all modeled galaxies in the Appendix).In the left column we show the results for an arbitrarily chosen galaxy from the simba simulation, while in the right we show the results from the smuggle model.The orange lines show the median best fit SFH from prospector, while the shaded region shows the inter-quartile dispersion.The simulated SFH is shown in solid blue, and in dashed lighter blue, we bin the simulated SFH at the same time resolution.As is clear none of the imposed SFH priors adequately reproduces the early Universe SFH in either galaxy formation model, though the models that allow for the most dramatic star SFH variations performs the best.
In detail, the Dirichlet prior constrains the fraction of mass formed in each time bin (in the model SFH) such that the fractional sSFR in each time bin follows the Dirichlet distribution.While some variation in the fractional contribution to the total stellar mass from a given SFH bin is allowed, these do not tend to be dramatic with the Dirichlet prior (Leja et al. 2017).As a result, the early SFH for each model is relatively low-level, constant SFR, with the dominant changes happening in the latest time bins, where the majority of the light contributes to the SED.The Continuity prior is similar to the Dirichlet one, though as discussed, favors less dramatic variations in the SFH.The default Continuity prior was tuned by Leja et al. (2019a) to the Illustris-TNG star formation histories which are relatively smooth (owing to the effective equation of state pressurizing the interstellar medium).As a result, the default Continuity prior (as we employ here) will therefore demonstrate even less power on short timescales (Leja et al. 2019a;Lower et al. 2020).The Tacchella et al. (2022) prior modifies this Continuity prior to allow for a burstier SFH, though we still find that these modifications are insufficient to result in the level of burstiness seen in the simba and smuggle simulations.Because of this, the low level nearly constant SFH at early times dominates the stellar mass buildup with the Dirichlet, Continuity, and Tacchella et al. (2022) modifications to the Continuity priors, giving somewhat similar derived stellar masses for all model galaxies.
The rising SFH, developed by Wang et al. (2023), favors rising star formation histories.Generally, the SFHs for both the simba and smuggle models rise with time, with the burstiness of the SFH as the main variant between the two.As a result, the Wang et al. (2023) rising SFH prior performs quite well on the simba galaxy formation model, which have relatively minor bursts, and moderately well on the burstier smuggle SFH.Finally, the Suess et al. ( 2022) SFH prior is a 2-component prior designed to fit the SFH of post starburst galaxies: in particular, it allows for sharp changes in the recent SFH.As a result, while the Suess et al. (2022) prior is designed for post starburst galaxies, it performs reasonably well for our model high- galaxies owing to their highly time-variable SFHs.
We show the impact of these SED fits on the derived stellar masses in Figure 6, where we compare the SED fit  * to the true  * for each galaxy, with panels ordered akin to Figure 5.We color code the galaxies with a measure of their SFH 'burstiness', here parameterized as the standard deviation of the SFR over the history of the Universe.We additionally provide summary statistics for each model via the uncertainty (quantified as the standard deviation of  * (fit)/ * (true)), and the bias (quantified as the median of the same ratio).Smaller uncertainties, and biases that are closer to unity demonstrate higher accuracy in the derived stellar masses.
Generally, no model performs particularly well, with uncertainties including both over-estimates and under-estimates.The lack of early star formation in the SED fits owes to 'outshining': the most recent burst of star formation dominates the SED, and therefore has a significantly outsized impact on the resulting fit.The SFH priors that allow for rising SFHs, as well as those that allow for the most dramatic variations in the SFH quantitatively perform the best (bottom two rows of Figure 6), though neither fully capture the very recent star formation variability.
Whether an SED fit over-predicts or under-predicts the true stellar mass depends in large part on the amount of early star formation that the SED fit predicts.In Figure 7, we show the corner plot for an arbitrarily chosen simba galaxy using the Dirichlet prior.There are significant uncertainties in the derived physical properties as well as co-variances between them, specifically between stellar mass and SFR and the dust attenuation parameters.These degenerate solutions make inferring the true galaxy properties difficult when the available data does not have enough constraining power.
We caution that these results are particular to the bands employed here (NIRCAM+MIRI), as well as the particulars of the SFHs modeled, and are therefore intended to demonstrate the range of uncertainty rather than the specific direction of uncertainty in  * estimates at high-.The take-away from this analysis is not that a particular SFH prior tends to under predict or over predict stellar masses at high-redshift, but rather a demonstration of the relative level of uncertainty in fitting the SEDs of galaxies whose stellar masses are dominated by early star formation that is being outshined by current star formation.

General Discussion
We have, thus far, seen that when observing galaxies at high redshift (here, we model  ≳ 7), SED fitting techniques have difficulty in correctly inferring the stellar masses of galaxies owing to the substantial contribution of stellar mass build up by individual star formation episodes at earlier times.This is true for a range of star formation histories, ranging from the smoother SFHs seen in traditional cosmological simula-tions to more bursty SFHs in explicit feedback models.The fundamental issue is that current star formation at the time of detection is able to outshine the prior stellar mass buildup, making it difficult for SED fitting software to infer its presence.This situation will be most extreme for very bursty systems, though at least within the context of the range of models explored here, appears to be a generic problem.Whether SFHs at high- are bursty in the early Universe remains an open question, though recent observations and models appear to suggest at least some level of burstiness (e.g.Dressler et al. 2023;Endsley et al. 2023;Shen et al. 2023).
As these systems evolve to lower redshift, it is expected that traditional SED fitting techniques will perform substantially better as the ratio of flux from current star formation to that from the integrated stellar mass buildup decreases.Indeed, this was demonstrated explicitly by Lower et al. (2020), who performed experiments similar to those presented in this paper for  = 0 model galaxies in the simba simulation, and found that SED fitting with non-parametric SFHs models could accurately determine galaxy stellar masses.
This said, we note some caution when interpreting the results presented here.The simulations here encompass two different models for stellar feedback that result in varied SFHs, for galaxies drawn from a relatively small box (25/ℎ) Mpc on a side).While these models span a diverse range of SFHs for similar mass galaxies, they are not necessarily comprehensive.It is possible that some model forms of star formation history (i.e.where most of the mass is formed at the time of observation) could result in highly accurate derivations of the stellar mass.Similarly, the modest size of our cosmological boxes excludes the most massive and rare systems.These results should not be taken to mean that a particular SFH prior will always overpredict, or always underpredict stellar masses: the exact relationship between the modeled stellar masses of galaxies and true ones depend on a wide range of choices in SED modeling (e.g.Conroy 2013;Lower et al. 2020;Gilda et al. 2021).Instead, these results are to simply intended to reflect the uncertainty associated with SED modeling early Universe galaxies.

Relationship to other Studies
While we have explicitly demonstrated the issues of deriving stellar masses in high- galaxies via direct numerical simulation, this issue has been hypothesized in the observational literature in a wide range of contexts.For example, Tacchella et al. (2022) studied the stellar populations of 11 galaxies from 9 <  < 11, and found that the inferred stellar ages were significantly impacted by the assumed SFH prior, and noted that multiple priors were able to fit the data equivalently.Similarly, Topping et al. (2022) modeled the stellar masses of UV selected galaxies at  ∼ 7 − 8, and found that the derived stellar masses can be uncertain by up to an order of magnitude owing to the outshining of older stellar populations by a current burst.Whitler et al. (2023) modeled the stellar ages of UV-bright  ∼ 7 galaxies, and found a potential tension between the relatively young inferred stellar ages of their sample, and the detection rate of higher-redshift sources by JWST.This tension can be alleviated if the redshift  ∼ 7 galaxies have older stellar components formed in earlier bursts, as in the simulations presented here.Giménez-Arteaga et al. (2023) fit spatially resolved measurements of 5 <  < 9 galaxies in the JWST SMACS 0723 field with bagpipes (Carnall et al. 2018), and found evidence for both a bursty SFH, as well as potential for an older stellar population (outshined by current star formation).When taking this outshining into account, Giménez-Arteaga et al. 2023 found reduced inferred stellar masses by 0.5 − 1 dex.Finally, Iyer et al. (2019) note, in their development of the dense-basis methodology for nonparametric SFH reconstruction from observed galaxy SEDs, that the ability to model older stellar populations is priordominated rather than likelihood-dominated, and that sharp variations in the SFH may be difficult to infer.
In addition to studying the impact of outshining on stellar mass inference of high- galaxies, we have additionally analyzed the impact of star formation history priors on the derived masses.The inference of star formation histories has been demonstrated by a number of authors to be strongly dependent on the assumed priors for the SFH.For example, Leja et al. (2019a) have derived stellar masses 0.1 − 0.3 dex larger using non-parameric star formation histories as opposed to traditional parametric models in the 3D-HST galaxy catalog.Similarly, studying galaxies at  > 7, Tacchella et al. (2022) and Whitler et al. (2023) found inferred stellar ages and masses to be highly sensitive to form of the assumed SFH in the SED fitting.Finally, Suess et al. (2022) showed that the assumed SFH priors can substantially impact the inferred ages of post starburst galaxies via synthetic stellar population modeling.

Possible Ways Forward
The primary outcome of this Letter is to demonstrate the uncertainty of the measurement of high- stellar masses.We advocate for investment by the community in the development of methods to reduce the bias and uncertainty in these measurements.
One possibility is through the development of new star formation history priors in SED fitting codes.Already we have seen in Figure 6 that priors that allow rapid transitions in star formation history perform reasonably well.In a similar vein, including information from spectral line features that trace star formation over different time scales may help to quantify the burstiness, and inform the modeling of the SFH (e.g.Iyer et al. 2022).
As an alternative to traditional SED fitting techniques, machine learning methods may hold promise.Gilda et al. (2021) demonstrated the efficacy of using mock SEDs from large cosmological simulations as a training set for machine learningbased SED fitting software.Gilda et al. (2021) demonstrated that in some circumstances, these techniques can far outperform traditional SED fitting.

SUMMARY AND OUTLOOK
In this paper, we have employed cosmological galaxy evolution simulations in order to investigate the ability for modern SED fitting techniques to recover the stellar masses of highredshift galaxies.Owing to the relatively young ages of  ≳ 7 galaxies, early bursts of star formation constitute a significant fraction of the formed mass at the time of SED modeling.As a result, if the stellar light from this early star formation (c.f. Figure 5) is outshined by late-time star formation at the time of detection, then the recovered stellar masses can be incorrect by nearly an order of magnitude (Figure 6).The impact of these uncertainties will decrease at lower redshifts, as the ratio of the light emitted from star formation at the time of detection to the the light from the underlying stellar mass decreases.As JWST pushes the frontier of high-redshift science to increasingly early times, we encourage the development of new techniques in order to accurately derive the physical properties of these first galaxies.

APPENDIX
Here, we show the model SFH for all of our simulated galaxies, as well as the prospector SED fits (for every modeled SFH prior) for each galaxy.The galaxies are ordered by mass (i.e. for an individual plot, the simba galaxy and smuggle galaxy each represent the  th most massive galaxy in the cosmological volume.

Figure 1 .
Figure 1.Histogram of the stellar masses at  ∼ 7.5 for the 20 most massive galaxies in the simba cosmological simulation (red) and the smuggle cosmological box (blue).Both boxes are 25/ℎ Mpc on a side.

Figure 2 .
Figure2.Example star formation histories of 5 arbitrarily chosen galaxies from the simba and smuggle cosmological boxes.The left panel shows simba, while the right panel shows smuggle.The former has an artificially pressurized ISM (as do most cosmological galaxy formation models), while the right represents a class of explicit feedback models.The two models highlight the diversity of modeled star formation history burstiness possible in galaxy formation models.We use both of these models in our analysis to model a range of possibilities when investigating the efficacy of SED fitting in young galaxies.Note that the sudden drop off in the SFH at  < 7 is an artifact of our constructing the SFHs at  = 7, and not physical.
Figure 2 includes the spatially varying dust extinction curves.The simba model uses a Milky Way template for PAH emission (scaled for the local energy deposited), while the smuggle model uses the Draine et al. (2021) model, based on the local grain populations (Narayanan et al. 2023), resulting in fairly different mid-IR emission features.Therefore, as in modeling the SFHs, the two dust models that we include here bracket a reasonable range of modeled dust extinction and obscuration.In the bottom row of Figure 3, we show the same SEDs, but this time include Charlot & Fall (2000) "birthclouds" around young stars.These birthclouds are included in a subresolution fashion such that any star particle with age  age < 10 7 yr experiences an extra attenuation of  V =1.We show observational comparisons from the Endsley et al. (2022)2 and Labbé et al. (2023) surveys in the top and bottom, respectively (these are split in the top and bottom for clarity).

Figure 3 .
Figure 3. Modeled SEDs for the simba (left) and smuggle (right) models at  ∼ 7.5.The top row shows the results from powderday radiative transfer when only the diffuse dust is included, while the bottom row includes Charlot & Fall (2000) birthclouds with an  V = 1 screen in front of all  < 10 Myr star particles.The green points show observations from Endsley et al. (2022), while grey points show observations fromLabbé et al. (2023).We find that birthclouds dominate the obscuration in the UV/optical.

Figure 4 .
Figure 4. Comparison of SED fit to input SED.The input mock galaxy SED generated by powderday is in black, while the 16−84 th percentile confidence intervals in the fit are denoted by the blue shaded region.The orange lines show the NIRCAM and MIRI filter locations for the fit, and the SEDs are redshifted to  ≈ 7.

Figure 5 .
Figure 5. Star formation histories for an arbitrarily chosen galaxies from the simba simulation (left), and smuggle simulation (right), with the best fit derived star formation histories from SED modeling overlaid (orange).The simulated SFH binned at the same time resolution as the SED fits is shown in a lighter, dashed-blue line.Owing to outshining by recently formed stellar populations, early stellar mass buildup can be missed by SED modeling.This will have a significant impact in the derived stellar masses (c.f.Figure6).
3. The Bursty Continuity prior, developed byTacchella et al. (2022) adjusts the parameters in the Student's t-distribution to  = 1 and  = 2, resulting in a burstier SFH.4.Rising SFH prior developed byWang et al. (2023) to preferentially favor rising star formation histories.5.PSB prior developed bySuess et al. (2022) for post starburst galaxies which have sharp changes in their recent SFHs.

Figure 6 .
Figure 6.Comparison of the best fit derived  * from SED fitting to the true  * for the simba model (left) and the smuggle model (right).The result of poorly fit star formation histories (Figure5) results in up to ∼order of magnitude uncertainties in the derived stellar masses of high- galaxies with systematic trends in over(under) estimating the  * depending on the nature of the exact prior used.The color coding shows  SFR which is the standard deviation in the star formation rates over the lifetime of the galaxy, and is a measure of the 'burstiness' of the system.The uncertainty is measured as the standard deviation of  * (fit)/ * (true), while the bias is measured as the median of  * (fit)/ * (true).Generally, the simba simulation exhibits less burstiness (c.f.Figure2).The uncertainties and biases vary dramatically, though the priors that favor rising SFHs (e.g.Wang et al. 2023), as well as those that allow for more dramatic SFH variations (e.g.Suess et al. 2022) tend to perform the best.

Figure 7 .
Figure 7. Example prospector model parameter posteriors for an arbitrarily chosen simba galaxy fit with the Dirichlet prior SFH model.The true physical quantities are: log(M * /M ⊙ )=9.2; log(SFR/M ⊙ yr −1 )=0.85; log(Z/Z ⊙ )) = -1.8;The stellar mass is moderately degenerate with the dust attenuation parameters while the SFR is highly degenerate with the birth cloud dust component.The stellar metallicity posterior is bimodal, so is not well constrained by the data.

Figure 8 .
Figure 8. Star formation histories for the 1 st most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 10 .
Figure 10.Star formation histories for the 3 rd most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 11 .
Figure 11.Star formation histories for the 4 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 12 .
Figure 12.Star formation histories for the 5 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 13 .
Figure 13.Star formation histories for the 6 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 14 .
Figure 14.Star formation histories for the 7 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 15 .
Figure 15.Star formation histories for the 8 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 16 .
Figure 16.Star formation histories for the 9 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 17 .
Figure 17.Star formation histories for the 10 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 18 .
Figure 18.Star formation histories for the 11 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 19 .
Figure 19.Star formation histories for the 12 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 20 .
Figure 20.Star formation histories for the 13 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 21 .
Figure 21.Star formation histories for the 14 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 22 .
Figure 22.Star formation histories for the 15 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 23 .
Figure 23.Star formation histories for the 16 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 24 .
Figure 24.Star formation histories for the 17 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 25 .
Figure 25.Star formation histories for the 18 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 26 .
Figure 26.Star formation histories for the 19 th most massive galaxy in the simba (left) and smuggle (right) simulation.

Figure 27 .
Figure 27.Star formation histories for the 20 th most massive galaxy in the simba (left) and smuggle (right) simulation.