Has JWST Already Falsified Dark-matter-driven Galaxy Formation?

The formation of the first galaxies in the observed universe is a key question in modern astrophysics and one of the most important science goals of the recently launched James Webb Space Telescope (JWST). The Near Infrared Camera instrument (NIRCam; Rieke et al. 2005) of JWST observes the universe in the ≈0.6–5 μm regime. This allows the detection of objects at redshifts z ≳ 12, thus, revealing the evolutionary stage of galaxies in the early universe just ≈400 Myr after the Big Bang. The most distant confirmed galaxy is GN-z11 with a stellar mass of M_* ≈ 10⁹ M_⊙ at a spectroscopic redshift of ${z}_{\mathrm{spec}}={11.09}_{-0.12}^{+0.08}$ detected with the Hubble Space Telescope (Oesch et al. 2016). Recently, Naidu et al. (2022b) reported the discovery of two luminous galaxy candidates labeled as GLASS-z11 (hereafter GL-z11) and GLASS-z13 (GL-z13) with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })={9.4}_{-0.3}^{+0.3}$ and ${9.0}_{-0.4}^{+0.3}$ located at photometric redshifts of ${z}_{\mathrm{phot}}={10.9}_{-0.4}^{+0.5}$ and ${z}_{\mathrm{phot}}={13.1}_{-0.7}^{+0.8}$, respectively (see their Table 3). In a subsequent publication, Naidu et al. (2022a) presented the galaxy candidate CEERS-1749 with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })={9.6}_{-0.2}^{+0.2}$ most likely located at ${z}_{\mathrm{phot}}={16.0}_{-0.6}^{+0.6}$ but a secondary redshift solution of z ≈ 5 cannot be excluded (see their Figure 1 and Table 3).

Further luminous high-redshift galaxy candidates have been discovered by JWST (e.g., Adams et al. 2022; Atek et al. 2022; Furtak et al. 2022; Harikane et al. 2022; Yan et al. 2022; Labbe et al. 2022). For example, Adams et al. (2022) studied the properties of candidates over a redshift range of 9 < z_phot < 12 and the source with ID 1514 has ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })={9.8}_{-0.2}^{+0.2}$ at ${z}_{\mathrm{phot}}={9.85}_{-0.12}^{+0.18}$ (see their Tables 3 and 5). Labbe et al. (2022) identified candidates with M_* > 10¹⁰ M_⊙ over 7 < z_phot < 11 from which the two most massive are ID 14924 with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=10.93$ at z_phot = 9.92 and ID 38094 with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=11.16$ at z_phot = 7.56 (see their Figure 3). The spectroscopical confirmation of these objects is still outstanding because the high photometric redshifts of these galaxy candidates can, e.g., artificially emerge due to dust attenuation (Naidu et al. 2022a, 2022b; Zavala et al. 2022).

In this contribution, we aim to investigate if ID 1514, ID 14924, GL-z11, GL-z13, and CEERS-1749 are consistent with the hierarchical buildup of stellar mass as predicted by the ΛCDM paradigm (Efstathiou et al. 1990; Ostriker & Steinhardt 1995) using the Illustris The Next Generation (TNG; Pillepich et al. 2018a; Nelson et al. 2019b; Pillepich et al. 2019) and Evolution and Assembly of GaLaxies and their Environments (EAGLE; Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016) projects.

The IllustrisTNG (Pillepich et al. 2018a, 2018b; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019a, 2019b; Springel et al. 2018; Pillepich et al. 2019) and EAGLE (Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016) projects consist of a suite of hydrodynamical cosmological simulation runs conducted in the ΛCDM framework.

