The formation of the first cosmic structures in hybrid Milgromian cosmology

Introduction

Over the past three decades, thanks to large galaxy redshift surveys, cosmology has been transformed into a quantitative branch of astrophysics. The extensive and high-precision observational legacy has enhanced our understanding about the distribution of matter on extragalactic scales. Extracting the encoded cosmological information enables us to constrain fundamental parameters and to examine the viability of various theoretical models. Projects such as the Baryon Oscillation Spectroscopic Survey (Dawson et al. 2013, BOSS), being part of the Sloan Digital Sky Survey (Greisel et al. 2013, SDSS), probe the clustering of matter using 1.5 million luminous galaxies by measuring the Baryon Acoustic Oscillations (BAO). Built upon it, the extended Baryon Oscillation Spectroscopic Survey (Dawson et al. 2016, eBOSS) program, using the optical 2.5m Sloan Foundation Telescope at the Apache Point Observatory in New Mexico (Gunn et al. 2006), expanded the selection of luminous galaxies in the redshift range $0.7 < z < 0.9$, investigating the history of dark energy.

The two-point correlation function $\xi(r)$ quantifies the spatial correlation of matter-density fluctuations at a separation $r$, relative to an unclustered distribution, while its Fourier transform, the matter power spectrum $P(k)$, describes the amplitude of these fluctuations as a function of wavenumber $k$. According to the conventional cosmological model, primordial density fluctuations began as (almost) scale-invariant and adiabatic in the early Universe (Harrison 1970; Zeldovich 1972; Peebles & Yu [80]) and resulted from an exponentially rapid expansion phase during cosmic inflation. The initial tiny perturbations are characterized by a power law, $P(k) \propto k^{n_s}$, $n_s \approx1$ and $k$ being the comoving wavenumber. The amplitude grows as a function of time through gravitational effects, leading to non-linear evolution, giving birth to structures like galaxies or groups we observe today.

By measuring temperature anisotropies in the $10^{-5}$ level of the non-uniform blackbody spectrum, all-sky surveys such as WMAP (Bennett et al. 2013) and Planck (Tristram et al. 2024) have captured the density inhomogeneities imprinted in the surface of the last scattering on behalf of the cosmic microwave background (CMB) radiation. The widely accepted $\Lambda$CDM cosmological model (Efstathiou et al. 1990; Ostriker & Steinhardt 1995) now relies heavily on the CMB power spectrum to explain the clustering of the large-scale structure (but see Gjergo & Kroupa (2025); Kroupa et al. (2026)). Despite the successes, significant tensions against the concordance model persist. Matter bulk flows (Watkins et al. 2023) are found to be highly improbable in the $\Lambda$CDM model, as well as an immense under-dense region of $\approx600$ Mpc (Haslbauer et al. 2020; Mazurenko et al. 2024; Mazurenko et al. 2025; Banik et al. 2026). Claims against the statistically isotropic CMB, which is used as initial conditions for cosmological simulations, are raised as well (Jones et al. 2023).

The $\Lambda$CDM architecture relies primarily on the existence of Cold Dark Matter (CDM) particles and their early decoupling from the primordial plasma, so that vast gravitational potential wells (CDM halos) constitute the sites where the baryonic gas cools down and turns into stars, forming a galaxy. Nevertheless, its particle detection remains elusive (Bao et al. 2026) and astrophysical experiments dispute its presence (Kroupa et al. 2023).

On the other hand, Milgromian Dynamics, or Modified Newtonian Dynamics (Milgrom 1983a,b,c, MOND) was proposed as alternative to explain the flat galactic rotation curves without invoking the existence of unseen dark matter. The MOND phenomenology has undergone several tests across various scales. These tests range from wide binary stars (Hernandez et al. 2024; Hernandez & Kroupa 2025; Chae & Yoon 2026) to star clusters (Pflamm-Altenburg 2025), and they cover a broad mass spectrum, from ultra-faint dwarfs to groups and clusters (McGaugh et al. 2026). In numerous instances, such as in galaxy clusters (Zhang et al. 2026b) and potentially the Bullet cluster (Zhang et al. 2026a; Hernandez 2026), the assertions regarding the virial mass discrepancy (Sanders 1999) have been mitigated, and a covariant formulation is already available (Skordis & Zlosnik 2022); however, MOND cosmological frameworks often strive to replicate the concordance $\Lambda$CDM model, thereby ultimately encountering identical challenges on cosmological scales (missing satellites (Klypin et al. 1999), the too-big-to-fail conundrum (Boylan-Kolchin et al. 2011), and early galaxy formation as observed by JWST (Naidu et al. 2026)). The MOND community continues to expand (the White Paper is forthcoming), however the foundational cosmology is still in a rudimentary phase (Milgrom 1999). Since the devoted simulations of relativistic MOND theories like AeST (Aether scalar-tensor, Bataki et al. 2024) and Khronon (Blanchet & Skordis 2024) have not yet been achieved, Milgormian cosmology settles down in a less coherent but functional route.

The hybrid cosmologies exist between the two seemingly incompatible frameworks of the $\Lambda$CDM model and MOND (McGaugh 2015). The combination of a light sterile neutrino, which acts as a HDM particle, and Milgromian gravitation, provides a valuable approach for exploring cosmology, utilizing the strengths of both phenomenological models. The $\nu$HDM cosmological model originates from Sanders’ study (Sanders 1999), who first identified the discrepancy between the estimated dynamical MOND mass and the observed masses of gas-rich galaxy clusters. While dark baryons (Milgrom 2008; Kelleher & Lelli 2024) or stellar remnants (Zhang et al. 2026b) can help reconcile this discrepancy, sterile neutrinos also play a role in fitting the third peak of the CMB power spectrum (Samaras et al. 2025) without disrupting the dominance of MOND at the scale of galaxies. $N$-body simulations (Russell et al. 2026) and hydrodynamical simulations (Wittenburg et al. 2023; Samaras 2025) have addressed the subsequent structure formation; however, several questions remain unresolved. Specifically, MOND cosmology has been further supported by Sanders (1998); Sanders (2001), who argued that the enhanced gravity compared to Newtonian gravity would facilitate the formation of bound objects especially in the early Universe. Nevertheless, direct comparisons on the matter power spectrum are currently absent from the literature. Blanchet & Skordis (2025) have successfully reproduced the $\Lambda$CDM power spectrum, $P(k)$, using purely baryonic MOND approaches.

