Excess of high-z galaxies as a test for bumpy power spectrum of density perturbations

Modiﬁed matter power spectra with approximately Gaussian bump on sub-Mpc scales can be a result of a complex inﬂation. We consider ﬁve spectra with diﬀerent amplitudes A and locations k 0 and run N-body simulations in a cube (5 Mpc/h ) 3 at z > 8 to reveal the halo mass functions and their evolution with redshift. We have found that the ST formula provides a good approximation to a such kind of matter spectra. In the considered models the dark matter halo formation starts much more earlier than in ΛCDM, which in turn can result in an earlier star formation and a nuclear activity in galaxies. At z = 0 the halo mass functions are hardly distinguishable from the standard ΛCDM, therefore the models with the bumpy spectra can be identiﬁed in observations by their excess in number of bright sources at high redshift only.


INTRODUCTION
It is commonly believed that an inflation predicts a spectrum of density perturbations similar to a power law or which can be approximated by a power law over the wide range of scales up to the current horizon value.Indeed, a smooth spectrum with barely discernible deviation is a common feature of many inflationary models based on the single scalar field theory.For example, the density perturbation spectrum in chaotic inflation with the single massive scalar field looks like in terms of q-scalar (see Lukash (1980), Lukash (2010)), where k is a wave number, k1 is the wave number at the end of inflation, and m is a mass of inflaton.The running parameter of this spectrum cannot be detected at the available accuracy measurements.
During the inflationary stage of the Universe evolution, processes that disrupt the smooth spectrum of density perturbations could have taken place.Many such models have been proposed in the context of inflation with two fields, ⋆ mtkachev@asc.rssi.ru† spilipenko@asc.rssi.ru‡ helen@asc.rssi.ru§ lukash@asc.rssi.ruresonant amplification with oscillatory features, and others (see Inomata et al. (2023), numerous references wherein, and Kawasaki et al. (2006); Saito et al. (2008); Ballesteros & Taoso (2018); Braglia et al. (2020); Iacconi et al. (2022); Cai et al. (2020); Zhou et al. (2020); Mishra & Sahni (2020)).The motivation for this research has been greatly enhanced by recent observational arguments in favour of primordial black holes as viable dark matter candidates.These arguments are the detection of gravitational waves from merging black holes (see Abbott et al. (2016)) and no detection of dark matter particle candidates of the Standard Model of particle physics (see Aleksandrov et al. (2021), Nevzorov (2023)).To generate primordial black holes it is necessary to amplify the density perturbation spectrum at some scales, which determines the mass spectrum of black holes.Possible mechanisms of the amplification and matter power spectra can be found in aforementioned references.Unless one assumes the fine-tuning of related inflationary quantities, it follows that a peak-like features (also called a bump) can appear at other scales too (see, e.g.Fig. 2 of Ivanov et al. (1994)).If a spectral feature appears on the scales k < 10 3 h/Mpc, it will affect nonlinear large structure formation, in particular, abundance of galaxies, their mass growth histories and clustering.The exact shape of the spectral feature depends on the inflation model, dozens of examples can be found in the references listed above.Instead of analysing nonlinear structure formation for every particular variant of such a peak, we assume that a single general shape can reproduce the most significant consequences of the introduction of the peak.We choose a Gaussian as such a general shape.In the Appendix A we demonstrate how our results could change if one uses the same method to analyse a particular example of a spectral feature predicted by one of the inflation models.It is important to remark, that the location(s) of the break(s) in the inflaton potential and, consequently, the peak(s) in the matter spectrum cannot be selected a priori (they are free parameters of any such a theory).For example, they may arise at sub-Mpc scales and be responsible for some unusual feature in halo mass function at low and high redshifts.This case is under consideration in this paper.
The power spectrum of density perturbations is now measured at comoving wavenumbers k < 1 h/Mpc from the Cosmic Microwave Background (CMB), large scale structure of the Universe, Ly-α forest (Chabanier et al. 2019).No deviations from the power-law primordial spectrum have been found.The deviations at smaller scales, 1 < k < 10 3 h/Mpc, would affect galaxy formation.A simple and straightforward estimate using Press-Schechter (PS) formalism (Press & Schechter 1974) shows that addition of a peak-like feature at k > 1 h/Mpc will result in amplification of the number density of galaxy or dwarf-sized haloes at high redshifts, while this amplification will be almost smoothed away at redshift z = 01 .This means that it is difficult to find a trace of the bump in the power spectrum at z = 0, and observations at the Epoch of Reionization (EoR) are more suitable for that.
Recently surprisingly high number of massive galaxies at z ≥ 9 has been discovered in James Webb Space Telescope (JWST) data (Naidu et al. 2022b;Castellano et al. 2022;Finkelstein et al. 2022;Donnan et al. 2023;Labbé et al. 2023).It is actively debated now, whether these findings are in tension with the abundance of haloes predicted by the ΛCDM (Boylan-Kolchin 2023; Lovell et al. 2022;Chen et al. 2023;Prada et al. 2023;Shen et al. 2023).A number of modifications of the cosmological model have been proposed, including modifications of the power spectrum with a tilt at blue end (Parashari & Laha 2023;Hütsi et al. 2023) and a bump (Padmanabhan & Loeb 2023).Our aim is to present the expected deviations from the ΛCDM halo abundance for cosmological models with the bump in the power spectrum.We demonstrate the use of bumpy power spectrum models by applying the Extreme Value Statistics (EVS) (Gumbel 1958) to compare observations with theoretical predictions.Recently Lovell et al. (2022) have shown using EVS that JWST observations can be explained by ΛCDM mass function only if one adopts very efficient star formation at the EoR: about 100% of baryons in haloes must transform into stars.
The halo mass function, which is needed to construct the EVS, can be obtained with modifications of the PS formalism, e.g. the more precise Sheth-Tormen (ST) formula (Sheth & Tormen 1999) but there is a couple of possible caveats here.As shown by Klypin et al. (2011), ST approximation gives a fairly accurate fit for z = 0 for a wide range of halo masses, however, it tends to overpredict the abundance of haloes at higher redshifts.This downside can be at least partially circumvented by adjusting the radius R of the real space spherical top-hat window function Ŵ (kR) in the ST formula that we used (for more accurate improvements over the ST approximation see e.g.Tinker et al. (2008), Behroozi et al. (2013) and Despali et al. (2016)).Furthermore, while ST is known to reproduce the ΛCDM mass function quite well (Jenkins et al. 2001;Reed et al. 2003), it is not so good for cosmological models with more complicated spectra, e.g.Warm Dark Matter which has a sharp cut-off at small scales (Angulo et al. 2013).So, before using ST or any other theoretical or empirical mass function for calculating the number density of haloes in cosmological models with bumpy spectra, we need to check that this technique is applicable.For this, we run a set of N-body cosmological simulations with different bumpy spectra, as well as a standard ΛCDM spectrum, and compare the simulated mass functions with the theoretical one.
Some investigations of halo formation in cosmological models with bumpy power spectra have been done in the literature.Authors of Knebe et al. (2001) have analysed the impact of a positive or negative bump on the galaxy cluster mass function.They have found that PS formalism predicts the mass function quite well, but they have simulated different part of the power spectrum, k < 1 h/Mpc.Since the slope of the power spectrum transfer function changes with k and become more steep at higher k, it is not clear whether the results of Knebe et al. (2001) are applicable for the range of k we are interested in.Later Bagla & Prasad (2009) have simulated a model with the bump at k ≈ 1 h/Mpc, but the simulation resolution was not very high and the change of the mass function they found is moderate in comparison with the noise.
The paper is organised as follows: in Section 2 we describe the considered models; in Section 3 we describe the details of our numerical simulations and provide the list of parameters used to run the N-body simulations.In Section 4 we present the ST mass functions adopted for modified matter spectra.In Section 5 we discuss the early halo formation realising in considered cosmological models.Finally, we leave the Section 6 for discussion of various implications.

