Early Structure Formation from Primordial Density Fluctuations with a Blue, Tilted Power Spectrum: High-redshift Galaxies

Recent observations by the James Webb Space Telescope (JWST) discovered unexpectedly abundant luminous galaxies at high redshift, posing possibly a severe challenge to popular galaxy formation models. We study early structure formation in a cosmological model with a blue, tilted power spectrum (BTPS) given by P(k)∝kms with m s > 1 at small length scales. We run a set of cosmological N-body simulations and derive the abundance of dark matter halos and galaxies under simplified assumptions on star formation efficiency. The enhanced small-scale power allows rapid nonlinear structure formation at z > 7, and galaxies with stellar mass exceeding 1010 M ⊙ can be formed by z = 9. Because of frequent mergers, the structure of galaxies and galaxy groups appears clumpy. The BTPS model reproduces the observed stellar mass density at z = 7–9, and thus eases the claimed tension between galaxy formation theory and recent JWST observations. The large-scale structure of the present-day Universe is largely unaffected by the modification of the small-scale power spectrum. We conduct a systematic study by varying the slope of the small-scale power spectrum to derive constraints on the BTPS model from a set of observations of high-redshift galaxies.


INTRODUCTION
The so-called Λ Cold Dark Matter (ΛCDM) cosmology successfully reproduces a broad range of observations of the large-scale structure of the Universe, and thus has been established as the standard cosmological model.A vital element of the standard model is the primordial density fluctuations generated in the very early universe with a nearly scale-independent power spectrum.While the large-scale density fluctuations have been observationally probed to k ∼ 1 Mpc −1 , the amplitude and the shape of the power spectrum on smaller, (sub-)galactic length scales are poorly constrained (e.g., Hlozek et al. 2012;Bullock & Boylan-Kolchin 2017).Hence, theoretical studies on galaxy formation often rely on significant Corresponding author: Shingo Hirano shingo-hirano@kanagawa-u.ac.jp extrapolation of the assumed scale-invariant primordial power spectrum (PPS).
Various possibilities have been proposed from the physics of the early universe that posit deviations from scale-invariance.Blue, tilted, or enhanced power spectra can arise in beyond-standard cosmological models (e.g., Covi & Lyth 1999;Martin & Brandenberger 2001;Gong & Sasaki 2011;Inman & Kohri 2023), and yield several interesting cosmological and astrophysical consequences (e.g., Clesse & García-Bellido 2015;Germani & Prokopec 2017).Earlier in Hirano et al. (2015), we studied early structure formation in models with a blue, tilted power spectrum (BTPS).It was shown that the enhanced small-scale density fluctuations drive the formation of nonlinear structure and Population III stars at very early epochs.
In this paper, we study the formation and abundance of the first galaxies in the BTPS model in light of re-cent observations by the James Webb Space Telescope (JWST).Some galaxies (candidates) with unexpectedly high stellar masses have been discovered (e.g., Finkelstein et al. 2022;Labbé et al. 2023).Boylan-Kolchin (2023) concludes that the inferred high stellar masses of the observed galaxy candidates require an extremely high star formation efficiency far exceeding the plausible values of ϵ ≲ 0.3 suggested by popular galaxy formation models (Gribel et al. 2017;Tacchella et al. 2018;Behroozi et al. 2020).The challenge brought by recent JWST observations motivates us to reconsider the detailed physics of galaxy formation in the early universe or the standard cosmology model.
Modification of PPS may provide a viable solution by promoting early structure formation (Parashari & Laha 2023;Padmanabhan & Loeb 2023).Parashari & Laha (2023) compute the cumulative comoving stellar mass density (CCSMD) by adopting a modified form of PPS similar to Hirano et al. (2015).Their model can successfully reproduce the observed CCSMD without assuming an unrealistically high star formation efficiency.They further argue that observations of high-redshift galaxies can provide invaluable information on small-scale density fluctuations that cannot be directly probed.It is crucial and timely to study the formation of highredshift galaxies by performing cosmological simulations for the BTPS model.By calculating the halo mass functions for a wide range of models with varying the initial power spectrum, we are able to derive constraints on the PPS from the observed CCSMD.We also compute the ultra-violet luminosity function (UVLF) of galaxies in the BTPS model and compare them with the Hubble Space Telescope observations.
Throughout the present paper, we adopt the cosmological parameters with total matter density Ω m = 0.3153, baryon density Ω b = 0.0493 in units of the critical density, a Hubble constant H 0 = 67.36km s −1 Mpc −1 , the root-mean-square matter fluctuation averaged over a sphere of radius 8 h −1 Mpc σ 8 = 0.8111, and primordial index n s = 0.9649 (Planck Collaboration et al. 2020).