In this study, we perform hydrodynamical simulations and develop our own numerical techniques described in § 2. We evaluate the resulting power spectrum $P(k)$ from the numerical runs and compare it to existing observations (§ 3). We present the corresponding findings in § 4, and in § 5, we conclude our study.

Simulations

The nuHDM cosmological model and its variant opt-nuHDM

The $\nu$HDM model is based on a conservative ansatz. The Big Bang singularity, the cosmic inflation, and the expansion history, dictated by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, are kept in place as in the $\Lambda$CDM picture. The crucial distinction between the $\nu$HDM models and the fiducial $\Lambda$CDM one is the substitution of the CDM particle with a light, neutral, right-handed sterile neutrino ($m_{\nu_s}^{\nu\mathrm{HDM}} = 11$ eV/$c^2$ or $m_{\nu_s}^{\mathrm{opt\text{-}\nu\mathrm{HDM}}} \approx13$ eV/$c^2$) as well as the MONDian gravity instead of the Newtonian one. Our hybrid scheme assumes that the cosmological background follows the General Relativity (GR) Friedmann equations (Angus 2009; Haslbauer et al. 2020; Wittenburg et al. 2023), while the gravitational dynamics recover the GR Poisson equation in the high-acceleration (Newtonian) branch of the quasi-static regime (Zlosnik et al. 2007; Thomas et al. 2023). At low accelerations, the quasi-static gravitational dynamics are described by a MOND-type modified Poisson equation sourced by density perturbations $\delta\rho= \rho- \bar{\rho}$ (Sanders 2001; Nusser 2002; Llinares et al. 2008; Angus & Diaferio 2011; Katz et al. 2013; Candlish 2016; Wittenburg et al. 2023). The generalized Einstein-Aether (GEA) theory of Zlosnik et al. (2007) provides a possible realization of these assumptions. For coefficients $c_1=-c_3$ and $c_2=0$, so that $c_1+3c_2+c_3=0$, and a free function satisfying $\mathcal{F}(0)=0$, the vector field features the homogeneous expansion background (their sec. V). In the static weak-field limit, they obtained a modified Poisson equation for a single gravitational potential, of the form proposed by Bekenstein & Milgrom (1984), assuming that the spatial components of the vector field are of higher than first order in the perturbations (their sec. III). For the particular free function chosen in their sec. VI, this equation reproduces the ‘simple’ MOND interpolating function. On an expanding background, Thomas et al. (2023, secs. 4.1–4.3) showed that, within the GEA models and approximately they considered, cosmological MOND behavior can arise when $c_1+3c_2+c_3=0$. For an appropriate choice of the free function, the quasi-static gravitational operations recover the GR Poisson equation in the Newtonian branch, including in the corresponding linear perturbative regime, and yield a MOND-type Poisson equation sourced by density perturbations in the deep-MOND regime. In contrast, for $c_1+3c_2+c_3\ne0$, their analysis does not produce MONDian behavior on cosmological scales.

The MOND-type equations derived in these GEA limits have the nonlinear modified-Poisson structure associated with AQUAl (Bekenstein & Milgrom 1984), whereas our simulations employ the quasi-linear formulation (QUMOND; Milgrom 2010). The two formulations produce identical gravitational fields in configurations with one-dimensional symmetry, such as spherical systems, when their interpolating functions are appropriately related, but they are not generally equivalent for arbitrary three-dimensional mass distributions. Nevertheless, cosmological simulations using both formulations have demonstrated qualitatively similar structure-formation behavior relative to $\Lambda$CDM for the models and initial conditions considered (Candlish 2016). Thomas et al. (2023, sec. 4.3) also noted that existing MOND $N$-body codes can be used to investigate structure formation in suitable GEA cosmologies, provided that the matter content and simulation parameters are chosen consistently, but cautioned that the initial conditions used in previous studies may not correspond to those required by a fully consistent GEA cosmology. In the present work, we adopt GR-based linear transfer functions and initial conditions as a phenomenological approximation, motivated by the expectation that MOND modifications are small in the early Universe (Sanders 1998; Haslbauer et al. 2020, sec. 3.1.3). However, establishing their consistency with a particular GEA theory would require a full treatment of the linear cosmological perturbations, including the additional vector field and the sterile-neutrino component, which is beyond the scope of this work. Furthermore, the QUMOND formulation adopted here has not been derived directly from the particular relativistic GEA model considered above. Our simulations should therefore be regarded as a phenomenological study, motivated by the theoretical framework discussed above, rather than an exact realization of a specific relativistic theory.

Our results may nevertheless be relevant to relativistic completions admitting a similar cosmological MOND limit. Alternative theories, such as AeST (Skordis & Złośnik 2021) and the Khronon theory (Blanchet & Skordis 2024, 2025), contain additional gravitational degrees of freedom that can provide effective dark-matter-like contributions to cosmological evolution, including approximately dust-like behavior for suitable choices of the model functions and parameters, without requiring the sterile-neutrino component adopted here. These theories have different underlying cosmological dynamics and are not directly investigated in the present work.