COSMOLOGICAL MODEL
To study the influence of the peak-like features of density perturbations spectrum on the evolution of dark matter haloes we focused on the power spectra, which can be defined as a product of the standard ΛCDM spectrum and a Gaussian bump of the following form: where k is a wave number, A, k0, and σ k are bump parameters.We assume an equal value of σ k = 0.1 in all models.Numerical values of the parameters can be found in the Table 1, whereas the Fig. 1 illustrates the matter spectra.One can note that the shape (1) is not predicted directly by simple modifications of the inflation model (see, e.g., Inomata et al. (2023)).We consider the Gaussian shape as a universal approximation of the peak shape in various models.
The width of the peak, σ k , does not play an important On both panels grey dashed lines represent the Nyquist frequency for (5 M pc/h)3 cube and 512 3 particles resolution.The parameters of each individual spectrum can be found in the Table 1.
role in our studies: it is easy to show using theoretical halo mass functions (e.g., PS or ST) that the mass function does not change when Aσ k = const for σ k ≪ 1.So there is enough to vary only one parameter, in our case the amplitude A. We do not consider bumps with σ k ≳ 1 as this would made the study more complicated.
As one can see from the Table 1, all modified spectra from gauss 1 to gauss 4 have Gaussian bump maxima at different values of k0.The last spectrum, gauss 5, is the exception, and has the same k0 as gauss 4, but smaller amplitude A. This was done to study the case, where the amplitude of the Gaussian bump has small impact on the resulting halo mass function.