NUMERICAL SIMULATIONS
We largely follow the method of our previous study (Hirano et al. 2015) to perform cosmological N -body simulations for three sets of models with different PPS.
whereas the one with enhancement at small scales is given by We consider models with different PPS to investigate the dependence of CCSMD and UVLF on the parameters (k p and m s ): (1) We first fix the pivot scale at k p = 1 h cMpc −1 (h comoving Mpc −1 ) and adopt two values m s = 1.5 and 2.0, which we call "soft-tilt" and "hardtilt", respectively.We note that our soft-tilt model is in marginally acceptable parameter space from available observational constraints (see Figure 2 in Parashari & Laha 2023).
(2) Next, we examine the influence of the pivot scale by adopting k p = 0.3 and 5.0, which we call shallow model and deep model, respectively.The latter is motivated by the recent study on the structure of Milky Way satellite galaxies (see Figure 8 in Esteban et al. 2023).
For these BTPS models, we adjust the normalization σ 8 to ensure that the amplitude of the fluctuations at wavelength larger than the pivot scale is the same as in the standard ΛCDM model ("standard model").Figure 1 shows the resulting matter power spectra at z = 1089.

Cosmological simulations
We use the public code MUSIC (Hahn & Abel 2011) to generate cosmological initial conditions.We employ 512 3 dark matter particles in comoving cubes of L box = 10, 25, and 50 cMpc h −1 .Three different simulation volumes are adopted to accurately derive the statistics of nonlinear structure over a wide range of length and mass scales.The particle mass is 6.52 × 10 5 , 1.02 × 10 7 , and 8.15 × 10 7 M ⊙ , respectively.Dark matter halos with M halo = 10 8 M ⊙ are resolved by more than 150 particles in our highest-resolution runs.Table 1 summarizes the basic simulation parameters.We use the parallel N -body code GADGET-2 (Springel 2005) to follow structure formation from redshift z = 99 to 0.

Cumulatice comoving stellar mass density
We follow Boylan-Kolchin (2023) and Parashari & Laha (2023) to compute the cumulative comoving stellar mass density (CCSMD) from the simulation outputs.We determine the required star formation efficiency (ϵ) to reconcile the recent JWST observation.
We first run a Friends-of-Friends group finder with linking length b = 0.164(m/ρ(z)) −3 , where m is the Nbody particle mass and ρ(z) is the mean density of the universe at redshift z, to obtain the halo mass function dn(M, z)/dM .We combine the halo mass functions from our simulations with different volumes (particle masses) to construct the halo mass function for each model in a wide mass range.Next, we calculate the cumulative comoving mass density of halos Finally, we compute the CCSMD with stellar mass larger than 156, and also assume that the star formation efficiency ϵ (≤ 1) is constant over the redshift range we consider.

STELLAR MASS DENSITY
We focus on the cumulative comoving stellar density (CCSMD) as a primary statistic to compare our simulation results and recent JWST observations.We first discuss nonlinear structure formation in the BTPS models with different tilt, m s .Then we study the dependence of CCSMD on the other parameter k p .We show that our "shallow" models with k p < 1.0 h cMpc −1 reproduce the observed CCSMD.We identify regions in the k p −m s parameter space where the observed CCSMD is reconciled with a reasonable choice of star formation efficiency of ϵ = 0.1 − 0.3.