We work with a particular family of $\nu$HDM variant models, derived from the optimal fit of the Planck CMB power spectrum. Firstly, the $\nu$HDM model leaves the rest of the parameters identical to the $\Lambda$CDM model ones, except the $\Omega_{\mathrm{CDM}}$, which is forced to zero, and that part of the matter budget is attributed to the $\Omega_{\nu} \approx0.26$. Secondly, in the optimized version ($H_0 \approx55.6$ km/s/Mpc, $\Omega_{\mathrm{m}} \approx0.49$, $\Omega_{\nu} \approx0.42$), there is an important lag on the onset (confirmed by using CAMB) of the reionization era ($z_{\mathrm{reion}}^{\Lambda\mathrm{CDM}} = 7.686 = z_{\mathrm{reion}}^{\nu\mathrm{HDM}} \ne z_{\mathrm{reion}}^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} = 5.941$) and an enhancement of 0.8 Gyr in the age prediction ($t_0^{\Lambda\mathrm{CDM}} = 13.7$ Gyr $\ne t_0^{\nu\mathrm{HDM}} \ne t_0^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} = 14.9$ Gyr¹), due to inflation of the total matter budget. For more detailed descriptions, we refer the reader to Samaras et al. (2025); Samaras (2025), and Samaras & Kroupa (2026) (under revision) concerning the functionalities of this particular hybrid framework.

From the particle physics point of view, developments of a Lagrangian accommodating a sterile neutrino hypothetical particle (Asaka et al. 2007; Salvio 2015) as a Neutrino Minimal Standard Model ($\nu$MSM) extension have been achieved too. Evidence for neutrino oscillations emanates from solar neutrinos (Cleveland et al. 1998) or reactor experiments (Eguchi et al. 2003), but a detection of a sterile neutrino via experiments like KATRIN (Mohanty 2023) and NO$\nu$A (Acero et al. 2025) incorporates many experimental degeneracies (mass range, mixing angles, etc.) and model dependence (3 active massless neutrinos and a sterile one). We adopt the standard production mechanism by Dodelson & Widrow (1994) via oscillations in the early Universe.

The main goal of this study is to explore the first emergent structures in the $\nu$HDM alternative models. MOND has long been proposed (Sanders 2008) as a potential resolution for the required enhancement of gravity on large scales in order to explain observations that are presumably incompatible with a scenario based on the Newtonian gravitational collapse inside the dark matter halos (standard $\Lambda$CDM model). The sterile neutrino compromise was a necessary adjustment to fit the CMB power spectrum and a way to silence the early claims for a mismatch in the dynamical mass of galaxy clusters (Sanders 1999). The combination of a HDM particle and the baryonic gas collapsing under MONDian gravity added a complication that has been sparsely explored. The distinctive interplay between the MONDian baryons and the light sterile neutrino allows the exploration of the "Dark Sector" and it appears as a tempting candidate to many cosmological puzzles.

The James Webb Space Telescope (JWST) has uncovered unusually bright and rapidly growing galaxies in the early Universe (Naidu et al. 2026; Heintz et al. 2026), contradicting most existing galaxy formation models. Using the HMF package (Murray et al. 2013, Halo Mass Function calculator), we first estimate the halo mass function with respect to redshift. Assuming a baryon fraction of $f_b^{\Lambda\mathrm{CDM}} = 0.157$ [84], $f_b^{\nu\mathrm{HDM}} \approx0.157$ and $f_b^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} \approx0.15$, we use the transfer functions from Samaras et al. (2025) for the hybrid models and convert the halo mass to the stellar mass one, comparing it to the $\Lambda$CDM model (Fig. 1). The collection of high-$z$ galaxy observational data is taken from Laurson (2023, fig. 6). Assuming a star-formation efficiency (SFE) of 100% (blue line), 10% (orange line), and any value below which is likely in the $\Lambda$CDM model (green shaded region) for the $\Lambda$CDM model, the shaded regions in Figure 1 show various thresholds of estimated maximum stellar masses for a cumulative number density of 1 object in a volume of $10^5$ Mpc$^3$. Compared to the local Universe, the 100% SFE is ten times more efficient. The ability of galaxies in the early $\Lambda$CDM universe to convert all of their gas into stars appears to have no clear explanation. If the survey volume was larger, the number density of galaxies would be smaller and in accord with the $\Lambda$CDM scenario based on abundance matching (Behroozi & Silk 2018). Conversely, larger number densities are problematic for the $\Lambda$CDM model.

The halo mass function is approximated by a Press-Schechter (Press & Schechter 1974) function, which dictates a hierarchical structure formation scenario [107]. However, both hybrid MOND models are forming gravitationally bound objects in a top-down scheme because of the sterile neutrino free-streaming property, where the Zeldovich pancakes are fragmented to smaller bodies. For the $\nu$HDM variant models, the HMF software allows the use of the transfer functions at $z = 199.0$ from Samaras et al. (2025) as inputs. The late structure formation is evident for both hybrid models. The red line corresponds to the progenitor $\nu$HDM model and the green one to the opt-$\nu$HDM cosmology. The two lines are not parallel to each other (either to the $\Lambda$CDM line) because of the difference in the total matter budget ($\Omega_{\mathrm{m}}^{\nu\mathrm{HDM}} = 0.3 \ne\Omega_{\mathrm{m}}^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} \approx0.49$) and the Hubble constant ($H_0^{\nu\mathrm{HDM}} = 67.64$ km/s/Mpc $\ne H_0^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} = 55.6$ km/s/Mpc). In early times ($z > 10$), both hybrid models produce the same maximum stellar mass estimates, but in later times ($z \lesssim8$), the opt-$\nu$HDM is expected to produce more massive stellar objects. Nevertheless, the resolution (horizontal gray line in Figure 1) of hydrodynamical simulations (§ 2.3) we have conducted is not sufficiently high to detect those stellar particles. None of the earlier studies by Wittenburg et al. (2023); Russell et al. (2026); Samaras (2025) focused on stars; therefore, this serves as further confirmation of the delayed formation of bound objects within this hybrid MOND framework.