The IllustrisTNG project assumes a Planck-2015 (Planck Collaboration XIII 2016) cosmology with the cosmological parameters being H₀ = 67.74 km s⁻¹ Mpc⁻¹, Ω_b,0 = 0.0486, Ω_m,0 = 0.3089, Ω_Λ,0 = 0.6911, σ₈ = 0.8159, and n_s = 0.9667. The simulation runs are based on the moving-mesh code AREPO (Springel 2010) and self-consistently evolve the gas cells, stellar, black hole, and dark matter particles from redshift z = 127 up to present time. The Subfind Subhalos ⁵ and Group ⁶ catalogs are available for 100 different time steps (snapshots) in the redshift range of 0 ≤ z ≤ 20.05. Here, we analyze the snapshot at redshift z = 14.99 (corresponding to an age of the universe of t = 0.271 Gyr and a snapshot number of snapnum = 1), at z = 11.98 (t = 0.370 Gyr, snapnum = 2), at z = 10.98 (t = 0.418 Gyr, snapnum = 3), and at z = 10.00 (t = 0.475 Gyr, snapnum = 4) of the high-resolution realization TNG50-1 and TNG100-1. The former has a box side of 35 h⁻¹ = 51.7 comoving Mpc (cMpc), respectively, where h is the present Hubble constant H₀ in units of 100 km s⁻¹ Mpc⁻¹, a baryonic element mass of m_b = 8.5 × 10⁴ M_⊙, and a dark matter particle mass of m_dm = 4.5 × 10⁵ M_⊙. TNG100-1 has, with 75 h⁻¹ = 110.7 cMpc, a larger box size and with m_b =1.4 × 10⁶ M_⊙, and m_dm = 7.5 × 10⁶ M_⊙, a lower resolution than TNG50-1 (see, e.g., also Table 1 of Nelson et al. 2019b).

The EAGLE project is consistent with the Planck-2013 (Planck Collaboration I 2014) cosmology being H₀ =67.77 km s⁻¹ Mpc⁻¹, Ω_b,0 = 0.04825, Ω_m,0 = 0.307, Ω_Λ,0 =0.693, σ₈ = 0.8288, and n_s = 0.9611 (see also Table 1 of Schaye et al. 2015). Its simulations run with a modification of the GADGET-3 smoothed particle hydrodynamics code (e.g., Springel 2005) starting also at z = 127 and self-consistently evolving the baryonic and dark matter particles up to the present day. The publicly available subhalo catalogs (Table B.1 of McAlpine et al. 2016) are recorded for 29 snapshots in the redshift range of 0 ≤ z ≤ 20.00 (see Table C.1 of McAlpine et al. 2016). We use the two high-resolution realization runs RefL0025N0752 and RecalL0025N0752, and the two lower resolution runs RefL0050N0752 and RefL0100N1504 at z =15.13 (snapnum = 1) and z = 9.99 (snapnum = 2). The two high-resolution runs have a box size of 25 cMpc with an initial baryonic particle mass of m_b = 2.26 × 10⁵ M_⊙, and a dark matter particle mass of m_dm = 1.21 × 10⁶ M_⊙. RefL0050N0752 and RefL0100N1504 have a size of 50 cMpc and 100 cMpc, respectively, and both have an initial baryonic particle mass of m_b = 1.81 × 10⁶ M_⊙ and a dark matter particle mass of m_dm = 9.70 × 10⁶ M_⊙ (Table 2 of Schaye et al. 2015).

Table 1. Comparison of Several Observed Galaxy Candidates with ΛCDM Simulations

	Redshift	Stellar Mass
	(z)	${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })$
CEERS-1749	${16.0}_{-0.6}^{+0.6}$	${9.6}_{-0.2}^{+0.2}$ (8.53) a
GL-z13	${13.1}_{-0.7}^{+0.8}$	${9.0}_{-0.4}^{+0.3}$ (8.01) a
GL-z11	${10.9}_{-0.4}^{+0.5}$	${9.4}_{-0.3}^{+0.3}$ (8.13) a
ID 1514	${9.85}_{-0.12}^{+0.18}$	${9.8}_{-0.2}^{+0.2}$ (8.53) a
TNG50-1	14.99	6.87
TNG50-1	11.98	7.64
TNG50-1	10.98	7.97
TNG50-1	10.00	8.33
TNG100-1	14.99	7.32
TNG100-1	11.98	8.07
TNG100-1	10.98	8.43
TNG100-1	10.00	8.73
RefL0025N0752	15.13	7.05
RefL0025N0752	9.99	7.59
RecalL0025N0752	15.13	6.91
RecalL0025N0752	9.99	7.78
RefL0050N0752	15.13	7.58
RefL0050N0752	9.99	8.76
RefL0100N1504	15.13	7.70
RefL0100N1504	9.99	9.06

Notes. The first four rows show the photometric redshifts and stellar masses of the observed galaxy candidates ID 1514 (see Tables 3 and 5 of Adams et al. 2022), GL-z11, GL-z13 (see Table 3 of Naidu et al. 2022b), and CEERS-1749 (see Table 3 of Naidu et al. 2022a). The other rows list the maximum stellar mass of a subhalo at a given redshift for the TNG and EAGLE runs.

