Can Cosmological Simulations Reproduce the Spectroscopically Confirmed Galaxies Seen at $z\geq 10$?

Recent photometric detections of extreme $(z>10)$ redshift galaxies from the JWST have been shown to be in strong tension with existing simulation models for galaxy formation, and in the most acute case, in tension with $\Lambda CDM$ itself. These results, however, all rest on the confirmation of these distances by spectroscopy. Recently, the JADES survey has detected the most distant galaxies with spectroscopically confirmed redshifts, with four galaxies found with redshifts between $z=10.38$ and $z=13.2$. In this paper, we compare simulation predictions from four large cosmological volumes and two zoom-in protoclusters with the JADES observations to determine whether these spectroscopically confirmed galaxy detections are in tension with existing models for galaxy formation, or with $\Lambda CDM$ more broadly. We find that existing models for cosmological galaxy formation can generally reproduce the observations for JADES, in terms of galaxy stellar masses, star formation rates, and the number density of galaxies at $z>10$.


INTRODUCTION
The successful deployment of the JWST has already produced observations of the highest-redshift galaxies detected to date. The first sets of detections reported Labbe et al. 2022;Naidu et al. 2022a;Adams et al. 2023) have found galaxies with M * > 10 8 M at z phot > 10. The most extreme of these examples, with M * > 10 10 M at z phot > 10 have already been shown to be in strong tension with ΛCDM (Haslbauer et al. 2022;Boylan-Kolchin 2022). These tensions have, however, been clouded by the large uncertainty in fitting photometric redshifts at such extreme distances (Bouwens et al. 2022;Naidu et al. 2022b;Kaasinen et al. 2022). Naidu et al. (2022b) and Zavala et al. (2022) found that breaks in the Spectral Energy Distribution (SED) produced by dust obscuration at z ∼ 5 can masquerade as a Lyman break at z ∼ 17 for recent JWST observations (see also (Fujimoto et al. 2022) for a similar effect in ALMA observations).

bkeller1@memphis.edu
As Boylan-Kolchin (2022) showed, if the comoving number density of galaxies at z = 10 with M * > 10 10 M is as high as can be estimated from Labbe et al. (2022), it would be a significant challenge for ΛCDM : analogous to Haldane's "fossil rabbits in the Precambrian" (Harvey et al. 1996). Haslbauer et al. (2022) has shown that the observations of galaxies with high stellar mass (M * > 10 9 ) at high redshift z 10 from Adams et al. (2023), Labbe et al. (2022), Naidu et al. (2022a), and Naidu et al. (2022b) are in extreme tension with the simulation predictions of EA-GLE , TNG50, and TNG100 (Pillepich et al. 2018a;Nelson et al. 2019a). These authors warn however that spectroscopic confirmation of redshifts is needed before final conclusions can be reached: these tensions all rest on the estimations of stellar masses provided by SED fitting and, critically, on the distance estimates themselves. Without spectroscopic confirmation of the redshifts reported in these observations, these potential tensions may be illusory: artifacts of overestimated distances, and thus overestimated intrinsic luminosities. Indeed, as Behroozi et al. (2020) has shown, semi-empirical modelling of galaxy for-mation in ΛCDM predicts JWST-detectable galaxies with M * > 10 7 M to at least z ∼ 13.5.
In this letter, we compare the observed JADES galaxies to predictions from an array of large cosmological hydrodynamical simulations. These simulations reproduce observed galaxy population statistics (such as the observed galaxy stellar mass function and fundamental plane of star formation) at low redshift. With JADES, we are able to test whether those same models for galaxy formation fail to reproduce these new observations of high redshift galaxy formation.