N-BODY SIMULATIONS
We have run a series of 6 dark matter only simulations in a box size of (5 M pc/h) 3 , 512 3 particles each.5 of simulations correspond to different modified matter power spectra, and one corresponds to the standard ΛCDM one.
For our simulation suite we use the publicly available Nbody code GADGET-22 (Springel 2005), which is widely used for cosmological simulations.This code uses a combined Tree + Particle Mesh (TreeMP) algorithm to estimate the gravitational accelerations for each particle by decomposing the gravitational forces into a long range term, computed from Particle-Mesh methods, and short scale interactions from the nearest neighbours using Tree methods, and can be used with periodic boundary conditions in the comoving frame.The code is designed to be MPI parallel, enabling it to efficiently utilise distributed computing resources by dividing the computational tasks among multiple processors, resulting in faster execution and scalability (i.e.O(N log N )), so it can handle a large number of particles with reasonable computational resources.
All simulations with modified matter spectra start at z = 1000 in order to account for potential early formation of virialised structures, while the simulation with ΛCDM spectrum starts at z = 300.The final redshift of the simulations was set to z = 8, in order to reduce possible artefacts due to space periodicity of initial conditions in the relatively small-sized box.
We generate initial conditions for our simulations using publicly available code ginnungagap 3 , the matter power spectrum is defined for each simulation individually by applying the transform function (1), and where ΛCDM power spectrum is generated with publicly available code CLASS (Blas et al. 2011).We use the same initial random seed number for all our simulations, so they differ only by the amplitude of the power spectrum.When the initial conditions were generated, the amplitude of the longest wavelength mode is within 20% of the theoretical one, so the cosmic variance does not significantly affects the high-mass end of our mass functions.
A total of 100 snapshots for each simulation were stored at redshift intervals equally-spaced in logarithmic scale, starting from z = 35 to z = 8.Halo analysis was performed with publicly available code AHF4 (Knollmann & Knebe 2009) with the assumption that each halo consists of no less than 50 particles, while virial overdensity criterion was assumed 200 ρcrit.

HALO MASS FUNCTION
For each simulation we plotted a halo mass function (HMF) at two different redshifts: z ≃ 13 and z ≃ 8, where the former was chosen as a typical high-redshift value for the most distant galaxies found today.While z ≃ 8 is the final redshift, up to which our simulation was calculated.
The comparison of simulated mass functions with the theoretical (ST) ones is shown in Fig. 2. Each panel (from left to right and from top to bottom) includes HMFs for ΛCDM (black lines) spectrum and one of the modified spectra (see legend for the specific colour).Solid lines represent HMF for z ≃ 8, while dashed lines do HMF with z ≃ 13.As expected, an increase in the power spectrum amplitude at certain wave numbers results in a corresponding increase in the abundance of haloes with masses associated with those wave numbers, as demonstrated by the Fig. 2. In addition, each panel includes the corresponding HMFs, obtained from ST approximation5 .One can see from Fig. 2 that HMFs obtained from simulations with modified power spectrum appear to be in good accordance with the HMFs, obtained via ST approximation for the corresponding spectra.We compute ST HMF using the real space spherical top-hat window function.We have checked that it gives better coincidence with simulated HMFs than using Fourierspace top-hat or Gaussian window function for the bumpy power spectra.
Considering that our simulations final redshift is z = 8, in Fig. 3 we also compare the theoretical (ST) HMFs with the standard ΛCDM spectrum and the modified spectrum gauss 1 at redshifts z = 9, 6, and 0. As one can see, in case of z = 6 and z = 0 lines, as redshift decreases, the difference between HMFs for ΛCDM spectrum and HMFs for modified spectrum diminishes.This trend suggests that the impact of the modification on the HMF becomes less significant as the Universe evolves to lower redshifts.
Additionally, we had some concerns regarding the statistical adequacy of our simulation suite due to a relatively small size of simulation boxes (5 M pc/h) 3 .Thus, we considered several snapshots from the simulation Extremely Small MultiDark Planck (ESMD) carried out by an international consortium within the framework of the MultiDark project6 , and available at CosmoSim database7 .ESMD has a much larger box size, (64M pc/h) 3 .As one can see, HMFs derived from the ESMD simulation, represented by the dashed black line on the Fig. 3, closely aligns with the theoretical ST HMFs for the ΛCDM spectrum, which indicates that ST HMFs that we plotted does not noticeably overestimate the number of haloes in the case of our simulation suite as well.
Overall, the close agreement between the ST HMFs and the simulated data in case of both ΛCDM spectrum and the modified spectra validates the utility of the ST approximation in capturing the halo mass distribution in cosmological simulations with modified spectra.