^a Lowest possible stellar mass value for an invariant canonical IMF as quantified in Section 3.1.

Download table as:  ASCIITypeset image

The galaxy stellar mass function (GSMF) at redshifts z = 14.99, 11.98, 10.98, and 10.00 in the TNG runs and at z = 15.13 and 9.99 in the EAGLE runs are presented in Figure 1. The global peak of the distribution depends on the resolution and/or box size of the simulation runs such that the formation of low massive galaxies depends on the particle resolution but also because small simulation boxes lack large-scale density fluctuations. As a consequence, not-large-enough simulation boxes would not allow the formation of large galaxy clusters, therefore hampering the growth of central (but also noncentral) galaxies. Thus, we mainly focus on the larger simulation boxes TNG100-1 and RefL0100N1504.

Zoom In Zoom Out Reset image size

In the following, we compare the stellar mass buildup as predicted by the ΛCDM simulations with the masses of the observed high-redshift galaxy candidates ID 1514 (Adams et al. 2022), ID 14924 (Labbe et al. 2022) ⁷ , GL-z11, GL-z13 (Naidu et al. 2022b), and CEERS-1749 (Naidu et al. 2022a). For the TNG runs, ID 1514 located at ${z}_{\mathrm{phot}}={9.85}_{-0.12}^{+0.18}$ and ID 14924 at z_phot = 9.92 are compared with simulated galaxies at z = 10.00. GL-z11 at ${10.9}_{-0.4}^{+0.5}$, GL-z13 at ${13.1}_{-0.7}^{+0.8}$, and CEERS-1749 at ${16.0}_{-0.6}^{+0.6}$ are compared with the simulations at z = 10.98, 11.98, and 14.99, respectively. The comparison with GL-z13 and CEERS-1749 is therewith more conservative because snapshots corresponding to lower redshifts than observed are addressed allowing the galaxies to grow in stars for a longer time span than in the observed cases.

In TNG100-1, the GSMF at z = 14.99 reaches a maximum value of ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=7.32$ being therewith ≈ 191 times lower than the stellar mass of CEERS-1749 with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\,={9.6}_{-0.2}^{+0.2}$. At z = 11.98 and z = 10.98, the maximum stellar masses of subhalos are ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=8.07$ and 8.43, respectively, which is less massive than GL-z13 with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })={9.0}_{-0.4}^{+0.3}$ and GL-z11 with ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })={9.4}_{-0.3}^{+0.3}$. The maximum value of ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=8.73$ at z = 10.00 is significantly lower than the stellar mass of ID 1514, which has ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })={9.8}_{-0.2}^{+0.2}$. The discrepancy becomes even more significant for ID 14924, which has ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=10.93$ at z_phot = 9.92. The above stellar masses refer to all stellar particles bound to the considered subhalo depending therewith on the subhalo-finding algorithm. The identification of subhalos can be disturbed, e.g., by merger events, which frequently occur especially at high redshifts. Therefore, we also assess the maximum stellar masses of halos, which accounts for both the fact that the subhalo finder can split a galaxy in clumps underestimating the total mass of the galaxy, and for the inclusion of observationally unresolved satellite galaxies. The maximum stellar masses of halos are ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })=7.32$ (7.05), 8.10 (7.88), 8.43 (8.24), and 8.78 (8.59) at z = 14.99, 11.98, 10.98, and 10.00, in the TNG100-1 (TNG50-1) run, respectively. Thus, using the most massive halo instead of the most massive subhalo in terms of its stellar mass does not significantly effect the results of the TNG runs.

In the EAGLE runs, ID 1514/ID 14924 and CEERS-1749 are compared with the GSMF at z = 9.99 and 15.13, respectively. Unfortunately, the EAGLE database does not list snapshots that match the observed redshifts of GL-z11 and GLz-13. Thus, the EAGLE analysis only focuses on ID 1514, ID 14924, and CEERS-1749.

The RefL0050N0752 and RefL0100N1504 snapshots contain galaxies reaching up to ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\approx 7.70$ at z = 15.13 and ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\approx 9.06$ at z = 9.99, which is ≈79 and ≈5.5 times lower than the observed stellar mass of CEERS-1749 and ID 1514, respectively. The stellar mass of ID 14924 is 74 times higher than the most massive simulated galaxy at z = 9.99.