Initial conditions

The initial conditions (IC) were generated at $z = 199.0$ by jointly employing the camb (Lewis et al. 2000) and the music (Hahn & Abel 2011) codes. The three cosmologies studied are flat ($\Omega_k = 0$), and their expansion histories conform to the FLRW metric. At high redshift ($z \approx200$)², the typical accelerations in the CMB frame exceed the MOND acceleration $a_0$. As a result, structure formation is expected to be minimally affected by MOND until $z \lesssim50$ (Haslbauer et al. 2020, sec 3.1.3). For the standard $\Lambda$CDM model, we utilize the cosmological parameters consistent with Planck Collaboration et al. [82]: $(\Omega_\Lambda, \Omega_m, \Omega_b, \sigma_8, n_s, h, Y) = (0.692, 0.308, 0.0482, 0.829, 0.961, 0.678, 0.24)$. In the case of the $\nu$HDM and opt-$\nu$HDM frameworks, the cosmological parameters are sourced from Wittenburg et al. (2023) for $\nu$HDM and from Samaras et al. (2025) for opt-$\nu$HDM, respectively. Samaras (2025) and Samaras & Kroupa (2026, under revision) have outlined the analytical method for generating density and velocity fields in the early Universe regarding these hybrid MOND cosmologies. Therefore, we present a summary of the cosmological parameters in Table 1.

Model	$\nu$HDM	opt-$\nu$HDM
$\omega_\nu= \Omega_\nu h^2$	0.1205	0.1307 $\pm$ 0.001
$\omega_c = \Omega_c h^2$	0.000645	0.00001
$\omega_b = \Omega_b h^2$	0.022383	0.022383
$\Omega_m$	0.3	0.495 $\pm$ 0.009
$H_0$ [km s$^{-1}$ Mpc$^{-1}$]	67.64	55.64 $\pm$ 0.32
$A_{\mathrm{Pl}}$	1.0	1.0 $\pm$ 0.0025
$n_s$	0.966	0.872 $\pm$ 0.0037
$\tau$	0.0543	0.0375 $\pm$ 0.0069
$\ln(10^{10}A_s)$	3.0448	2.987 $\pm$ 0.015

Hydrodynamical simulation

We perform the $\Lambda$CDM simulations utilizing the Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002) and the MONDian ones with the customized version of the RAMSES code called PANTOM of RAMSES (Lüghausen et al. 2015; Nagesh et al. 2021) (hereafter PoR). We chose the minimal refinement parameter, which sets the size of the coarse grid to $l_{\min} = 8$, whereas the maximum level of refinement was set to $l_{\max} = 9$.³ One should note that we used relatively low resolution; however, simulations suites for cosmological studies like Quijote (Villaescusa-Navarro et al. 2020) or Aemulus (DeRose et al. 2019) have smaller or comparable resolutions to our simulations. The runs commence at $z = 199.0$ for a $200\,h^{-1}\,\mathrm{Mpc}$ box, and star-formation is not activated. Importantly, we do not follow the Wittenburg et al. (2023) formalism in identifying bound structures in order to understand the binding procedure (§ 3.4).

Analysis of simulations

Matter power spectra

We calculate the total matter power spectrum by taking the fast Fourier transform of the dimensionless field $\delta_{\mathrm{m}}$, which is usually called the density contrast as

Considering that the total matter is composed of collisionless matter (dark matter in case of $\Lambda$CDM and sterile neutrinos in $\nu$HDM variant models) and baryons, we firstly prepare a 3D density array by interpolating the particles onto the $N = 256^3$ uniform grid via cloud-in-cell (CIC) method. Then, we sum those components as

to have one total matter density cube $\rho_{\mathrm{m}}$.

Estimation of the density field on a grid

We construct 2D projections of matter over-density $\Delta_i$ at the grid point $i$ given by:

where $\rho_i$ is density at the grid point $i$ and $\overline{\rho}$ is the average density. For particle distributions (CDM for the $\Lambda$CDM model and sterile neutrinos for the hybrid MOND models), we prepare a density cube by interpolating the particles onto the $N = 256^3$ uniform grid via the CIC method ((Hockney & Eastwood 1988; Angulo & Hahn 2022)), which is a common numerical method for particle-grid interpolation. Note that the choice of interpolation scheme was a compromise between accuracy and computational expense. In case of the hydro variables (e.g. pressure and density), we calculated the Cartesian cube using the internal RAMSES routine amr2cube. The structure and evolution of cosmic web can be characterized by the probability distribution function (PDF) of the matter density field (Einasto et al. 2021).

Estimation of the temperature field on a grid

Due to our computational resources, we do not include radiative heating and cooling, and therefore we do not have information about the ionization fractions. However each RAMSES and PoR output snapshot contains the information about the thermodynamic pressure $P_{\mathrm{cell}}$. The gas internal energy can be derived from pressure as

where $\gamma= 5/3$ for monoatomic gas composition. Then the gas temperature can be derived from internal energy as

where $k_{\mathrm{B}}$ is the Boltzmann constant, $\mu$ denotes the mean molecular weight and $m_{\mathrm{p}}$ is the mass of a proton. Note that this derivation requires the mean molecular mass of the gas to be known. This depends on the ionization state, which we do not store. However, we can calculate it using following very crude approximations

and

where $Y = 0.24$ is the helium mass fraction ((Gaikwad et al. 2017; Galárraga-Espinosa et al. 2021)).

Finding gravitationally bound structures

The cosmogonic action of the $\nu$HDM variant models is intricate due to the cooperation of the sterile neutrinos and the MONDian gravity. The cosmic web displays a well-defined filamentary structure, with large galaxy clusters situated at the nodes. The overall abundance aligns with observations (Samaras & Kroupa 2026, under revision); however, earlier studies ((Angus & Diaferio 2011; Katz et al. 2013)) indicated late structure formation and the emergence of excessively massive clusters, with total masses around $M_{\mathrm{tot}} \simeq10^{17}M_{\odot}$.