Dependence on the tilt, m s
Figure 2 shows the projected density distributions for the models with different m s .The enhanced small-scale density fluctuations yield stronger density contrast in the BTPS models than in the standard model; numerous nonlinear objects (halos) are already formed by z = 9.We compare the halo mass functions in Figure 3, which will be used later to discuss the cumulative stellar mass.The number density of halos with mass 10 9 < M/M ⊙ <   10 11 differs by an order of magnitude between the hardtilt model and the standard model.
The middle panels of Figure 2 show that massive galaxies (halos) are strongly clustered in the BTPS models.The relative galaxy bias to the underlying dark matter distribution will provide crucial information on the nature of the early galaxy population.Future widefield observations of galaxy distribution and clustering by JWST or the Roman Space Telescope will enable us to discriminate theoretical models of galaxy formation (Muñoz et al. 2023).
The bottom panels of Figure 2  Figure 4 shows the CCSMDs for the three models at z = 9 and 7.5.There we assume ϵ = 0.1 (solid lines) and 0.3 (dashed lines) with the latter corresponding to the plausible upper limit suggested by recent galaxy formation models (Gribel et al. 2017;Tacchella et al. 2018;Behroozi et al. 2020).The green regions indicate the JWST observations adopted from Parashari & Laha (2023), who calculated the CCSMD from the observation of Labbé et al. (2023) with the corresponding spectroscopic updates in two redshift bins z ∈ [7, 8.5] and [8.5, 10] using the three most massive galaxies.
As has been suggested in the recent literature, a large value of ϵ ∼ 0.3 is required for the standard model to reproduce the mean value of the observed CCSMD at z = 9 (Figure 4a).Even with the "maximal" ϵ, the CC-SMD does not reach the lower limit of the uncertainty range of the observations at z = 7.5 (Figure 4b).The CCSMD of the soft-tilt model is about ten and three times larger than the standard model at z = 9 and 7.5, respectively.Thus the soft-tilt model requires a moderate star formation efficiency of ϵ ∼ 0.1−0.3 to match the observed CCSMDs at z = 9 and 7.5.Table 1 summarizes the star formation efficiency required for the CCSMD of each model to exceed the observationally inferred lower limit (ϵ min ), mean (ϵ mean ), and upper limit (ϵ max ).

Dependence on the pivot scale, k p
We have set the fiducial pivot scale at k p = 1 somewhat arbitrarily.To examine the effect of the pivot scale on the CCSMD, we consider "shallow" and "deep" models with k p = 0.3 and 5, respectively.Figure 5 shows that the resulting CCSMDs are nearly the same as that of the standard model at the massive end, and the values are slightly larger at small masses.Thus the deep model requires a very high star formation efficiency to reconcile with the JWST observations, similar to the standard model.
Interestingly, our shallow model successfully increases the CCSMD to match the observation, thus relaxing the star formation efficiency requirement.Even though the enhancement on the small-scale power appears small compared to the other models (soft, hard, deep), the shallow model effectively increases the abundance of galactic dark halos at the relevant mass and redshift.Note that a small bump in the power spectrum around k ∼ 1 could also explain the apparently rapid galaxy formation in the early universe (Padmanabhan & Loeb 2023).

Parameter survey for low-k p
Motivated by the success of our shallow model, we explore a wider range of models with low k p and derive an observationally allowed parameter range of the power spectrum.Figure 6 summarizes the results on the k p -m s plane.We accept a model if the CCSMD value exceeds the central value of the JWST observations (shaded regions in Figures 4 and 5).We judge the models at two different redshifts (z = 7.5 and 9) and with two different star formation efficiencies (ϵ = 0.1 and 0.3) independently.In Figure 6, models on and above the respective boundary lines are "acceptable".For instance, models with m s = 1.1 and with k p < 0.2 successfully reproduce the observed CCSMD and hence "accepted".
We argue that a slight enhancement of the PPS in the previously unexplored range of k p = 0.1 − 0.2 could explain the rapid formation of galaxies in the early Universe.Among low-k p BTPS models, we select one model We calculate the CCSMDs for each model at z = 9 (blue) and 7.5 (red) with ϵ = 0.1 (solid) and 0.3 (dashed).We divide the "accepted" and constrained regions by the dashed and solid lines.Above the line for each model, the CCSMD matches the mean of the observed range (green regions in Figure 4).
with k p = 0.2 h cMpc −1 and m s = 1.3 as the "minimal shallow" model that could explain available observations with the modest change of the PPS.