The evolution of the stellar mass growth is summarized in Figure 2 and Table 1 by showing the maximum stellar mass of a subhalo in dependence of redshift for different simulation runs. The observed high-redshift galaxy candidates are, by more than 1 order of magnitude, more massive than the most massive simulated galaxies in the ΛCDM framework.

Zoom In Zoom Out Reset image size

The inferred stellar mass of observed galaxies is sensitive to the adopted SFH and initial mass function (IMF). In the following sections, we first investigate if the tension reported here of the stellar mass buildup in the early universe can be resolved if different SFHs of the observed galaxy candidates are assumed. Second, the effect of a varying IMF on the observed stellar masses is discussed.

3.1. The Minimum Inferred Galaxy Masses for Different Star Formation Histories

In order to calculate the minimum possible mass that the observed high-redshift galaxy candidates CEERS-1749, GL-z11, GL-z13, and ID 1514 can have for an invariant IMF, different sets of models of galaxies with different star formation histories (SFHs) are constructed. Using stellar population synthesis models, we let the age of the modeled galaxies vary in the range of $\left[\approx 4,400\,\mathrm{Myr}\right]$ to investigate their UV-band (1500 Å) stellar mass (including remnants)-to-light ratio, M_*/L_UV, for an invariant canonical IMF (Kroupa 2001; Kroupa et al. 2013). The lower limit is set by the implemented stellar evolution tracks of the Padova group (Marigo & Girardi 2007; Marigo et al. 2008; see also Zonoozi et al. 2019) and is roughly comparable to the mean stellar age (≈1–20 Myr) of observed high-redshift galaxies (5 ≲ z_spec ≲ 8) as found by Carnall et al. (2022).

Since the spectral energy distribution fitting analysis of high-redshift star-forming galaxies shows more consistency with increasing SFHs, here we adopt the delayed-τ model (e.g., Kroupa et al. 2020a)

$$\begin{eqnarray}&&{\psi }_{\mathrm{del}}(t)={\psi }_{0}t\exp (-t/\tau ),\end{eqnarray} \tag{ 1 }$$

and an exponentially increasing SFH

$$\begin{eqnarray}&&{\psi }_{\exp }(t)={\psi }_{0}\exp (t/\tau ),\end{eqnarray} \tag{ 2 }$$

where ψ(t) is the star formation rate (SFR), t is the age since star formation started, ψ₀ is the normalization parameter, and τ is the e-folding timescale. The mass and light of a galaxy are calculated by an integral over the SFR. Note that, using the invariant IMF, at a given t, the mass-to-light ratio is independent of the total mass that is converted into stars. This is because of cancellation of the normalization parameter, ψ₀.

The effect of the SFHs on the M_*/L_UV ratio of galaxies by adopting different values of τ is shown in the left panels of Figure 3. We set that galaxies start forming stars 200 Myr after the Big Bang, with an averaged metallicity of [Fe/H] = −2, and we assume that the mass loss from galaxies is only through stellar evolution in the form of ejected gas. Since the stellar loss due to dynamical evolution is significant only for systems with initial stellar mass less than M_* = 10⁶ M_⊙, no stars are lost by the dynamical evolution of galaxies. As can be seen, the minimum mass-to-light ratio that can be considered for these galaxies assuming different SFHs is M_*/L_UV ≈ 3.2 × 10⁻⁸ M_⊙/L_⊙.

Zoom In Zoom Out Reset image size

According to the estimated UV absolute magnitude of these objects, M_{UV,CEERS−1749} = −22.0, M_UV,GL−z13 = −20.7, M_UV,GL−z11 = −21.0, and M_{UV,ID 1514} = −22.0 mag, we obtain the lowest possible stellar masses of ${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })\,=8.53,8.01,8.13$, and 8.53 for CEERS-1749, GL-z13, GL-z11, and ID 1514 in the case of an exponential SFH, respectively, as visualized in the right panels of Figure 3.

These lower stellar mass limits just resolve the discrepancy for GL-z11, GL-z13, and ID 1514 (see the vertical dashed–dotted lines in Figure 1). In the case of CEERS-1749, the maximum stellar mass obtained in the ΛCDM simulation is ≈6.8 lower than its inferred lower limit.