Unlike earlier studies, we do not use Amiga Halo Finder ((Knollmann & Knebe 2009)), AHF) to detect gravitationally bound objects. Bias against low-mass structures is reported ((Wittenburg et al. 2023), sec 2.5). Furthermore, according to the Halo-Finder Comparison Project (Knebe et al. 2013), agreement was obtained between different codes regarding the total, dark and stellar components of the bound objects. However, the diffuse gas content of the haloes shows great disparity, especially for low-mass satellite galaxies. Hence, we decide to adopt two different halo finding techniques⁴:

• The Friends-of-Friends (FoF) method ((Efstathiou et al. 1985; Turk et al. 2011)): the basic principle of this method is that we firstly calculate the average inter-particle spacing from the total number of particles and the volume of the region. Then, the pairs of particles closer than the inter-particle distance (compared to overall average) $l_{\mathrm{FoF}} = 0.2$ are linked together. Particles can be paired with more than one other particle. The final groups are formed the networks of particles linked together by friends.

• The HOP algorithm ((Eisenstein & Hut 1998; Turk et al. 2011)): We begin by estimating the local density at each particle using a smoothing kernel. Then the algorithm builds chains of linked particles by hopping from one particle to its densest neighbor. A particle which is its own densest neighbor is the end of the chain. All chains that share the same densest particle are grouped together. Groups are included, linked together, or discarded depending on the over-density threshold parameter $\delta= 160$.

Note that the chosen values of $l_{\mathrm{FoF}}$ and $\delta$ in case of FoF and HOP method, respectively, are typical for such applications (see (Turk et al. 2011; Smith et al. 2025)).

Results and discussion

Matter power spectra

We highlight the total matter power spectra of the three cosmologies at $z = 0$, comparing them with a compilation of data⁵ from Planck (Planck Collaboration et al. 2020a) - DES Y1 (Troxel et al. 2018) - SDSS (Reid et al. 2010) - eBOSS (Abolfathi et al. 2018) in Figure 2. This analysis underscores the enhanced power of the hybrid MOND cosmogony. The simulated surplus of power arises from the interplay between sterile neutrinos and MOND gravity. The overall shape of the power spectrum $P(k)$ deviates from that of the $\Lambda$CDM model across all corresponding scales. In both the $\nu$HDM and opt-$\nu$HDM models, the power spectrum exhibits a "knee" in the range of $10^{-1}\,h\,\mathrm{Mpc}^{-1} < k < 10^{0}\,h\,\mathrm{Mpc}^{-1}$, followed by a steeper decline at smaller scales ($k > 10^{0}\,h\,\mathrm{Mpc}^{-1}$). At even smaller scales ($k \gtrsim10^{0.5}\,h\,\mathrm{Mpc}^{-1}$), the hybrid cosmologies show a significant suppression of power, attributed to the free-streaming of sterile neutrinos and the subsequent erasure of the evolving cosmic structure. The $\Lambda$CDM model, which is based on a collisionless cold dark matter (CDM) component, exhibits significantly more power at the scales of $1\,h\,\mathrm{Mpc}^{-1} < k < 6\,h\,\mathrm{Mpc}^{-1}$. This overabundance of power—approximately a factor $\approx20 - 30$—occurs at the scales where structure formation takes place, specifically in the range $10^{-2}\,h\,\mathrm{Mpc}^{-1} < k < 10^{0}\,h\,\mathrm{Mpc}^{-1}$. This phenomenon helps to account for the presence of massive bound objects and the relative emptiness of void regions in hybrid scenarios. However, it is important to note that the observations are fitted to $P(k)$ templates of the $\Lambda$CDM model, making them not entirely independent (Tegmark & Zaldarriaga 2002).

The matter power spectra from hydrodynamical simulations evolved with RAMSES and PoR codes are illustrated in Figure 3. In this plot, we use blue to represent the $\Lambda$CDM model, red for the $\nu$HDM model, and green for the opt-$\nu$HDM model, focusing on the redshift range of $0 < z < 5$. At early times ($2 < z < 5$), a distinct pattern emerges: at very large scales $k \lesssim10^{-1} h\ \mathrm{Mpc}^{-1}$, the hybrid MOND models exhibit more power than the $\Lambda$CDM model due to their top-down structure formation approach. The emerging objects form through monolithic collapse and subsequent fragmentation. At lower scales ($k \gtrsim10^{0} h\ \mathrm{Mpc}^{-1}$), the effects of sterile neutrinos become apparent, as structure formation experiences significant delays compared to standard cosmology. The free-streaming nature of the right-handed HDM component suppresses the enhanced gravitational effects of MOND, resulting in the absence of bound structures at redshifts higher than 3. This suggests that the hybrid MOND cosmologies could prove to be unproductive of high-z galaxies.

On the other hand, for later times ($0 \leq z \leq2$), the pattern of $P(k)$ differs in the upper row of Figure 3. The hybrid MOND cosmologies continue to exhibit more power on large scales $k \lesssim10^{-1} h\ \mathrm{Mpc}^{-1}$. However, on smaller scales, the evolution of cosmic structure becomes more complex. The three distinct cosmologies show roughly the same power at $z = 2$. As time progresses, fragmentation gives way to gravitational collapse under MOND gravity. Sterile neutrinos, which are no longer relativistic and can now cluster, work together with baryonic gas to begin forming bound objects. This process creates over-densities that increase the power compared to the $\Lambda$CDM model. It is rather unclear what the power spectrum would be in the absence of sterile neutrinos, but their combination with MOND exceed the $\Lambda$CDM model estimates at low-$z$. At very small scales ($k \approx10^{0.5} h\ \mathrm{Mpc}^{-1}$), the sterile neutrino effects become more apparent, suppressing the power to a lower than the $\Lambda$CDM level. This is a late-time effect of their free-streaming inside the simulation volume. Note that we also plot the redshift-dependent total matter power spectra in Appendix § 8 for all the studied cosmologies.