SIMULATION DATA
In order to probe the predictions of current galaxy formation models, we examine simulation data from the EAGLE Schaye et al. 2015;McAlpine et al. 2016), Illustris (Vogelsberger et al. 2014a,b;Genel et al. 2014;Nelson et al. 2015), TNG100 (Pillepich et al. 2018b;Springel et al. 2018;Nelson et al. 2018;Naiman et al. 2018;Marinacci et al. 2018;Nelson et al. 2019b), and Simba (Davé et al. 2019) cosmological volumes. Each of these simulations has a volume of ∼ 10 6 cMpc 3 , with baryonic mass resolution of 10 6 M − 10 7 M , which allows them to (marginally) resolve the formation of M * = 10 8 M − 10 9 M galaxies (an M * = 5 × 10 8 M galaxy in EAGLE, Illustris, TNG100, and Simba will contain 276, 384, 357, and 27 star particles respectively). Each of these projects includes models for gas cooling, star formation, and feedback from supernovae (SNe) and active galactic nuclei (AGN) that are tuned to reproduce the z ∼ 0 stellar mass function (among other low-redshift population statistics). In Table 1 we list the cosmological parameters used for each simulation, as well as the size of the simulation volume, the range and number of snapshots with redshift between z ∼ 10 and z ∼ 14, and the mass resolution for dark matter (DM) and baryonic particles. We rely on the public releases of the halo catalogs from each simulation to determine galaxy stellar masses and star formation rates in these simulation data sets. We have also included data from the cosmological zoomin simulations OBELISK (Trebitsch et al. 2021) and Romu-lusC (Tremmel et al. 2019). OBELISK is a zoom-in of the most massive halo at z = 2 in the HORIZON-AGN volume (Dubois et al. 2014). This region is selected to include all DM particles within 4R vir of the halo center at z = 2. At z = 0, this cluster has a halo mass of M halo ∼ 6.6×10 14 M . RomulusC is a zoom-in simulation of a smaller galaxy cluster, with halo mass M halo = 1.5 × 10 14 M at z = 0. It is drawn from a 50 3 cMpc 3 DM-only volume, and applies the same modelling approach for cooling, star formation, supermassive black hole (SMBH) evolution, and feedback as the Romulus volume (Tremmel et al. 2017). Unlike the other data sets we examine here, these zooms do not include a full resolution sample of a large volume; instead they offer an higher-resolution picture of early collapsing overdensities.

RESULTS
We begin by simply showing the distribution of galaxy stellar masses in each simulation volume as a function of redshift, shown in Figure 1. Not all of the simulations we examine here have well-sampled snapshots above z = 10 (EAGLE in particular only includes one snapshot between z ∼ 10 and z ∼ 14). However, as can be seen, all of the simulations produce a large number of galaxies with M * > 10 8 M at z ∼ 10. At all redshifts, Simba and OBELISK produce the most massive galaxies of the simulations we examine, in part due to the larger volume they simulate/are drawn from (∼ 2.4 times the comoving volume of TNG100, the next largest vol-  Table 1. Simulation data compared in this study. Each cosmological volume has a box of side length ∼ 100 cMpc. We show the choice of cosmological parameters, box size, redshift range between z ∼ 10 − 14, number of snapshot outputs in that range, and resolution for DM and baryons. The ( † ) denotes that these are zoom-in simulations of an individual galaxy protocluster, drawn from a lower-resolution volume. OBELISK is an Eulerian simulation, and thus the baryonic resolution reported here is given as the typical mass of a star particle and the mass of a gas cell with ∆x = 35 pc at a density of n > 10 cm −3 . ume), as well as differences in the choice of subgrid physics models. As we move to higher redshift, only Simba and OBELISK produce even a single galaxy above the best-fit stellar masses of JADES-GS-z11-0, JADES-GS-z12-0, and JADES-GS-z13-0. Illustris and TNG100 are unable to produce galaxies at z ∼ 11.5 − 12 which reach even the lower estimate for the stellar mass of JADES-GS-z11-0.
As Robertson et al. (2022) has also provided estimates for the star formation rates (SFR) in the 4 JADES observations, we show in Figure 2 the SFR-M * plane for simulated galaxies at z ∼ 10. All simulation volumes produce galaxies at z ∼ 10 with sSF R = 10 Gyr −1 , except for RomulusC, which has an sSF R roughly half this value. An sSF R of 10 Gyr −1 is approximately 4 times higher than the star forming main sequence at z = 2 (Rodighiero et al. 2011), consistent with the observed trend of increasing sSFR towards higher z ∼ 6 redshift (Santini et al. 2017). The banding seen in the lower-mass Simba galaxies is simply a function of resolution, as these galaxies contain only a handful of star particles. Interestingly, there appears to be no noticeable trend in the SFR as a function of stellar mass for the JADES observations. However, drawing strong conclusions as to the nature of the SFR−M * relation at z ≥ 10 is ill-advised given the small number of observations, their relatively large uncertainties, and the intrinsic scatter in the SFR-M * plane that the simulations predict. Beyond these, the NIRSpec observations of the JADES galaxies are selected from photometric observations, which are subject to a continuum flux (and therefore SFR) limit, introducing a potential observational bias. Overall, the SFRs measured from the simulations match the JADES observations reasonably well, though with perhaps higher SFRs for simulated galaxies at masses above M * ∼ 5 × 10 8 M . Further observations will help reveal if the flat SF R − M * curve is a real feature of z ≥ 10 galaxy formation, or simply an artifact of the small sample size of the spectroscopically confirmed JADES galaxies and/or uncertainties in the estimation of M * or SF R.
In order to make a more quantitative estimate of potential tension between simulated cosmological volumes and the Robertson et al. (2022) observations, we now look to the comoving number density of galaxies at each of the redshifts probed by JADES. The deep observing area of JADES, using the JWST NIRCam, (9.7 square arcmin), yields a volume of V ∼ (9 × 9 × 494) cMpc 3 ∼ 4 × 10 4 cMpc 3 bounded by the comoving distance of 494 Mpc between the JADES-GS-z10-0 at z = 10.38 and JADES-GS-z13-0 at z = 13.2. These candidates were, however, originally selected photometrically from a wider field of 65 arcmin 2 , which yields a volume of V ∼ 2.7 × 10 5 cMpc 3 . We can thus estimate the number density of the JADES galaxy observations as n ∼ 3.7 × 10 −6 cMpc −3 . In Figure 3, we show the comoving number density of galaxies above a given stellar mass for simulation snapshots nearest in redshift to each of the JADES observations (if the nearest snapshot is separated by more than 0.5 in redshift, we omit it from the plot). Each of the JADES observations above z > 10.38 implies For z ∼ 10, we also show the independent measurements for the galaxy stellar mass function measured by the GREATS ALMA survey by Stefanon et al. (2021). The highest predicted number density for massive halos at these redshifts is from Simba, and only JADES-GS-z11-0 and JADES-GS-z12-0 imply number densities higher than in Simba for galaxies at their stellar mass. Even for the most extreme case, JADES-GS-z11-0, the probability of finding a galaxy with stellar mass above M * = 10 8.9 M at z ∼ 11, in a random volume of 2.7 × 10 5 cMpc 3 in Simba, is 8%. If we instead take the lower estimate for the stellar masses from JADES-GS-z11-0 (M * = 10 8.5 M ), this probability rises to 92%. There does not appear to be any significant tension in the density of galaxies with M * ∼ 10 8 M at z > 10 inferred from JADES and the number density produced by at least one existing cosmological simulation (Simba).