3.2. Galaxy Masses for a Varying IMF

In the previous section we applied an invariant IMF but recent observations (e.g., Schneider et al. 2018; Zhang et al. 2018; Senchyna et al. 2021) suggest that the mass distribution of a stellar population may depend on its local star-forming environment. Especially metal-poor Population III stars are expected to follow a top-heavy IMF.

A theoretical framework to describe the stellar population of an entire galaxy is the integrated galactic initial mass function (IGIMF) theory, which adds up all the IMFs of star-forming regions (embedded clusters) within a galaxy (Kroupa & Weidner 2003; Weidner & Kroupa 2006). The resulting galaxy-wide IMF (gwIMF) systematically varies with the global SFR and averaged metallicity of the galaxy, i.e., the gwIMF becomes top-heavy for galaxies with SFR ≳ 1 M_⊙ yr⁻¹ and metallicities [Fe/H] < 0 (see Figure 2 of Jeřábková et al. 2018) compared to the canonical IMF. In order to study the effect of a varying IMF on the high-redshift galaxy candidates, we calculate the stellar population of galaxies assuming the latest IGIMF formalism by Yan et al. (2021) and using an IGIMF Fortran code developed by Akram Hasani Zonoozi. For simplification, we assume that all galaxies start to form stars 200 Myr after the Big Bang, with a constant SFR over time, and an average metallicity of [Fe/H] = −2. Since the gwIMF is time dependent, we require that realistic IGIMF models have to match the observed M_UV within ±1 mag in the 1σ interval of the observed redshift of the corresponding galaxy candidate. The left panel of Figure 4 shows the time evolution of the absolute UV-band magnitude for IGIMF models that fulfill these constraints for observed galaxy candidates CEERS-1749, GL-z13, GL-z11, and ID 1514. These models have constant SFRs in the range of ≈2–30 M_⊙ yr⁻¹ and thus a top-heavy IMF. The stellar masses at the observed redshifts of the galaxy candidates are shown in the right panel of Figure 4 and are systematically lower than the derived stellar masses of Adams et al. (2022), Naidu et al. (2022b), and Naidu et al. (2022a) because of the top-heavy gwIMF compared to the canonical IMF.

Zoom In Zoom Out Reset image size

While the redshifts of the galaxy candidates need to be spectroscopically verified, we use these JWST observations to quantify how quickly galaxies form in the currently most advanced cosmological simulations. Using state-of-the-art ΛCDM simulations of the IllustrisTNG and EAGLE project, we showed that the stellar mass buildup is much more efficient in the early universe than predicted by these ΛCDM models (see also, e.g., Boylan-Kolchin 2022; Lovell et al. 2022). In particular, the stellar masses of ID 1514 (Adams et al. 2022), ID 14924 (Labbe et al. 2022), GL-z11, GL-z13 (Naidu et al. 2022b), and CEERS-1749 (Naidu et al. 2022a) analyzed in Section 3 are higher by about 1 order of magnitude than the most massive galaxies formed in these simulations. In particular, massive high-redshift candidates appear more frequent at z ≳ 10 than expected in the ΛCDM framework. For example, Boylan-Kolchin (2022) argued that a volume of ≈10⁸ cMpc³ is required to explain ID 14924 with ${\mathrm{log}}_{10}(M/{M}_{\odot })=10.93$ at z = 9.92 (Labbe et al. 2022). However, the survey covers ≈10⁵ cMpc³ at z = 10 ± 1 (see Section 3 of Boylan-Kolchin 2022). The TNG100-1 and RefL0100N1504 simulations have a box volume of ≈10⁶ cMpc³ suggesting that the absence of massive galaxies in these runs is not because of a too small simulation volume.

The discrepancy between the observed and simulated stellar mass buildup could be caused by several reasons. First of all, high photometric redshifts can emerge due to dust reddening. For example, Zavala et al. (2022) demonstrated that Lyman-break galaxy candidates at z_phot ≳ 12 can resemble dusty star-forming galaxies at z ≲ 6–7. Second, it could be that the high observed stellar masses are caused by an erroneous calibration of JWST.