Density Projections

Figure 4 presents projections of the cold dark matter (CDM) component for the standard $\Lambda$CDM model, projections of sterile neutrinos over-densities for the hybrid models. The left column represents the standard $\Lambda$CDM cosmology, the middle column depicts the $\nu$HDM model, and the right column illustrates the opt-$\nu$HDM cosmology. Cosmic evolution progresses downwards to lower redshifts. A visual examination of the snapshots at $z = 3$ indicates that the fiducial model is more fine-grained compared to the coarser hybrid counterparts. At early times, the $\nu$HDM model demonstrates the highest level of homogeneity, displaying minimal over-density values. The opt-$\nu$HDM model appears to exhibit the largest regions of emptiness. As time advances, the $\Lambda$CDM cosmology becomes increasingly detailed, revealing numerous micro-structures. In contrast, the hybrid MOND models develop distinct structures relatively early. At $z = 2$, the $\nu$HDM over-densities are widespread throughout the simulation volume, indicating a potential for significant cosmogonic activity. The opt-$\nu$HDM model at intermediate redshifts ($1 \lesssim z \lesssim2$) appears more homogeneous than its predecessor, featuring a higher matter budget, with $\Omega_{\mathrm{m}}^{\nu\mathrm{HDM}} = 0.3$ compared to $\Omega_{\mathrm{m}}^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} = 0.49$. At $z = 0$, the $\Lambda$CDM model exhibits a neat and smooth development of the cosmic web, matching large-scale observations (Jones et al. 2015). The sharp hybrid cosmologies lead to the formation of massive bound objects and a distinct cosmic web, characterized by long, notable filaments and smaller galactic systems arranged along the spines. Additionally, the extended under-dense regions observed in the opt-$\nu$HDM simulation at $z = 0$ align with theoretical predictions of MOND (Famaey & McGaugh 2012, sec. 9). Future research will explore the alignment of galactic spins within the filamentary sheets of MOND simulations.

Furthermore, Figure 5 illustrates the over-densities of baryonic matter as a function of redshift. In general, baryons tend to follow the distribution of dark matter, regardless of the type. The three cosmological models examined demonstrate significantly different patterns of gaseous evolution. The $\Lambda$CDM model generates structure on small scales, primarily through the merging process. On large scales, the structure in the $\Lambda$CDM model appears fairly uniform at redshifts $z < 2$. The cosmic baryonic web develops gradually into larger structures, characterized by small voids and less prominent filaments. Similar to the distribution of sterile neutrinos, the emergence of cosmic bodies in the MONDian gaseous model is significantly delayed compared to the $\Lambda$CDM model. Fragmentation of large clouds takes place around $z \approx2$ (see also § 4.4). Hybrid MOND models develop massive bound structures earlier than those in the fiducial cosmology, particularly evident in the panels at $z = 1$. Additionally, the overall scale of homogeneity is achieved at higher scales in the MOND models, which lack the finer and more intricate structure characteristic of the $\Lambda$CDM cosmology. The opt-$\nu$HDM model (shown in the right column) appears to be the most affected by sterile neutrino streaming, as indicated by its increased sterile neutrino budget, with $\Omega_{\nu}^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} \approx0.42$ compared to $\Omega_{\nu}^{\nu\mathrm{HDM}} \approx0.26$. The baryonic gas in the hybrid MONDian model appears to be similar, with only a few specific locations where matter accumulates, leaving the rest of the cosmological model quiet and empty compared to the $\Lambda$CDM cosmic web, with extended filaments.

Temperature Projections

We calculate the gas temperature and the gas over-densities in the $T-\Delta$ (log) plane (Figure 6). The three examined cosmologies display distinct temperature maps (upper panel of Figure 6). The $\Lambda$CDM model (left column) features numerous dispersed temperature peaks located at the nodes of cosmic filaments, indicating the presence of massive clusters of galactic systems. The temperature range is relatively narrow compared to that of the hybrid MOND models ($-2.5 \lesssim\log(T^{\mathrm{ADM}}/\mathrm{K}) \lesssim7.5$). The $\nu$HDM model (middle panel) is the least dense of the three models, featuring a few isolated regions with higher temperatures. The filamentary structure appears less extensive compared to that of the $\Lambda$CDM model. In contrast, the opt-$\nu$HDM model (right column) exhibits long, thick regions of elevated temperature. Additionally, in several instances, the temperature maps are not oriented along the typical axis of the filaments. The highly over-dense regions exhibit higher temperatures due to hydrodynamical processes; we will refer to this heating as shock heating resulting from structure formation. However, it is important to note that this numerical setup does not encompass all relevant physics, including the photo-heating of low-density gas by the ultraviolet background (UVB) to higher temperatures, either any cooling processes (Ondro et al. 2025). Simulations that incorporate additional heating and cooling processes (Hui & Gnedin 1997) indicate that temperature and density maintain a reasonably tight correlation for mildly non-linear densities, specifically for $\Delta\leq10$ (Gaikwad et al. 2017, 2018). This tight correlation is also evident in the $\nu$HDM and opt-$\nu$HDM models. Therefore, we intend to investigate the properties of the gas which surrounds cosmic filaments in the context of the $\nu$HDM and opt-$\nu$HDM models in future research.

The lower panels of Figure 6 show the dispersion of temperature values regarding the three cosmological models. The standard cosmology develops its structure in a hierarchical manner, resulting in a limited range of density contrast on the $x$-axis. In contrast, the hybrid MOND cosmologies, which combine sterile neutrinos with the MONDian enhancement of gravity, permit a broader spectrum of temperature and density inhomogeneities.

Evolution of bound objects