EXTREME VALUE STATISTICS WITH BUMPY POWER SPECTRA
Within this Section, we discuss how the effects of modified cosmological power spectrum might affect the observations of halo formation in the early Universe.We test how well different power spectra fit the observational data on z > 8 galaxies with EVS approach.The idea behind the EVS is to test if the most extreme observed value of some property is indeed an outlier.In our case, we test the maximal halo mass of galaxies observed by JWST at high redshifts.To infer halo masses from observations, we use the stellar fraction f * and we find for each observed galaxy the minimal value of f * for which this galaxy is not a 3σ outlier from the theoretical maximal halo mass distribution.
The data we use are presented in the Table 2 and also displayed on the Fig. 4. Filled cyan dots represent the last two sources from Table 2, while filled pink dots represent the rest of the sources from the Table 2.When comparing observations to the theory, one needs to take into account the Eddington Bias (Eddington 1913).Due to the steepness of the HMF and galaxy stellar mass function, measurements of low mass objects can be easily upscattered.White-space dots with coloured border represent masses, corrected for Eddington Bias as follows: where ϵ is the local slope of the underlying halo mass function, and σ ln M is the uncertainty in the halo/stellar mass estimate.
For the given area of the survey we compute the PDF of the most massive halo as a function of z for N observations (see Lovell et al. (2022)) as: where f (m) and F (m) for a fixed fraction of the sky f sky between redshifts zmin and zmax are calculated as follows: where ntot is a normalisation factor giving the total number density of haloes: Using the aforementioned equations along with the Eq. ( 3), while considering an observational survey volume with an Left panel presents an approximation of the masses of the most massive haloes using the theoretical (ST) mass function, depicted with thin solid lines.
Right panel depicts the masses of the most massive halo in our numerical simulations at a given redshift, represented by the thick solid lines.
The observational data from Table 2 are shown as filled cyan and pink dots with error bars.The filled light-grey and light-red areas represent the 3σ ranges for ΛCDM and gauss 1 mass approximations derived from EVS, and they are depicted in both the left and right panels.
Bottom panel represents the difference between the simulation and ST masses respectively.
the other hand, the observational data at pre-reionization epoch (z > 10) is missing, while some theoretical models predict high star formation efficiencies, up to f * = 0.8 (Susa & Umemura 2004).
To elucidate the relationship between the bump and the value of f * , we calculate the values of this quality in ΛCDM and gauss 1 for each of high-z sources.The Table 2 shows the minimal values of f * for ΛCDM and gauss 1 respectively, which are required for the given galaxy to fall into 3σ range.As can be seen, with the assumption of the bumpy power spectrum, most of the galaxies need several times smaller value of f * .
As can be seen from Fig. 4, the masses of the most massive haloes in the simulations with modified spectra differ significantly from those in the ΛCDM model above some redshift, which depends on the position of the bump.At smaller redshifts this quantity converges to the ΛCDM model.In particular, our model gauss 1 fits the observations much better than the ΛCDM, allowing smaller, more reasonable values of the stellar fraction in high-z galaxies.
Both masses approximated from ST and the halo masses from simulations with modified spectra exhibit a similar trend with the deviations less than 30% for the majority of spectra up to z ≈ 12, although some discrepancy might be observed at higher redshift values (see the bottom panel of Fig 4).This implies that ST approximation can be used for calculating the EVS for bumpy power spectra.