Furthermore, it has been argued that star formation could be much more efficient in the early universe (e.g., Naidu et al. 2022b; Harikane et al. 2022; Mason et al. 2022). Assuming that the ΛCDM model is the correct description of the universe, this would mean that the underlying galaxy formation and evolution interstellar medium models of the EAGLE and IllustrisTNG runs must be improved in order to reproduce such galaxies. For example, these simulations assume that gas above a given density threshold (e.g., Schaye et al. 2015; Nelson et al. 2019b) is able to form stars, which is likely a too simplified implementation especially for describing high-redshift galaxies. Boylan-Kolchin (2022) showed that even a 100% star formation efficiency in ΛCDM would not be enough to explain the stellar mass density measured by Labbe et al. (2022).

Another possibility is that the IMF systematically varies with the galactic properties. The IllustrisTNG and EAGLE simulations and the analysis of the Sections 3 and 3.1 assume an invariant IMF, but it is expected that metal-poor star-forming stellar populations follow a top-heavy IMF. Using the IGIMF theory we calculated the gwIMF of the observed galaxy candidates in dependence of the metallicity and SFR of forming galaxies resulting in lower stellar masses compared to an invariant canonical IMF. For this, the IGIMF theory has to be included in cosmological simulations (see, e.g., Ploeckinger et al. 2014) in order to make a firm conclusion if a top-heavy IMF can resolve the reported tension.

Finally, the present findings can also imply that structure formation is much more efficient and/or that the observed universe is even older than predicted by ΛCDM. The existence of these massive galaxies ≈300–400 Myr after the Big Bang also questions the hierarchical (bottom-up) structure formation suggesting that late-type galaxies begin to form early through the initial monolithic collapse of rotating post-Big-Bang gas clouds (Wittenburg et al. 2020) while early-type massive galaxies and associated formation of supermassive black halos form by the monolithic collapse of post-Big-Bang gas clouds with little net rotation (e.g., Kroupa et al. 2020b; Wittenburg et al. 2020; Yan et al. 2021; Eappen et al. 2022).

Evidence for an enhanced growth of structures has been reported at different astrophysical scales and redshift ranges in the observed universe. For example, Steinhardt et al. (2016) showed that the observed number density of luminous galaxies at 5 ≲ z ≲ 10 is much higher than predicted by the ΛCDM model (see their Figure 1). However, their analysis relies on the stellar-to-halo mass relation from Leauthaud et al. (2012) measured only at z = 0.2 − 1, while, e.g., Behroozi et al. (2019) suggest that there is a strong evolution at z ≳ 5. Furthermore, the existence of the massive interacting galaxy cluster El Gordo (ACT-CL J0102-4915; Marriage et al. 2011) at z = 0.87 and the Keenan–Barger–Cowie void (Keenan et al. 2013) at z ≲ 0.07 both individually falsify the hierarchical ΛCDM structure formation with more than 5σ (Haslbauer et al. 2020; Asencio et al. 2021).

An enhanced growth of structure compared to the ΛCDM paradigm is expected in Milgromian dynamics (Milgrom 1983; Angus 2009; Malekjani et al. 2009; Famaey & McGaugh 2012; Kroupa et al. 2012; Haslbauer et al. 2020; Banik & Zhao 2022). Assuming the cosmic microwave background as the z = 1100 boundary condition and because of the reduced power on <1 Mpc scales compared to ΛCDM (Angus & Diaferio 2011) due to the missing cold dark matter (CDM) component, it may be impossible to form galaxies in the early universe.

This work indicates that the currently available most advanced ΛCDM simulations cannot form galaxies as massive as observed at z_phot ≳ 10. This tension needs to be readdressed for extreme SFHs and/or if the gwIMF was top-heavy, which would reduce the stellar mass buildup through more intense feedback.

Upcoming ultradeep and wider-area JWST observations will reveal more light on the number density of such luminous high-redshift galaxies over redshift required to evaluate the significance of the here-reported tension of the stellar mass buildup of high-redshift galaxies in more detail.

We thank an anonymous referee for helpful comments that significantly improved the manuscript. This project was largely conducted at the Charles University in Prague and we acknowledge the "DAAD-Eastern-European" exchange program for financing visits at the Charles University in Prague during which we discussed several aspects of the formation and evolution of high-redshift galaxy candidates in the ΛCDM and Milgromian frameworks. We are also grateful to Nils Wittenburg, Nick Samaras, Ingo Thies, and Elena Asencio for helpful discussions on structure formation in Milgromian dynamics.