Wittenburg et al. (2023); Samaras (2025) have detected objects at higher redshifts, $z \approx4.5$ as spherical over-densities. Contrary to previous studies, we investigate the time evolution of the number of bound objects, as shown in Figure 7 applying different numerical methods. It is important to note that our resolution was insufficient for identifying subhaloes; therefore, only individual haloes are analyzed. We utilize two distinct algorithms: the FoF in the upper panels and the HOP one in the lower panels. At the left of Figure 7 is the $\Lambda$CDM model; the middle panel displays the $\nu$HDM model, while the right panel corresponds to the opt-$\nu$HDM model. The $\Lambda$CDM model exhibits a greater number of low-mass objects compared to the hybrid models, which aligns with expectations based on the hierarchical nature of structure formation. The hybrid MOND models exhibit slow growth in a top-down scenario; however, after $z \approx2$, galaxy formation becomes quite intense for the hybrid MOND scenarios, resulting in a shallower distribution of objects. When the sterile neutrinos cease to suppress power, the growth of structure becomes very energetic. Between the two hybrid models, the FoF algorithm indicates a lower number of low-mass bound objects for the opt-$\nu$HDM model, while the counts of massive clusters are similar.

The HOP technique indicates that low-mass objects ($< 10^{12}M_\odot$) are absent in all three cosmogonic scenarios. The $\Lambda$CDM model reveals an earlier emergence of massive objects at around $z \approx5$, accompanied by a sharp increase in their quantity at lower redshifts. In contrast, the hybrid MOND models generate massive objects rapidly, and their quantities remain stable throughout cosmic evolution, resulting in a shallower distribution. In the low-resolution setup, none of the three cosmological models examined succeed in forming objects early enough to correspond with JWST observations (Yan et al. 2022; Harikane et al. 2025; Naidu et al. 2026). Future research will focus on increasing the particle count and enhancing resolution to further explore high-redshift objects.

Last, we report the most massive halos in the simulations. Based on the FoF algorithm:

• $\max M_{\mathrm{Tot}}^{\Lambda\mathrm{CDM}} = 2.25 \times10^{15}M_\odot$

• $\max M_{\mathrm{Tot}}^{\nu\mathrm{HDM}} = 7.6 \times10^{16}M_\odot$

• $\max M_{\mathrm{Tot}}^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} = 1.07 \times10^{17}M_\odot$

Based on the HOP technique:

• $\max M_{\mathrm{Tot}}^{\Lambda\mathrm{CDM}} = 1.75 \times10^{15}M_\odot$

• $\max M_{\mathrm{Tot}}^{\nu\mathrm{HDM}} = 5.6 \times10^{16}M_\odot$

• $\max M_{\mathrm{Tot}}^{\mathrm{opt}\text{-}\nu\mathrm{HDM}} = 7.0 \times10^{16}M_\odot$

The mass estimates of the most massive bound structure are enhanced when using the FoF technique. We have revisited the well-known issue of the $\nu$HDM model, which produces unrealistically massive clusters (Katz et al. 2013; Wittenburg et al. 2023). However, Samaras (2025, chap. 3) noted that by employing more complicated star-formation recipes and incorporating supernova feedback, the total mass of the most massive halos decreases by $\approx30\%$.

Over-density & Velocity distributions at $z = 2$

This section connects the previous one to the velocity distribution of three cosmological models by examining the probability density function (PDF) of the gas over-densities at $z = 2$ for the three studied cosmologies collected from smooth coherent filaments.

The $\Lambda$CDM model exhibits an approximately log-normal distribution (Kofman et al. 1994), with its maximum at $\Delta_{\max}^{\Lambda\mathrm{CDM}} \approx0.45$. In contrast, we identify a bimodal distribution for the hybrid MOND models: in both models, the first peak is at $\Delta\approx0.65$–$0.7$ and the second at $\Delta\approx0.85$, where it appears as a shoulder in the histograms of Figure 8. All three PDFs are shown for $0 < \Delta< 2.5$ with the same binning. This may indicate the co-existence of sterile neutrinos and MONDian gravity.

Katz et al. (2013, fig. 4) study the galaxy cluster bulk flows resulting from $N$-body MOND simulations. They discovered significant MOND velocities. In our analysis, we confirm their findings by comparing them to the bulk velocities of the $\Lambda$CDM model, as shown in the middle and right panels of Figure 8. The first row represents the $\Lambda$CDM model, the middle row corresponds to the $\nu$HDM model, and the third row illustrates the opt-$\nu$HDM model. The velocity contrasts for the hybrid MOND models display distinct characteristics. The middle panel of Figure 8 indicates that the velocity behavior (taken from a $\approx0.8\,h^{-1}\,\mathrm{Mpc}$ slice from a $200\,h^{-1}\,\mathrm{Mpc}$ simulation box along the x-axis ($v_x$)) of the three cosmologies mirrors that of the density distributions shown in Figure 5. Similarly to the density projection, the velocity projection for the hybrid MOND models reveals extended high-velocity fields adjacent to low-velocity regions. The velocity field of the $\Lambda$CDM model appears more homogeneous, with less pronounced peaks in velocity and a broader range of velocity contrasts. In contrast, MOND is predicted to generate higher velocities than the $\Lambda$CDM model, which raises concerns regarding consistency with standard cosmology (Tully et al. 2023; Watkins et al. 2023). The outflows can be associated with a large local void (Haslbauer et al. 2020; Mazurenko et al. 2024). This void creates a potential hill, leading to galaxies flowing away from it and potentially addressing the Hubble tension (Banik & Samaras 2024; Di Valentino et al. 2025).