DISCUSSION AND CONCLUSIONS
In the paper we studied a cosmological model with modified power spectrum of density perturbations.In this model a matter power spectrum has a Gaussian bump with a certain amplitude A and a scale location k0 (see eq.( 1)).We considered five different models with parameter values summarised in the Table 1.All models were studied with both numerical N-body simulations and analytical ST approximation.The latter supported the results of the simulations and was used to relate the observational data with our analysis using the EVS approach.
We have found that models with a bumpy spectra have higher abundance of haloes which could host first galaxies observed at z > 10.At smaller redshifts the difference in the mass function with respect to the ΛCDM model vanishes.Models with the bumpy power spectrum can better explain the results recently obtained from the observations of JWST (Naidu et al. 2022b;Labbé et al. 2023), by allowing reasonable star formation efficiency, f * < 0.2.At the same time, in the case of ΛCDM power spectrum, several of the sources detected require f * > 0.3, and one source requires f * ≈ 1.In particular, a bump with the amplitude A = 20 and position at k0 = 7 h/Mpc can explain the presence of objects with the mass ∼ 10 10 M⊙ at z > 15.It is possible to fine-tune the bump for the particular value of the star formation efficiency, but we haven't done this exercise because the data on halo masses, to our opinion, is yet not enough reliable to derive the exact parameters of the bump from it.Our goal was to show that the models with the bumpy spectrum in principle can explain the overabundance of haloes at high redshifts and converge to the ΛCDM at low redshifts.
Nowadays a number of studies propose different modifications to the power spectrum to explain the abundance of high-z galaxies (Parashari & Laha 2023;Padmanabhan & Loeb 2023;Hütsi et al. 2023).An analysis of galaxy clustering (Muñoz et al. 2023) as well as sizes of haloes and galaxies (Demiański et al. 2023) could help to distinguish between different models, including ΛCDM.We plan to use our simulations of the bumpy spectra for such an analysis in future.A bump at k ∼ 7 − 10 h/Mpc in the power spectrum can have several observational consequences.As can be seen, e.g., from our Fig. 4, in the bumpy models haloes with masses below 10 10 M⊙ form earlier than in the ΛCDM.In particular, in the volume considered, haloes with M = 10 7 M⊙ appear at z ≈ 20 in the ΛCDM, and at z ≈ 35 in gauss 1 model.Thus, in the models with the bump the first stars would appear earlier, and there will be more time to seed and grow supermassive black holes.The next step of testing the models with the bump would be to predict the 21 cm signal from Reionization in these models (see, e.g., Muñoz et al. (2020)).
So, the considered class of cosmological models can be responsible for early dark matter halo formation and results in early star formation and galaxy nuclei activity.
Figure1.Top panel: Gaussian transform functions, applied to the ΛCDM spectrum.Bottom panel: resulting matter power spectra used in our simulation suite.On both panels grey dashed lines represent the Nyquist frequency for (5 M pc/h) 3 cube and 512 3 particles resolution.The parameters of each individual spectrum can be found in the Table1.

Figure 4 .
Figure 4.The mass of the most massive halo in the simulation cube, multiplied by the fixed baryon fraction f b = 0.157, as a function of redshift z for varying primordial matter power spectra (see the legend in the upper right corners).Left panel presents an approximation of the masses of the most massive haloes using the theoretical (ST) mass function, depicted with thin solid lines.Right panel depicts the masses of the most massive halo in our numerical simulations at a given redshift, represented by the thick solid lines.The observational data from Table2are shown as filled cyan and pink dots with error bars.The filled light-grey and light-red areas represent the 3σ ranges for ΛCDM and gauss 1 mass approximations derived from EVS, and they are depicted in both the left and right panels.Bottom panel represents the difference between the simulation and ST masses respectively.

Table 1 .
Most relevant parameters of the simulation suite.

Table 2 .
Comparison of several observed galaxy candidates with ΛCDM simulations, sorted by Redshift (z).