DISCUSSION
We have compared the recent detection of spectroscopically confirmed z > 10 galaxies from the JADES survey to the EAGLE, Illustris, TNG100, and Simba surveys. We show that all four simulation suites produce galaxies with M * ∼ 10 8.5 M at z ∼ 10, and that Simba in particular produces at least one galaxy with mass comparable to the JADES observations at each observed redshift. The JADES galaxies show a relatively flat SF R − M * trend, in contrast to the constant sSFR of ∼ 5 − 10 Gyr −1 for simulated galaxies at z ∼ 10. Comparing the estimated JADES galaxy stellar mass functions to the simulation predictions shows that these simulations are in reasonable agreement with the number density of galaxies at z ∼ 10 observed by JADES and the earlier ALMA measurements from Stefanon et al. (2021) (though Illustris and TNG100 appear to be slightly under-predicting the formation of galaxies at z ∼ 10). At z ∼ 11 − 12, JADES implies a slightly higher number density of galaxies at M > 10 8 M than are predicted by any of the simulation volumes we examine here. We find, however, that given the small number of objects confirmed in JADES (in a volume of ∼ 2.7 × 10 5 cMpc 3 ), none of the JADES observations is in greater than 2σ tension for the best-fit stellar mass, relative to the Simba simulation. JADES-GS-z11-0 implies a slightly higher density at z = 11 for galaxies with M * ∼ 10 9 M , but the difference between Simba and density inferred from JADES is only slight. A randomly chosen volume of 2.7 × 10 5 cMpc 3 from Simba at z = 11.4 will contain a galaxy with stellar mass greater than that of JADES-GS-z11-0 8% of the time. At the lower uncertainty for the JADES-GS-z11-0 stellar mass, this probability grows to 92%.
As future observations better constrain the stellar mass function at z ∼ 11, this tension may strengthen or disappear. If objects as massive as those seen by (Labbe et al. 2022;Adams et al. 2023;Naidu et al. 2022a) are confirmed to be common, with spectroscopic redshifts of z > 10, this would imply a serious problem with existing ΛCDM and galaxy formation theory (Boylan-Kolchin 2022; Haslbauer et al. 2022;Lovell et al. 2023). In order to be in tension with all models for galaxy formation in ΛCDM (by implying galaxy star formation efficiencies > 100%), the number density of M * = 10 9 M galaxies at z ∼ 11 would need to be more than 3 orders of magnitude higher than what is implied by the JADES-GS-z11-0 detection. This tension may be resolved by future refinements in the SED-fit stellar mass estimates lowering the stellar mass measured for JADES-GS-z11-0, or simply by JADES-GS-z11-0 residing in a moderately rare (P ∼ 8%) overdensity. For a mean cosmological mass density of ρ m ∼ 4×10 10 M cMpc −3 , a comoving volume of 4 × 10 4 cMpc 3 contains a total mass of ∼ 10 15 M . If the entire volume covered by the deepest JADES observations is a rare overdensity that will eventually collapse to form a z ∼ 0 cluster, it will be on the order of the most massive clusters detected (Lovell et al. 2023). It is unlikely that the entire JADES volume is probing a single large overdensity (given the fact that it is a pencil-beam volume with a comoving depth of 494 cMpc over an area on the sky of only (9 × 9) cMpc 2 ), but the possibility that one or more of the JADES galaxies above z > 10 resides within one or more protoclusters such as those simulated in OBELISK is still plausible.
One question that arises from this data is why the Simba simulations predict more M * > 10 8 M galaxies at z > 10 than Illustris and TNG100. While the Simba volume is somewhat larger than Illustris and TNG100 (Simba's volume is ∼ 2.4 times larger than TNG100), it also produces a higher density of galaxies with M * > 10 8 M at all redshifts we have examined here. Each of these simulations applies a different set of models for star formation and feedback by SNe and AGN. In Simba, SMBHs are seeded in galaxies with M * > 10 9.5 M , while in TNG100 SMBHs are seeded in halos with M halo > 7.4×10 10 M . This means that none of the galaxies we show in Figure 1 or in Figure 3 contain SMBHs, so the differences we find must be a function of the different star formation and feedback recipes used in TNG100 and Simba. Simba applies a tuned reduction in the SN-driven outflow mass loadings η above z > 3 to better match observations of the galaxy stellar mass function at z > 6, lowering the mass loadings by (a/0.25) 2 (Davé et al. 2019). It appears that lower mass loadings at high redshift are not only important to matching observations at z ∼ 6, but for z > 10 as well. Understanding how SNe regulate galaxy formation and drive outflows in these early epochs will be an important avenue of future simulation study.
While the results we have shown here show that there is no strong tension between at least one existing large-volume cosmological simulation and the spectroscopically confirmed galaxy detections from JADES, it is important to note that the cosmological volumes we have examined here are all relatively low resolution: even a 10 9 M galaxy only contains ∼ 50 star particles in Simba. These simulations are also tuned to reproduce z ∼ 0 galaxy properties. Star formation and feedback at low redshift are still relatively poorly understood processes (Naab & Ostriker 2017), a problem that becomes much more severe at epochs as early as z > 10 ( Visbal et al. 2020).
An obvious direction for future studies is to search for galaxies of similar mass in higher resolution, high redshift simulation volumes. In particular, simulations designed to study reionization such as Renaissance (O'Shea et al. 2015), OBELISK (Trebitsch et al. 2021), and SPHINX (Rosdahl et al. 2018), all achieve resolutions significantly higher than any of the volumes we have studied here. Many also feature more physically motivated models for stellar feedback, which are possible at these higher resolutions. Even at the same resolution as these (∼ 100Mpc) 3 volumes, rare regions can be probed by zooming in on overdensities from much larger volumes, a strategy applied by FLARES (Lovell et al. 2021). FLARES produces a stellar mass function at z ∼ 10 consistent with Stefanon et al. (2021), which we show here to be consistent with JADES. Zoom-in simulations of protoclusters such as OBELISK will be particularly fruitful, as the earliest bright galaxies to form are likely to be found in these environments. As we show in Figure 1, OBELISK contains a number of simulated galaxies with comparable stellar masses to each of the JADES detections at every redshift of the JADES sample, resolved with > 1000 star particles. Zoom-in simulations such as these will be a powerful tool for understanding the earliest phase of galaxy formation at z > 10.

CONCLUSIONS
We have compared the recent observations of the highest redshift galaxies with spectroscopically confirmed distances from the JADES survey (Robertson et al. 2022) to simulated galaxies from EAGLE, Illustris, TNG100, and Simba volumes and the OBELISK and RomulusC zoom-ins. In general, we find that each of these simulations produces galaxies with comparable stellar masses to the JADES galaxies by z ∼ 10. The most massive JADES galaxies have somewhat lower SFRs than simulated galaxies at z ∼ 10, but lie within the scatter of the simulations. The galaxy number density implied by the JADES galaxies at z ∼ 10 is consistent with both the simulations and past observations. At higher redshift, only Simba and OBELISK produce galaxies as massive as are found in JADES. The number density of galaxies inferred from JADES is slightly larger than what is predicted by Simba at z = 11 and z = 12, but at a low level of signifi-cance. Overall, there appears to be no strong tension between models for galaxy formation in cosmological hydrodynamic simulations and the most distant spectroscopically confirmed galaxies.