LUMINOSITY FUNCTION
The ultra-violet luminosity function (UVLF) may be another useful probe of the small-scale power, although the resulting constraints depend much on the adopted models of galaxy formation and mass-to-light conversion (Muñoz et al. 2023).
We follow Sabti et al. (2022a,b) to compute the UVLF from the simulation outputs.To simplify the model comparison, we adopt the best-fit parameter values in their fiducial model obtained from the Hubble Space Telescope data.We first calculate the halo mass accretion rate at redshift z by comparing halo masses of the same halo identified by the Friends-of-Friends algorithm at different redshift with dz = 0.1 as Ṁh (z) = (M h (z) − M h (z + dz))/dz.Next, we translate the halo mass accretion rate to the star formation rate (SFR) by assuming a double power-law relation, with the following forms, The best-fit values of six free parameters are {α * β * , ϵ s * , ϵ i * M s c , M i c } = {-0.594,0.611, -1.96, -2.17, 2.95, 12.1}.We relate the UV luminosity of a galaxy to the SFR as (Kennicutt 1998;Madau et al. 1998), with a popular choice of κ UV = 1.15 × 10 −28 M ⊙ s erg −1 yr −1 (Madau & Dickinson 2014).We then convert the UV luminosity to the absolute magnitude (Oke & Gunn 1983), (Kennicutt 1998;Madau et al. 1998), Using the above formulae, we calculate the absolute UV magnitudes (M UV ) from the halo mass accretion rates ( Ṁh ) obtained from the simulations.We calculate the UVLF for each model at different redshifts.
Figure 7 compares the recent observation (Bouwens et al. 2021) to the calculated UVLFs for the minimal shallow (k p = 0.2 h cMpc −1 and m s = 1.3) and standard models.Compared to the standard model, the minimal shallow model increases the UVLF, but does not show significant deviations from the observed range.One reason is that the BTPS model promotes formation of many low-mass galaxies but does not affect the formation of massive galaxies, which contribute to the UVLF in the relevant luminosity range.Interestingly, the minimal shallow model is closer to the observational data than the standard one with the model employed here.To the extent of this comparison, the current UVLF observations do not constrain the shallow tilt BTPS model.
We note that the star formation efficiencies adopted in the above UVLF model and our CCSMD model may not necessarily be consistent.The CCSMDs are calculated by assuming a fraction ϵ of the gas mass is converted into stars, whereas the UVLF model assumes a fraction f * of the accreted halo mass is converted to stars at a given instant.It is not straightforward to compare the two efficiencies directly, for one is related to the resulting stellar mass by a given epoch and the other is the instantaneous formation rate.Also, the UVLF model we consider here involves an effective mass-to-light ratio (Equation 10) which is uncertain and depends critically on the assumed stellar initial mass function.It is beyond the scope of the present paper to develop a consistent galaxy formation model, but it would be interesting to explore, for example, semi-analytic models of galaxy formation that treat the star formation process in a self-consistent manner.We also note the effective halo mass range probed by the JWST CCSMD is different from that of the HST UVLF function, the latter of which probes overall larger halos with M halo = 10 11 − 10 12 M ⊙ and hence the effective power at lower wavenumber (Muñoz et al. 2023).