We examine each velocity component and present $v_x,v_y,v_z$ in the right panel of Figure 8. The Gaussian $\Lambda$CDM probability density function (PDF) distribution is consistent across all $v_i$ components. In contrast, the PDF for the hybrid models is significantly smaller and demonstrates anisotropic behavior. It remains to be determined whether this discrepancy is an apparent phenomenon or a prediction of Modified Newtonian Dynamics (MOND). The tails of the hybrid model’s distribution exhibit a significantly wider range of values ($-2500 \lesssim v^{\nu\mathrm{HDM}}_{\mathrm{pec}}[\mathrm{km/s}] \lesssim2500$), which aligns with the expectations from MOND boost gravity. This discrepancy, in physical terms, results in the gas being drawn into the developing large-scale structure, thereby aiding in the clearing of voids.

Conclusions

We have conducted hydrodynamical simulations based on the standard model of cosmology, $\Lambda$CDM, with parameters $H_0 = 67.8$ km/s/Mpc, $\Omega_m \approx0.31$, and $\Omega_{\mathrm{CDM}} \approx0.26$. We compared these simu lations with MONDian hybrid models, specifically the $\nu$HDM, which has parameters $H_0 = 67.64$ km/s/Mpc, $\Omega_m = 0.3$, and $\Omega_\nu\approx0.26$, as well as the opt-$\nu$HDM model, characterized by $H_0 \approx55.6$ km/s/Mpc, $\Omega_m \approx0.49$, and $\Omega_\nu\approx0.42$. The overall conclusion of this work is that sterile neutrinos present significant challenges for the hybrid MOND models. Although they provide a means to fit realistic initial conditions, such as the Planck CMB power spectrum, they also impose limitations on the formation of the first bound structures (Figure 1) in cosmic evolution. The free-streaming phenomenon suppresses ongoing galaxy formation, acting as a hindrance to the enhanced gravitational effects predicted by MOND. The combined impact of these factors results in an excess of overall power across collapse scales, leading to an overestimation compared to the $\Lambda$CDM power spectrum, $P(k)$ (Figure 2).

While this issue alone is not problematic, the observational data from the total matter power spectrum (Figure 2) derived from Planck, SDSS, DES Y1 and eBOSS are all model-dependent, similar to the Planck CMB temperature fluctuation power spectrum (Samaras et al. 2025, sec. 5.5). Nonetheless, a significant concern regarding the hybrid models under investigation is the delayed emergence of bound objects (Figure 7). The recent discoveries of high-redshift galaxies by JWST (Harikane et al. 2025; Naidu et al. 2026) raise important questions for any coherent cosmological framework, although the resolution of our simulations is not sufficient to search for low-mass early Universe structures. Additionally, the sterile neutrino serves to counterbalance the gravitational effects enhanced by MOND. More broadly, the various $\nu$HDM models remain viable contenders due to the limitations in resolution present in this analysis. The density and velocity fields in the hybrid MOND models are significantly more

coarse-grained compared to the finer-grained fields in the $\Lambda$CDM models. Additionally, the temperature and over-density ranges in the $\nu$HDM variant models are broader, indicating an enhancement of the gravitational fields due to MOND. The hybrid models’ total matter power spectrum at $z \approx2$ surpasses the predicted power of the $\Lambda$CDM model at large scales. In light of this, we analyze each of the velocity components, revealing an anisotropy. We remind the reader that more advanced hybrid simulations are available and appear to refine the estimates of the most massive hybrid objects at $z = 0$. Future research will focus on star-formation techniques, proper photo-heating (by UV background) and resolution effects.

Acknowledgments

The authors would like to express their gratitude to the Czechoslovakian student conference series (BEZOVEC), where they met and initiated their collaboration. The authors would also like to thank Tom Złośník, Bharat Ratra, Chan-Gyung Park and Peter Laursen for the fruitful discussions. NS received support from the Charles University Grant Agency (GAUK) - 94224, as well as from the Ostrava Supercomputing Center for the project OPEN-32-11 #2885, where the $\nu$HDM simulations were conducted. NS received support from the Dynamix computer system at the Argelander Institute for Astronomy, Bonn University, as well as from the SPODYR group. Additionally, NS benefited from assistance through the Deutscher Akademischer Austauschdienst (DAAD) Bonn-Prague exchange program. Currently, NS is employed at the University of Ottawa in Ottawa, Ontario, Canada.

Data Availability

We acknowledge the use of the Planck Legacy Archive. Planck (http://www.esa.int/Planck) is a European Space Agency (ESA) science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. The Sloan Digital Sky Survey (SDSS) is a redshift survey using a dedicated 2.5-m wide-angle optical telescope at Apache Point Observatory in New Mexico, United States. We make use of the Year 1 data release of the Dark Energy Survey (DES). eBOSS is Spectroscopic Survey studying the Baryon Acoustic Oscillations. For reproducing this work, one will have to visit NS’s GitHub page and TO’s page.

Redshift evolution of $P(k)$

The evolution of the total matter power spectrum is shown in Figure A1. The hybrid scenarios surpass the $\Lambda$CDM model in power at very large scales ($k < 10^{-1} h\,\mathrm{Mpc}^{-1}$) throughout all phases of structure formation. In late times ($z < 2$), the total matter power spectra of the hybrid cosmologies depend on the interplay between MOND gravity and sterile neutrinos. Except for the very late times ($z = 0$), the power systematically increases for both the $\nu$HDM model and the opt-$\nu$HDM model. It is only at the final snapshot ($z = 0$) that the sterile neutrinos significantly suppress the power at scales $k \gtrsim10^{0.5} h\,\mathrm{Mpc}^{-1}$. The cause of this effect—whether it stems from numerical resolution, free-streaming activities of sterile neutrinos, or a combination of both—will be investigated through simulations with varying box sizes. Nonetheless, this effect is particularly pronounced in the hybrid cosmologies, as the $\Lambda$CDM model consistently enhances the matter power spectrum $P(k)$ as redshift $z$ decreases.

This paper has been typeset from a \TeX/\LaTeX file prepared by the author.