DISCUSSION
The result of our cosmological simulations is largely consistent with the estimate of Parashari & Laha (2023) based on an analytic halo mass function.Our simulations confirm that the standard ΛCDM model requires unrealistically high star formation efficiencies of ϵ > 0.3 to reconcile the observed CCSMD.Although there have already been proposals for galaxy formation physics that can realize a high star formation efficiency (e.g., Dekel et al. 2023), a slight modification of the PPS may be another promising solution that alleviates the need for significant deviation from the currently popular galaxy formation models.
Interestingly, recent JWST observations also show that there are many galaxies with clumpy structures and galaxies in the process of mergers in proto-cluster environments at z = 7 − 9 (Hashimoto et al. 2023;Hainline et al. 2023).The nonlinear structure forms early in the BTPS model and assembles rapidly via mergers.Thus, galaxies at z = 7 − 9 tend to appear clumpy, as seen in Figure 2.
The PPS strongly affects the formation epoch and the properties of the Population III stars.Hirano et al. (2015) run a set of hydrodynamics simulations starting from nearly the same BTPS and find that very massive stars with mass exceeding 100 M ⊙ are formed early.Recently, Xing et al. (2023) discovered a metal-poor star that shows definite chemical signatures of the so-called pair-instability supernova caused by a very massive star early in the formation history of the Milky Way.Formation efficiency and epoch of very massive stars may provide another hint on the small-scale density fluctuations.
The existence of super-massive black holes in the early universe generally suggests rapid growth of structure and formation of appropriate seed black holes at early epochs (e.g., Inayoshi et al. 2020).Larson et al. (2023) discovered a massive black hole candidate at z = 8.6.Stellar-mass black holes are formed as remnants of massive Population III stars as early as z ∼ 50 − 100 in the BTPS model, leaving enough time for the seeds to grow by mass accretion to massive black holes by z = 8.Frequent halo mergers realized in the BTPS model may also enhance the formation of massive black holes in the early universe (Wise et al. 2019;Regan 2023).
Tight constraints can be placed on the slope or the amplitude of PPS at sub-galactic length scales from observations of the history of cosmic reionization.Contribution to early reionization from individual Population III stars is severely limited by the measurement of the Thomson optical depth by Planck (Visbal et al. 2015).Minoda et al. (2023) derive constraints on the "running" of the PPS from CMB and recent observations of high-redshift galaxies, which still allow both negative (red) and positive (blue) running to a modest degree.Enhanced small-scale power promotes early star formation and also increases inter-galactic medium clumping (see Figure 2).The photon escape fraction of galaxies, an important quantity in the study of cosmic reionization, is uncertain and likely affected by the gas clumping and the abundance of substructure within early galaxies (Yeh et al. 2023).Also, radiative and mechanical feedback effects from Population III stars may actually suppress, at least temporarily, star formation in small galaxies that form later around z ∼ 10 − 15 (Bromm & Yoshida 2011).Overall, there still remains substantial uncertainty in both observations of reionization history (Hinshaw et al. 2013;Planck Collaboration et al. 2020;Forconi et al. 2023) and in astrophysical modelling (Fi-alkov 2022).It would be necessary to perform detailed numerical simulations of galaxy formation and reionization to place constraints on the BTPS model.
Structure of galactic halos may provide another clue on the small-scale power.Since low-mass satellite galaxies are too faint to be observed in the distant Universe, observations of the nearby Universe are often used to infer the small-scale density fluctuations (e.g., Esteban et al. 2023).Interestingly, Homma et al. (2023) recently reported that the abundance of ultrafaint dwarf galaxies in the Milky Way is higher than the theoretical prediction from the standard ΛCDM model with a single power-law PPS, posing a "too many satellites problem".Detailed galaxy formation models combined with a broad class of observations of sub-galactic objects are necessary to understand the formation and evolution of small-scale structures and to place tight constraints on the PPS.Tkachev et al. (2024) consider a modified matter power spectrum with an approximately Gaussian bump at sub-Mpc scales and study the effect on the halo mass function.They argue that observations of high-redshift galaxies are crucial to infer the feature of the small-scale power spectrum.
JWST has opened a new window into the distant universe.Future observations by JWST and other widefield cosmology surveys will reveal the physical properties and the formation history of galaxy populations at high redshift.Upcoming data may ultimately provide new insight into the early universe physics that generates the primordial density fluctuations from which the first galaxies are formed.
Numerical computations were performed on Cray XC50 at CfCA in National Astronomical Observatory of Japan and Yukawa-21 at YITP in Kyoto University.Numerical analyses were in part carried out on the analysis servers at CfCA in National Astronomical Observatory of Japan.This work was supported by JSPS KAKENHI Grant Numbers JP21K13960, JP21H01123, and JP22H01259 (S.H.), and MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Structure and Evolution of the Universe Unraveled by Fusion of Simulation and AI; Grant Number JPMXP1020230406, Project ID hp230204) (S.H. and N.Y).

Figure 2 .
Figure 2. Projected density distributions δ + 1 = ρ/ρ and DM halos with M halo ≥ 10 9 M⊙.The top and middle panels show the structure at z = 9 in a volume of side-length 10 cMpc h −1 , whereas the bottom panels show large-scale structure at z = 0 with a side-length of 50 cMpc h −1 .The left, center, and right panels are for standard, soft-tilt, and hard-tilt models.The circle size in the middle panels scales with the halo mass.The inset in panel (c) is a zoom-in image of one of the most massive halos.

Figure 3 .
Figure 3. Halo mass functions at z = 9, 7.5, and 0 from left to right.The solid lines show our simulation results and the dashed lines represent the analytical Sheth-Tormen mass functions (Sheth & Tormen 1999) calculated by the public code genmf (Reed et al. 2007).The gray, light blue, and dark blue lines are for the standard, soft-tilt, and hard-tilt models.We construct the simulation data by combining results for different box sizes, L box = 10, 25, and 50 cMpc h −1 .
suggest that the largescale structure of several tens Mpc in the present-day universe (z = 0) remains essentially the same.The enhancement of PPS at the small length scales does not ruin the success of the standard model.The same holds for the halo mass function at z = 0 plotted in Figure 3(c).

Figure 5 .
Figure 5. Same as Figure 4 but for the standard (gray), shallow (light red), and deep (red) models at z = 9 (panel a) and 7.5 (b).

Figure 6 .
Figure6.Summary of our simulation-based constraints in kp-ms space.We calculate the CCSMDs for each model at z = 9 (blue) and 7.5 (red) with ϵ = 0.1 (solid) and 0.3 (dashed).We divide the "accepted" and constrained regions by the dashed and solid lines.Above the line for each model, the CCSMD matches the mean of the observed range (green regions in Figure4).

Table 1 .
List of models with different primordial power spectrum