Hydrodynamical structure formation in Milgromian cosmology

ABSTRACT

We present the first hydrodynamical cosmological simulations in the νHDM framework based on Milgromian dynamics (MOND) with light (11 eV) sterile neutrinos. νHDM can explain the expansion history, CMB anisotropies, and galaxy cluster dynamics similarly to standard cosmology while preserving MOND’s successes on galaxy scales, making this the most conservative Milgromian framework. We generate initial conditions including sterile neutrinos using camb and music and modify the publicly available code phantom of ramses to run νHDM models. The simulations start at redshift |$z_e$| = 199, when the gravitational fields are stronger than |$a_{_0}$| provided this does not vary. We analyse the growth of structure and investigate the impact of resolution and box size, which is at most 600 comoving Mpc. Large density contrasts arise at late times, which may explain the KBC void and Hubble tension. We quantify the mass function of formed structures at different redshifts. We show that the sterile neutrino mass fraction in these structures is similar to the cosmic fraction at high masses (consistent with MOND dynamical analyses) but approaches zero at lower masses, as expected for galaxies. We also identify structures with a low peculiar velocity comparable to the Local Group, but these are rare. The onset of group/cluster-scale structure formation at |$z_e$| ≈ 4 appears to be in tension with observations of high redshift galaxies, which we discuss in comparison to prior analytical work in a MONDian framework. The formation of a cosmic web of filaments and voids demonstrates that this is not unique to standard Einstein/Newton-based cosmology.

1 INTRODUCTION

Galaxies, groups of galaxies, and galaxy clusters form and evolve under cosmological boundary conditions, which define their seed mass, interactions with each other, and accretion of gas from their surroundings. Most available self-consistent simulations of these processes have been made in the standard model of cosmology (SMoC), which assumes the universal validity of Einstein’s general theory of relativity. To match observations of the expanding Universe, the cosmic microwave background (CMB) anisotropies, and the properties of nearby galaxies, the SMoC needs to be augmented by the auxiliary hypotheses that dark matter particles and dark energy dominate the matter and energy content of the Universe. Although the latter can be accounted for by a cosmological constant Λ, which is already a part of General Relativity, modified gravity theories (e.g. Clifton et al. 2012; Baker et al. 2021) are usually considered in simulations only as an alternative to Λ (e.g. Li et al. 2012). Such models rarely question the hypothesis that most of the matter in the Universe and in galaxies is made of exotic dark matter particles which interact gravitationally but not electromagnetically and that are not accounted for in the standard model of particle physics. The exotic dark matter is generally assumed to consist of cold dark matter (CDM) particles that are non-relativistic at decoupling and interact only gravitationally with baryons, which combined with the cosmological constant and various other assumptions leads to the Lambda cold dark matter (ΛCDM) paradigm (Efstathiou, Sutherland & Maddox 1990; Ostriker & Steinhardt 1995).

The ΛCDM simulations that have been and are being performed include, among many others, Illustris (Vogelsberger et al. 2014; Nelson et al. 2015, 2019; Pillepich et al. 2018), EAGLE (Schaye et al. 2015; McAlpine et al. 2016), FIRE (Hopkins et al. 2014, 2018), HORIZON-AGN (Dubois et al. 2014), and NEWHORIZON (Dubois et al. 2021). These are based on different codes to treat hydrodynamics and different sub-grid algorithms to account for star formation and gas heating. The CDM part of the paradigm is sometimes altered by postulating different properties for the dark matter, assuming it to be e.g. fuzzy (FDM; Hui et al. 2017; May & Springel 2021), warm (WDM; e.g. Lovell et al. 2020), or self-interacting (SIDM; e.g. Rocha et al. 2013). While this alters the small-scale behaviour of the dark matter, constraints are often so tight that these alternative universes regularly end up looking like CDM models (Rogers & Peiris 2021). For example, FDM would need to behave like CDM in order to explain the large Newtonian dynamical mass-to-light ratios of the smallest ultra-faint galaxies, while on large scales the behaviour is similar to CDM such that the large-scale structure would be similar to the ΛCDM model.

Simulations like those mentioned above have allowed testing the SMoC against observations to unprecedented depth, thereby uncovering problems which pertain to small-scale issues (e.g. Kroupa et al. 2010, and references therein), including the planes of satellites problem (Kroupa, Theis & Boily 2005), which seems particularly robust to how the baryonic physics is treated (Pawlowski 2018, 2021a,b). There are also serious problems on Gpc scales (Di Valentino, Melchiorri & Silk 2020; Haslbauer, Banik & Kroupa 2020; Migkas et al. 2021; Mohayaee, Rameez & Sarkar 2021; Secrest et al. 2022; Bargiacchi et al. 2023; Lenart et al. 2023). The El Gordo interacting galaxy cluster at redshift z_e = 0.87 causes by itself a >6σ falsification of the SMoC (Asencio, Banik & Kroupa 2021), while the Hubble tension (Haslbauer et al. 2020; see also the review by Di Valentino 2021) is similarly serious and similarly immune to resolution through improved treatment of baryonic physics on galaxy scales. Furthermore, a direct test for the existence of CDM particles on galaxy scales using Chandrasekhar dynamical friction suggests these particles do not exist at 13σ confidence based on galaxy bars not slowing down as expected (Roshan et al. 2021), while the past orbital dynamics of the M81 group also precludes solutions with dark matter haloes (Oehm, Thies & Kroupa 2017). The observed perturbations to dwarf galaxies in the Fornax galaxy cluster and the lack of low surface brightness dwarfs towards its centre also lead to the same conclusion, this time using tidal stability arguments (Asencio et al. 2022). For a comprehensive recent review of problems faced by ΛCDM, we refer the reader to Kroupa (2015) and Perivolaropoulos & Skara (2022).

Given these tensions between theory and observations, it would be of great value to have available simulations of structure formation in a fundamentally different theoretical framework. Such simulations would allow us to better appreciate the above-mentioned tensions of the SMoC. Depending on the performance of such a model, we would obtain information as to whether the particular approach is promising or should be discarded altogether. Currently, the only other theoretical approach to cosmological structure formation that can plausibly address these issues and be coded in the form of particle equations of motion is based on Milgromian dynamics (MOND; Milgrom 1983; Bekenstein & Milgrom 1984; Famaey & McGaugh 2012; Merritt 2020; Banik & Zhao 2022). The basic idea of MOND is to not restrict oneself to the Solar System-scale constraints on the behaviour of gravity that allowed Newton and Einstein to formulate their theories, but to also consider more recent results from the rotation curves of galaxies, which reveal a much richer dynamics (Faber & Gallagher 1979, and references therein). To explain these observations using only the directly detected baryonic mass in stars and gas, Milgrom (1983) conjectured that the gravitational acceleration g in an isolated spherically symmetric system is asymptotically related to the Newtonian gravitational acceleration |$g_{_N}$| of the baryons alone according to

The key new ingredient is a fundamental acceleration scale |$a_{_0} = 1.2 \times 10^{-10}$| m s⁻² ≈3.9 pc Myr⁻² (this value is empirical but has remained stable for many decades; see Begeman, Broeils & Sanders 1991; Gentile, Famaey & de Blok 2011; McGaugh, Lelli & Schombert 2016). This behaviour follows from a Lagrangian, ensuring the usual symmetries and conservation laws with respect to the linear and angular momentum and the energy (Bekenstein & Milgrom 1984; Milgrom 2010). The least action principle then leads to a generalized Poisson equation that is non-linear in the mass distribution.

Though we use a non-relativistic formulation in this work assuming a standard background expansion history, relativistic formulations exist which are consistent with the CMB anisotropies, strong gravitational lensing, and the fact that gravitational and electromagnetic waves travel at the same speed to high precision (Skordis & Złośnik 2019, 2021, 2022). MOND may follow from properties of the quantum vacuum in the presence of a cosmological constant (Milgrom 1999).

MOND has been uncannily successful, with all a priori predictions on galactic scales made in Milgrom (1983) having been verified by subsequent observations (for reviews, see Famaey & McGaugh 2012 and Banik & Zhao 2022). The first cosmological structure formation simulations that show the formation of a cosmic web were done by Llinares, Knebe & Zhao (2008), where the authors assumed a purely baryonic universe with Ω_m = Ω_b = 0.04 and treated the matter as pressureless dust, neglecting hydrodynamics. The problems of MOND on galaxy cluster scales can be addressed without needing the specific features of the scalar and vector field of the Skordis & Złośnik (2021) model by postulating sterile neutrinos with a rest energy of 11 eV (Angus 2009; Angus, Diaferio & Kroupa 2011; Angus et al. 2013). Originally the idea of adding a collisionless component to account for the missing mass in MOND on cluster scales stems from Sanders (2003), who argued at the time that ordinary neutrinos could be a possible candidate (though this is no longer possible; see Katrin Collaboration 2019, 2022). These 11 eV sterile neutrinos naturally account for the CMB similarly to ΛCDM because they have the same relic mass density as the CDM in the ΛCDM paradigm, leading to a standard expansion history and very similar physics in the early universe (for a more detailed description, see section 3.1 of Haslbauer et al. 2020). We stress here that the collisionless sterile neutrino component is not related to the exotic CDM particles in the SMoC because the postulated light sterile neutrinos would not cluster on galaxy scales (Angus 2010), though the cosmic mean density of the two is about the same.

This neutrino hot dark matter (νHDM) model¹ was implemented numerically by Angus et al. (2011) and later by Katz et al. (2013), while Haslbauer et al. (2020) used semi-analytic methods to show that it mitigates the Hubble tension by means of galaxies falling towards the void edge during the formation of a large local supervoid, which is indeed observed as the KBC void (Keenan, Barger & Cowie 2013). Turning to overdensities, Asencio et al. (2021) used the N-body νHDM simulations of Katz et al. (2013) to show that this framework also accounts for the El Gordo galaxy cluster collision, thereby producing roughly the right frequency of very massive interacting galaxy clusters at high z_e, a test which ΛCDM fails with >6σ confidence. While the Llinares et al. (2008) and the Katz et al. (2013) simulations considered only large scales and thus treated the baryons as dissipationless dust, currently no equivalent modelling exists which takes into account the hydrodynamics in order to resolve scales down to individual galaxies.

We simulate 300 and 600 comoving Mpc (cMpc)-scale cosmological volumes with gas and sterile neutrinos, starting at z_e = 199 and evolved until the present time. This volume can reach galaxy group/cluster scales in moderately powerful in-house computer servers, thereby allowing a comparison with observed galaxy clusters. We introduce the publicly available method and codes² used here and analyse the mass function of the resulting galaxy cluster population at different epochs. An important aspect of our simulated massive galaxy clusters is that they are dominated by sterile neutrinos, thus alleviating concerns about the application of MOND to the Bullet Cluster (Angus et al. 2007). Subsequent publications will report more extensive analyses of the properties of galaxy groups and clusters down to individual galaxies, a scale which we are already able to reach in this work.

Section 2 introduces Milgromian gravitation and describes the νHDM cosmological framework. Since star formation is not taken into account in this first exploratory modelling, the method to find bound structures is described in Section 2.5. The models computed here are detailed in Section 2.6 and their results are discussed in Section 3. Section 4 concludes this contribution. Movies of the simulations are publicly available.³ We emphasize that these are the first ever hydrodynamic simulations performed in a MOND cosmological model. The here-studied νHDM model is conservative in the sense that it retains a very early inflationary phase, assumes the FLRW metric to be valid, and assumes dark energy, allowing it to reproduce the CMB with its high third peak and yield the same expansion history as the SMoC.

2 THEORETICAL BACKGROUND AND NUMERICAL METHODS

2.1 MOND

MOND (Milgrom 1983) can be formulated as a non-relativistic theory of gravity with a fully developed Lagrangian formalism (Bekenstein & Milgrom 1984; Famaey & McGaugh 2012). MOND was originally developed by Milgrom (1983) to explain the missing gravity problem on galaxy scales (evidenced in their flat outer rotation curves; see Bosma 1978; Rubin, Ford & Thonnard 1978; Faber & Gallagher 1979, and references therein) without invoking the existence of exotic dark matter particles. The underlying reason for the proposed departure from standard physics may be related to the quantum physics of the vacuum (Milgrom 1999; Pazy 2013; Smolin 2017; Senay, Mohammadi Sabet & Moradpour 2021). According to MOND, the gravitational accelerations depart from the Newtonian behaviour in the MOND regime of very weak gravitational fields (⁠|$g \ll a_{_0}$|⁠), but standard dynamics is recovered in the Newtonian or strong-field regime (⁠|$g \gg a_{_0}$|⁠). Equation (1) shows the asymptotic regimes, the transition between which occurs when g is comparable to Milgrom’s constant |$a_{_0}$|⁠. To transition between these regimes and handle more complicated geometries while also basing the theory on a Lagrangian to ensure the usual symmetries and conservation laws, Bekenstein & Milgrom (1984) developed the following MOND Poisson equation:

where Φ is the potential, N subscripts denote Newtonian quantities, |$\boldsymbol {r}$| is the position, and μ is the interpolating function with dimensionless argument |$y_\mathrm{int} \equiv \vert \nabla \Phi \vert /a_{_0}$|⁠, the gravitational field in units of |$a_{_0}$|⁠. For consistency with the empirical equation (1),

In spherical symmetry, |$\mu \boldsymbol {g} = \boldsymbol {g}_{_N}$|⁠. Equation (2) was commonly used in past simulations (e.g. see the publicly available numerical implementation in raymond; Candlish, Smith & Fellhauer 2015). Nowadays, a computationally less intensive formulation is typically implemented, which has the same asymptotic behaviour (see Section 2.3).

One of the most valuable attributes of MOND is its highly predictive nature on galaxy scales (Merritt 2020). In particular, MOND predicted a tight baryonic Tully–Fisher relation with slope 4 and anticipated that the relation would be tightest between the flat part of the rotation curve and the total baryonic mass (Lelli et al. 2019). Observations show that considering other velocity scales like the maximum of the rotation curve leads to a weaker correlation, while considering only the stellar mass rather than the total baryonic mass also leads to a breakdown of the relation at the low-mass end where the contribution of gas is important (McGaugh 2012). Beyond just the flat outer part of the rotation curve, MOND predicted a tight radial acceleration relation (RAR)⁴ between the observed g and |$g_{_N}$|⁠, which is very apparent in the observational data (Sanders 1990; McGaugh 2004, 2005, 2012; Wu & Kroupa 2015; Lelli et al. 2017), including lenticular galaxies (Shelest & Lelli 2020). For reviews of this strong empirical correlation, we refer the reader to Famaey & McGaugh (2012) and McGaugh (2020). While the focus in the past has been on rotating disc galaxies, MOND has also been successfully applied to predict the velocity dispersions of isolated dwarf galaxies in the Local Group (McGaugh et al. 2021). Indeed, it has been argued in an extensive recent review that the RAR is ‘not unique to rotationally supported systems or to rotational motion, nor is it confined to thin disc galaxies’ (section 3 of Banik & Zhao 2022). The RAR also extends to weak gravitational lensing, which demonstrates its validity down to |$g_{_N} \approx 10^{-5} \, a_{_0}$| (Brouwer et al. 2021).

2.2 The νHDM model

Despite the impressive successes of MOND on galaxy scales, it cannot easily account for observations on larger scales like the dynamics of galaxy clusters (Sanders 1999, 2003) and the CMB anisotropies (McGaugh 2004; Planck Collaboration VI 2020) if we restrict ourselves to a purely baryonic universe. One recent proposal to address this issue has been the new aether scalar tensor (AEST) relativistic theory of Skordis & Złośnik (2021) in which the action has a function that plays the role of the MOND interpolating function, but its value depends on both the spatial gradient squared of the scalar field and its temporal derivative. The time-dependent term behaves like gravitating dust, allowing AEST to reproduce the angular power spectrum of the CMB and the ΛCDM result for the matter power spectrum on large scales. One can speculate that it might also give rise to additional gravitational acceleration inside galaxy clusters, explaining the residual missing mass of MOND (Ettori et al. 2019). However, choosing AEST model parameters which might achieve this seems to cause problems matching the weak lensing RAR (Mistele, McGaugh & Hossenfelder 2023). In any case, it is clear that significant fine-tuning would be required to consistently explain the dynamics of both galaxies and galaxy clusters with just their luminous mass.

A more conservative approach (Section 1) is to consider that it is perfectly possible to have both MONDian gravity and a type of collisionless matter that does not affect the MOND fits to galaxy rotation curves. Angus (2009) developed such a model by postulating the existence of sterile neutrinos with a mass of m_νs = 11 eV c⁻², where c is the speed of light. Thermal sterile neutrinos with this mass would have the same relic mass density as the CDM particles in ΛCDM: if half of all quantum states are occupied (similarly to the active neutrinos), no fine-tuning of their abundance would be necessary (see Diaferio & Angus 2012, and references therein). Combined with the high accelerations around the epoch of recombination (⁠|$g \approx 20 \, a_{_0}$|⁠) and an almost standard expansion history, this would lead to the early universe behaving similarly to ΛCDM (for a more detailed discussion, see section 3.1 of Haslbauer et al. 2020). The primordial abundances of light elements would then work much the same in both paradigms, with only small differences. The dynamics of galaxy clusters would also be accounted for in νHDM (Angus, Famaey & Diaferio 2010), including the offset between the weak lensing and X-ray peaks in the Bullet Cluster (Angus et al. 2007).

From a particle physics perspective, the fact that ordinary neutrinos interchange their masses as they propagate suggests that there is an additional (sterile) neutrino in order to conserve kinetic energy and momentum. Sterile neutrinos are thus not equivalent to cold or fuzzy dark matter particles, but were in the first instance hypothesized to exist based on the observed properties of standard particles (Merle 2017). At a rest energy >10 eV, the free streaming length of sterile neutrinos would be short enough to not discernibly affect the CMB (section 6.4.3 of Planck Collaboration XIII 2016). The Tremaine–Gunn limit to their phase space density (Tremaine & Gunn 1979) would also be low enough that the dynamics of galaxies would not be affected by the sterile neutrinos, so galaxies would still be purely baryonic and MONDian (Angus 2010). This also implies that although νHDM contains a dominant collisionless matter component like in ΛCDM, it does not promote the formation of galaxies because it is unable to efficiently cluster on very small scales. We adopt this νHDM model and implement it into the simulation code phantom of ramses (por; Lüghausen, Famaey & Kroupa 2015; Nagesh et al. 2021).

2.3 Simulation code

por is a customized version of the hydrodynamical N-body adaptive mesh refinement (AMR) code ramses (Teyssier 2002), which implements the quasi-linear formulation of MOND (QUMOND; Milgrom 2010). The novel approach is to interpret the Milgromian potential |$\Phi \left(\boldsymbol {r} \right)$| as arising from the application of Newtonian gravity to some effective density distribution

where ρ_p is the ‘phantom dark matter’ (PDM) density defined such that ∇²Φ ≡ 4πGρ_eff. The PDM is merely a mathematical construction that does not correspond to real particles and therefore does not lead to Chandrasekhar dynamical friction. The advantage of thinking in this way is that whereas equation (2) involves a non-linear Poisson equation that is quite difficult to solve numerically, we may instead solve the following system of equations:

where |$\widetilde{\nu } \equiv \nu -1$| and ν is the QUMOND interpolating function, defined such that μν = 1 in spherical symmetry. In this case, |$\boldsymbol {g} = \nu \boldsymbol {g}_{_N}$|⁠, with ν taking the argument |$y_\mathrm{int} \equiv \vert \nabla \Phi _N \vert /a_{_0}$|⁠. This is already much more amenable to computer simulations because |$\boldsymbol {g}_{_N}$| and any simple function of it are readily calculable using standard techniques. This formulation of MOND also follows from a Lagrangian (Milgrom 2010), ensuring the usual symmetries and conservation laws with respect to the linear and angular momentum and the energy. For consistency with equation (3), ν must satisfy the following asymptotic limits:

por applies the simple interpolating function (Famaey & Binney 2005) because it agrees with recent observations on galaxy scales (Iocco, Pato & Bertone 2015; Lelli et al. 2017; Banik & Zhao 2018; Chae, Bernardi & Sheth 2018; Zhu et al. 2023)⁵

QUMOND can be solved for the total potential Φ using already existing Poisson solvers. This entails solving the Poisson equation twice, first to get Φ_N and its derivatives in order to compute the source term on the right-hand side of equation (5), and then to solve this equation. This method makes a direct comparison between the Newtonian and Milgromian cases easier (for further information about the technicalities of the code, see Wittenburg, Kroupa & Famaey 2020). A user manual for por has recently been written (Nagesh et al. 2021) and contains a summary of numerical MOND simulations that used it to explore isolated and interacting disc galaxies.⁶ Some of its interesting recent applications have been to the asymmetric tidal tails of open star clusters in the Solar neighbourhood (Kroupa et al. 2022), the tidal stability of dwarf galaxies in the Fornax Cluster (Asencio et al. 2022), the formation of the Local Group satellite planes from a past flyby encounter between the Milky Way and M31 galaxies (Bílek et al. 2018; Banik et al. 2022), and the star-formation rates and bar statistics of isolated disc galaxies initialized with a realistic size for their stellar mass (Nagesh et al. 2023).

In order to perform cosmological simulations with por, we had to make several changes based on the description in Candlish (2016). The main difference to the non-cosmological application of por is that to account for the cosmic expansion, we solve the Poisson equation in supercomoving coordinates (Martel & Shapiro 1998; Teyssier 2002,; Appendix  A):

with a_e ≡ 1/(1 + z_e) being the cosmic expansion factor at redshift z_e and a tilde denoting supercomoving coordinates except for |$\widetilde{\nu }$|⁠. The parameter |$\Omega _m \equiv \bar{\rho }/\rho _\mathrm{crit}$| is the fraction that the mean matter density |$\bar{\rho }$| currently comprises of the cosmic critical density

An important detail in a cosmological context is that the MOND Poisson equation is only applied to the density contrast relative to the average. This originates from the works of Sanders (2001) and Nusser (2002), which introduced the basic principle behind the growth of density perturbations in a Milgromian framework. It has recently been argued that our approach is valid in some relativistic MOND theories (Thomas, Mozaffari & Zlosnik 2023).

The hydrodynamical evolution of the gas is calculated by solving the Euler equations, which por does using a second-order Godunov scheme. We adopt the simple cooling and heating scheme (Wittenburg et al. 2020, and references therein), which makes use of look-up tables and assumes collisional excitation equilibrium (CIE). The simulations here do not allow stellar particles to form. The more rigorous scheme with radiative transfer and supernovae (see also Wittenburg et al. 2020) will have to be implemented in future higher resolution work. This should be straightforward as it is already coded and merely needs to be activated, though cosmological MOND simulations including star formation would be more computationally intensive. Such simulations have already been done in a non-cosmological context and give rather promising results (Renaud, Famaey & Kroupa 2016; Wittenburg et al. 2020; Eappen et al. 2022; Nagesh et al. 2023), with the latter work clarifying what sub-grid parameters would be appropriate in a MOND context and also that the particular choice does not much influence the results.

2.4 Initial conditions

The setup of a cosmological simulation involves two main steps:

• Specifying a power spectrum for the density fluctuations at the desired starting redshift; and

• Sampling this power spectrum using a random number generator to specify the initial density and velocity field.

To generate the power spectrum at the desired starting redshift, post-inflationary density fluctuations δ(k) are determined based on the cosmological parameters, with the wavenumber k being the inverse of the comoving length scale of the perturbation. The ratio between the density fluctuation at some early time just after the inflationary epoch and at some later epoch is known as the transfer function T(k), which is defined so that

where δ₀ refers to the early epoch and δ to the later epoch. The power spectrum P(k) ∝ δ²(k) is the Fourier transform of the two-point correlation function (e.g. Eisenstein & Hu 1998). The transfer function needs to evolve inflation-generated density perturbations through the radiation-dominated era, recombination, and part of the matter-dominated era up to the epoch when we wish to start the simulation. To limit the computational cost, linear perturbation theory is usually applied for as long as possible, so an important issue is when the structures become non-linear. In the SMoC, this occurs when |δ| ≪ 1 is no longer a good approximation. The linear regime breaks down earlier in νHDM because the underlying gravity law is non-linear. It was argued in section 3.1.3 of Haslbauer et al. (2020) that the typical gravitational field enters the MOND regime when |$z_e \lesssim 50$|⁠, so we start our simulations at z_e = 199 when typical accelerations are still in the Newtonian regime.

The transfer function depends on various physical processes like the Mészáros effect (Mészáros 1974), acoustic oscillations, Silk damping (Silk 1968), free streaming, and radiation drag (Peacock 2003, and references therein). In a classical CDM approach, there are no oscillations in the transfer function of the CDM component because it is collisionless and pressureless. The oscillations in the baryonic component stem from the coupling with radiation in the very early universe. After recombination, these oscillations are strongly damped as baryons are forced to follow the pattern of density fluctuations in the dominant CDM component, i.e. the baryons follow the gravitational potential of the CDM (black and blue lines in Fig. 1). In contrast, oscillations are evident for both the baryons and the massive neutrinos in the νHDM model (green and red lines in this figure). This is expected for the baryons because there is no CDM and thus little power on small scales at z_e = 199. Since the massive neutrinos have very little power on small scales due to free streaming effects, oscillations in the baryonic power spectrum are able to induce similar oscillations in the sterile neutrino component due to gravitational interactions between the baryons and neutrinos.

The presence of acoustic oscillations in the dominant mass component of the universe as late as z_e = 199 is a novel feature of the νHDM paradigm. One important consequence is that the overall transfer function can become negative, indicating that an initial overdensity becomes an underdensity and vice versa. In the linear regime, a purely gravitational problem cannot produce a negative transfer function. This arises in our model and in ΛCDM at higher redshifts due to the coupling with radiation, which causes the oscillatory behaviour evident in Fig. 1. Note, however, that the sign of T(k) does not change the power spectrum or the initial density distribution for our simulations because the initial phase of a density fluctuation is random, so only |T(k)| is relevant for our purposes.

The impact of MOND on the early universe was already discussed in McGaugh (1999) and references therein. Although the results were based on a purely baryonic model, we still expect the conclusions to be valid regarding MOND producing a more rapid growth of structure and pronounced baryon acoustic oscillations (BAOs) in the transfer function at epochs much later than recombination. Additionally, we expect to detect a difference in the absorption signal of the 21 cm spin–flip transition of neutral hydrogen seen against the CMB. McGaugh (2018) compared the 21 cm absorption feature in a purely baryonic universe to the ΛCDM prediction and concluded that the signal in a universe devoid of CDM should be significantly deeper. When comparing νHDM and ΛCDM, this discrepancy will certainly be less pronounced, but in general we expect a difference in the 21 cm absorption signal between these two models due to differences in the formation of galaxies and other ionizing sources. In particular, significant BAOs should be evident in the matter power spectrum because non-linear MOND effects (which only become important at |$z_e \lesssim 50$|⁠) would presumably take some time to wash out the strong BAOs through coupling between perturbations on different scales.

The transfer function for our simulations is calculated using the camb program (Lewis et al. 2000). We adopt the following parameters: Ω_m = 0.314595, Ω_b = 0.0492713, Ω_c = 0, |$\Omega _\Lambda = 0.685405$|⁠, σ₈ = 0.8101, n_spec = 0.9660499, |$H_0 = 67.4 \, \mathrm{km s^{-1}Mpc^{-1}}$| or h = 0.674 (this is consistent with the standard parameters from Planck Collaboration VI 2020), and Ω_ν = Ω_m − Ω_b = 0.2653237 for the neutrino component. The gas temperature at z_e = 199 is T = 370 K, which follows from the CMB conditions. The main difference is that we close the gap between Ω_b and Ω_m using 11 eV c⁻² sterile neutrinos instead of CDM particles. The sterile neutrino mass does not need to be specified in camb, but it does display the calculated mass on screen assuming thermalization in the early universe. To implement the sterile neutrinos, we follow the description in Angus (2009),⁷ which leads to the addition/alteration of the following lines of the camb ‘.ini’ namelist file:

 massless_neutrinos = 2.0293

 nu_mass_eigenstates = 2

 massive_neutrinos = 1 1

 nu_mass_degeneracies = 1.0147 1

 nu_mass_fractions = 0.0044 0.9956

 accurate_massive_neutrino_transfers = T

The left-hand panel of Fig. 2 shows the temperature fluctuations as CMB multipole moments in ΛCDM with standard cosmological parameters (Planck Collaboration VI 2020) and in the νHDM model simulated here. The νHDM model fits the observations quite well (as was already shown by Angus 2009). Further refinements to the cosmological parameters are possible, but we do not consider this here. Fig. 1 of Angus & Diaferio (2011) suggests that finding an optimal fit to the CMB in νHDM would slightly increase the inferred H₀, but not by enough to solve the Hubble tension − we are currently investigating this issue (Luchoomun et al., in preparation).

The right-hand panel of Fig. 2 shows the matter power spectrum in terms of wavenumber k for the ΛCDM and νHDM models at our starting redshift of z_e = 199. The models are practically identical for large scales (small k), but the power of the νHDM model decreases significantly at scales below 100 cMpc. This behaviour is to be expected as the sterile neutrinos do not cluster on small scales by construction. If they did, then the gravitational field on galactic scales would be influenced by them, but there is generally excellent agreement between the MOND dynamical mass of a galaxy and its baryonic mass (Section 1). Another interesting feature of the νHDM matter power spectrum on small scales is that it shows significant oscillations, which are also evident in the transfer function (Fig. 1). We conclude that these oscillations are driven by the well-known acoustic oscillations in the baryons, which have comparable or more power than the sterile neutrinos on these small scales and can therefore influence them significantly through gravity. The ΛCDM power spectrum does not show the same oscillations on small scales because the CDM is dominant on these scales and dampens the BAOs, which are a relic of the era when matter was tightly coupled to radiation.

Having specified the power spectrum, the density field is obtained using the publicly available code music (Hahn & Abel 2011), but with an externally supplied transfer function (see above), bypassing its inbuilt calculator. music works with the logarithm of the transfer function, which it implicitly assumes is positive. Since this assumption is violated in νHDM, we adjust the algorithm to work with |T(k)|, which as explained above does not alter the physics. We use three different resolution settings in which l_max = l_min + 1 in all cases and l_min = 7, 8, or 9. This is done in order to prevent the appearance of artificial overdensities, since the number and mass of the sterile neutrino particles is determined by l_min: the total number of sterile neutrino particles is |$N_\nu =2^{3 \, l_\mathrm{min}}$| in 3D and the mass of a single sterile neutrino particle is m_ν = Ω_νρ_critb³/N_ν, with b being the box size of the simulation volume. We also consider different comoving box sizes of b200 ≡ 200 cMpc h⁻¹ and b400 ≡ 400 cMpc h⁻¹, leading to a maximum spatial resolution of 195.3 ckpc h⁻¹ (b200l10), 390.6 ckpc h⁻¹ (b200l9 and b400l10), 781.2 ckpc h⁻¹ (b200l8 and b400l9), and 1562.4 ckpc h⁻¹ (b400l8), with the number at the end of each simulation code name referring to its maximum refinement level. music needs to generate initial conditions at all the refinement levels which may be used by the main por simulation.⁸

2.5 Finding gravitationally bound structures

We find gravitationally bound structures using the amiga halo finder (ahf; Gill, Knebe & Gibson 2004; Knollmann & Knebe 2009), partly because it comes with a converter for ramses data, i.e. for converting grid cells into ‘cell particles’ located at the cell centre with the same mass.⁹ Like ramses, ahf is an AMR-based code that uses the adaptive refinement strategy to identify subhaloes at a modest computational cost.

ahf identifies galaxies depending on the density in a certain volume and afterwards neglects gravitationally unbound structures. However, the criterion for identifying these is not valid in MOND, so we have to use the raw data and skip the latter step. As there is no halo finder currently available that also has the option to use MOND gravity and building one is beyond the scope of this project, we stress that this introduces a bias towards less massive objects. However, the reader should keep in mind that, due to the greater self-gravity, Milgromian structures are generally more stable than Newtonian ones, so many of the identified low mass structures are expected to be gravitationally bound. Additionally, we restrict our analysis to structures that contain baryons because low- and intermediate-mass structures and substructures that only contain sterile neutrinos are most likely chance alignments. Our structures are identified considering both baryons and neutrinos and must have at least 10 cell particles and/or sterile neutrinos. Our analysis of the simulations is still at a rather preliminary stage due to the difficulties created by the use of a non-linear gravity law for which relatively little prior work is available, especially on cosmological scales.

To compare our results with observations, we consider the edge of a structure to be at the radius R₁₈₀, the radius inside which the mean mass density is |$180 \, \rho _\mathrm{crit}$| (equation (9)).

with the mass considering both baryons and sterile neutrinos. The halo finder calculates the mass and radius iteratively starting from an initial guess.

There is an additional problem when comparing observations with MOND simulations: for a given gravitational field, the corresponding mass profile depends on the gravity law. Since mass is not a direct observable, the mass deduced in the context of an assumed gravity law is called the dynamical mass in that theory. Most works in the literature quote the Newtonian dynamical mass profile M_N(r), which is the enclosed Newtonian dynamical mass within radius r. By definition, M_N(r) ≡ r²g(r)/G for a spherical system, with the gravitational field g(r) generally being what the observations actually constrain through kinematical information. In MOND, the Newtonian gravitational acceleration is enhanced by a factor ν, but this enhancement can also be achieved within Newtonian gravity if we rescale the enclosed mass M(r) by this factor (Angus et al. 2011):

We include equation (12) in the halo finder, so every M₁₈₀ value reported here is actually the Newtonian dynamical mass obtained from the total mass (gas + neutrinos) in the simulation data by applying the above boost (this issue is discussed further in section 4.3 of Asencio et al. 2021). This allows for an easier comparison with observations as the literature usually gives the Newtonian dynamical mass rather than the MOND one.

2.6 Simulation parameters

Every simulation done for this work starts at z_e = 199 or a_e = 0.005 because the typical gravitational field from inhomogeneities enters the MOND regime only when |$z_e \lesssim 50$| (see section 3.1.3 of Haslbauer et al. 2020), so our choice achieves a good compromise between limiting the effect of MOND at later times and limiting the impact of radiation, which is dominant only at much earlier times. We run six models: b200l10, b400l10, b200l9, b400l9, b200l8, and b400l8. The labels ‘10’, ‘9’, and ‘8’ stand for the maximum refinement level used in the simulation (see Section 2.4), while ‘200’ and ‘400’ show the length of the simulation box in cMpc h⁻¹. These models are summarized in Table 1.

All simulations were performed on in-house SPODYR-group machines in Bonn, which have 64 CPU cores, 128 threads, and 512 GB of RAM. The high-resolution run (model b200l10) took a computational time of 5 d with 64 threads, corresponding to 7680 CPU hours.

3 RESULTS

Our simulations do not contain a key feature of standard simulations, namely CDM particles creating deep potential wells into which baryons fall after cooling (e.g. White & Rees 1978; Frenk & White 2012). Galaxies are never the less generally expected to form earlier in MOND than in the SMoC (Sanders 1998, 2001) because of the strong long-range enhancement to the gravitational acceleration. Resolving the formation of small structures is highly dependent on the maximum refinement level and the box size. Due to the lack of power on small scales (see the right-hand panel of Fig. 2) and the influence of structures on each other through tides and the external field effect (EFE; see sections 2.4 and 3.3 of Banik & Zhao 2022), it is crucial to simulate large enough boxes in order to properly capture structure formation in the νHDM model. This makes it very taxing to do high resolution simulations, especially hydrodynamical ones.

In the highest resolution model (b200l10), the first structures can only be identified when z_e = 4. This late formation of structures seems to be in discordance with observations: there are many very massive galaxies at |$z_e \gtrsim 4$|⁠, a fact which has become clearer with deeper and improved observations (Steinhardt et al. 2016; Marrone et al. 2018; Wang et al. 2019; Neeleman et al. 2020; Rennehan et al. 2020; Tsukui & Iguchi 2021), including very recently with data from the JWST (e.g. Furtak et al. 2022; Harikane et al. 2022; Labbe et al. 2022; Naidu et al. 2022; Haslbauer et al. 2022a; Adams et al. 2023; Atek et al. 2023; Yan et al. 2023). We discuss this issue in more detail in Section 3.2.

Fig. 3 shows the structures which form in our models by the present epoch, with the top panels showing the baryons and the bottom panels showing the sterile neutrinos.¹⁰ Comparing the left- and right-hand panels of Fig. 3 illustrates how the box size affects the growth of structure throughout cosmic time: model b400l10 (right-hand panels) has significantly more structure than b200l10 (left-hand panels). This demonstrates the difficulty of achieving ‘Copernican convergence’, i.e. having a sufficiently large volume such that the statistical properties do not change if a larger volume is considered. Copernican convergence is harder to achieve in MOND because the gravity law is non-linear, so the effect of a structure in the simulation on other structures declines much more slowly with distance. Moreover, structures on Mpc scales form in MOND in a top-down manner because there is very little power on small scales initially. Our models are not fully Copernican converged due to computational limitations.

3.1 The mass function

As our first quantitative analysis, we discuss the mass distribution of the formed structures and how this evolves with time. Fig. 4 shows the mass distribution of structures in all six models at different times (left-hand panel: b200 models, right-hand panel: b400 models, top: l_max = 8, middle: l_max = 9, bottom: l_max = 10). The mass of the most massive structure (top of each error bar) and the least massive structure (bottom of each error bar) grow over time, as is also the case for the mean and median mass (the filled points and black open points, respectively). The impact of resolution is apparent in that the first structures are identified earlier in the higher resolution models. The mass of the most massive structure increases when we consider a larger box, which is expected because larger structures can arise in a larger volume. These preliminary results are thus by no means numerically converged. However, they do clearly illustrate the impact of box size and resolution. For instance, the least massive structures in the b400 models (right-hand panels) are more massive than the least massive structures in the b200 models with the same l_max (left-hand panels) because the spatial resolution is higher in the latter case.

Fig. 5 shows the present-day (z_e = 0) differential mass function (baryons + neutrinos) of structures in all our simulations, with the panels organized in the same way as in Fig. 4. We place the structures into mass bins of width 0.5 dex and divide the number counts by the product of the bin width and the total comoving volume under consideration. The distributions are then fit by a power law:

where A is the normalization and α is the power-law index. The fitted values for these quantities are shown in Table 2 for every model.

Fig. 5 illustrates the mass distribution of the structures at the end of the simulations. When comparing simulations of the same box size (in the same column), the effect of resolution is clearly evident in that the lowest mass bin gets more populated at higher resolution. Additionally, the shape of the distribution seems to be changing from a simple power law closer to a Schechter function (Schechter 1976) with higher resolution, which is the observationally expected shape − this is especially evident for the b200 models. However, increasing the box size shifts the distribution to more massive structures. Not only are the most massive structures found in the larger box (as is evident from Fig. 4), but also the number density of such structures is higher in a larger box.

To estimate the cosmic variance in our simulations, we divide the whole b200 simulation cube into eight equally sized cubes (octants) with sides of length b = 100 cMpc h⁻¹, which is half that of the original simulation. Fig. 6 shows the mass function in each octant for model b200l10 at z_e = 1.0 (left-hand panel) and today (right-hand panel). The aforementioned existence of pronounced voids and overdense regions is also evident here in that there is more than one dex difference between the mass function in the least and most dense octants. Additionally, the figure shows that the mass of the most massive structure depends on the environment: the least dense octant forms structures with masses up to |$M_{180} = 10^{15.0} \, \mathrm{M}_\odot$|⁠, while the most dense sub-volumes contain structures with |$10^{17.0} \, \mathrm{M}_\odot \lt M_{180} \lt 10^{17.5} \, \mathrm{M}_\odot$| at z_e = 0.0. Qualitatively the same phenomenon is evident at z_e = 1.0 but shifted to lower masses, reaching |$10^{14.0} \, \mathrm{M}_\odot \lt M_{180} \lt 10^{14.5} \, \mathrm{M}_\odot$| for the least dense octant, while the two densest octants reach masses |$M_{180} \gt 10^{16} \, \mathrm{M}_\odot$|⁠.

The νHDM model seems to explain the existence of very massive galaxy clusters at high redshift like the El Gordo cluster (Menanteau et al. 2012): at z_e = 1, the νHDM mass function extends up to |$M_{180} \approx 10^{16} \, \mathrm{M}_\odot$|⁠. El Gordo has a total mass of about |$2 \times 10^{15} \, \mathrm{M}_\odot$| (Kim et al. 2021) and is composed of two interacting subclusters with a similar mass and an infall velocity of v_infall = 2500 km s⁻¹ at a redshift of z_e = 0.87 (Zhang, Yu & Lu 2015). Since El Gordo is observed after pericentre, the two progenitor clusters need to be spatially close at z_e = 1 and have a combined mass of |$\approx 2 \times 10^{15} \, \mathrm{M}_\odot$|⁠. Given also the limited sky area in which El Gordo was discovered, this is plausible only if the mass function at that time reaches somewhat higher masses − which it does in some octants (indicated by the vertical dashed line in the left-hand panel of Fig. 6). Our results are compatible with the work of Katz et al. (2013), who found El Gordo analogues in their νHDM simulations with box sizes of ≈750 cMpc (see their section 5.1). By comparing their result to the effective volume of the discovery survey, it was shown that νHDM produces approximately the correct frequency of El Gordo analogues (see section 4.3 of Asencio et al. 2021).

After analysing in detail the mass distribution of the structures in our simulations, we use Fig. 7 to compare them with the observed mass distribution of galaxy clusters. Wang et al. (2022b) obtained weak lensing masses for the clusters from the catalogue of Yang et al. (2021), which was generated by the halo-based group finder originally developed in Yang et al. (2005). The observations are based on the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys (Dey et al. 2019) and corrected for completeness. The cumulative mass function (CMF) is calculated by counting the number of clusters beyond certain mass thresholds and dividing by the respective volume under consideration for five different redshift slices with a width of 0.1, thus covering the range z_e = 0.1−0.6.

For an accurate comparison, we apply the same method to our simulations: we count the number of structures above the mass thresholds used in Wang et al. (2022b) and divide by the volume of the respective simulation for outputs that fall into the redshift range, i.e. we compare the observed redhsift bins only with outputs that are in the same redshift range. Fig. 7 shows the CMF for all five redshift bins and all six simulations. Every column depicts the same redshift bin for different models (top to bottom: b200l8, b200l9, b200l10, b400l8, b400l9, and b400l10); while every row shows the redshift evolution for one model, with redshift increasing from left to right.

The most obvious finding is that all models seemingly overproduce high mass structures in all redhsift bins, which we turn to momentarily. The simulated CMFs are also significantly flatter than observed. The two agree approximately at intermediate masses, but there is an obvious impact of resolution on the νHDM CMFs at the low mass end. The l8 models show a nearly horizontal CMF at lower masses, indicating that there are no structures at or below these masses (see also Fig. 5). As anticipated, the number density at any fixed mass increases with decreasing redshift for both observations and our models, which can be attributed to accretion and mergers (cf. Fig. 4). Taking the discrepancy between observations and our models at high masses at face value puts significant doubt on the theory. A similar conclusion was drawn by Angus et al. (2013) and Katz et al. (2013), who also did cosmological simulations in the νHDM framework. However, we argue below that their conclusions were premature.

Shortly after the publication of Angus et al. (2013) and Katz et al. (2013), it became apparent that for redshifts up to z_e = 0.07 (corresponding to a distance of ≈300 Mpc), a mismatch in the high mass range is to be expected given our location in a large underdensity known as the KBC void (Keenan et al. 2013). Their study and the more recent study of Wong et al. (2022) cover 90 per cent of the sky, indicating that the observational results on the KBC void are not due to incomplete coverage. They also probe a large portion of the galaxy luminosity function. It is therefore likely that the mass function of observed clusters in the Local Volume is not representative of an average region of the same size. Interestingly, Haslbauer et al. (2020) examined the KBC void and concluded that it naturally explains why the locally determined Hubble–Lemaitre constant H₀ exceeds that determined from CMB anisotropies (the Hubble tension is reviewed in Di Valentino 2021). Such a void was initially deemed to be incompatible with the supernova distance-redshift relation (Kenworthy, Scolnic & Riess 2019), but subsequent analyses have found that supernova data alone actually reveal a mild preference for a local void (Kazantzidis & Perivolaropoulos 2020; Luković, Haridasu & Vittorio 2020). Indeed, the results of Castello, Högås & Mörtsell (2022) indicate that the preferred void parameters are remarkably close to that of the KBC void, highlighting a consistency between results obtained from galaxy number counts and from supernovae. Further, Dainotti et al. (2021, 2022) and Jia, Hu & Wang (2022) investigate the redshift evolution of H₀ and find clear evidence for a smooth decrease from the high local value down to the lower value inferred from the CMB (Planck Collaboration VI 2020) out to z_e = 2. Schiavone, Giovanni & Bombacigno (2023) interpret these results via an effective Hubble constant that evolves with the redshift in an f(R) modified gravity scenario, thus indicating an intrinsically different local value of H₀. This is very interesting given how Haslbauer et al. (2020) argued that a local underdensity can change the measured H₀: the finding of a smooth transition might indicate that we are situated in a much larger underdensity, which in turn would lead to the local measurements not representing an average region of the Universe. Indeed, fig. 4 of Haslbauer et al. (2020) shows that their inference on the void radius has an extended tail that reaches very large values. Further support for the Universe being more inhomogeneous than assumed by the SMoC comes from Migkas et al. (2021), who analysed 570 clusters in the local Universe and built 10 different cluster scaling relations out of these to test the hypothesis of isotropy. They found a dipole-like anisotropy in the local value of H₀ with more than 5σ significance. This study underlines the question raised before: Is the local Universe representative of an average region of the Universe?

One way to test cosmological models despite a possible negative answer is to consider observations of the most massive clusters up to z_e = 1. El Gordo is not only more massive than the observed clusters closer to us shown in Wang, Hammer & Yang (2022a), but it is also further in the past. We would expect clusters today to reach an even higher mass. Asencio et al. (2021) analysed El Gordo in detail by looking at realistic hydrodynamical simulations of the pre-merger configuration and comparing the viable configurations to large-scale ΛCDM simulations. Those authors found that ‘detecting one pair with its mass and redshift rules out ΛCDM cosmology at 6.16σ’.

The two sub-clusters that merged into El Gordo need to have masses of about |$10^{15}\, \mathrm{M}_\odot$| at a redshift close to z_e = 1.0. From observations out to z_e = 0.6, it seems unlikely to find even one cluster with such a high mass, but to have two of them in close proximity at a larger redshift strongly indicates that there should be more structure on large scales than is generally expected and observed locally. Encouragingly, Asencio et al. (2021) argue in detail in their section 4.3 that the previous νHDM simulations done by Katz et al. (2013) show the correct frequency of El Gordo analogues. As a result, the apparent overproduction of massive structures (Angus et al. 2013) may well be incorrect as their results were based on comparison with low-redshift data sets. We therefore argue that the νHDM theory cannot be ruled out by the discrepancies at high masses evident in Fig. 7. In fact, the νHDM simulations reported here confirm that bound structures with an El Gordo-like mass readily appear at z_e = 1 (Fig. 4) and also show that these mass scales are only reached in dense regions of the simulation box at that redshift (see Fig. 6). Given that the mass of El Gordo is now known to about 10 per cent accuracy (Kim et al. 2021), any realistic cosmological model must be able to reproduce its properties, especially given the other potentially problematic clusters and cluster pairs discussed in the introduction to Asencio et al. (2021).

3.2 Time evolution of bound structures

One of the major processes affecting the galaxy population in the SMoC is mergers between galaxies, which are driven by dynamical friction between extended CDM haloes (Kroupa 2015). Since these are absent in MOND, the baryonic parts of galaxies would need to overlap for them to merge − they would simply fly past each other in the absence of dissipation (Renaud et al. 2016; Bílek et al. 2018; Banik et al. 2022). The collision cross-section between galaxies is thus much smaller in MOND, which may explain the high observed fraction of thin disc galaxies (Peebles 2020; Haslbauer et al. 2022b).

To explore the time evolution of the population of structures in MOND, we use Fig. 8. Its left-hand panel shows how the total number of structures (N_struc) evolves over time for every model, while the right-hand panel shows the fractional reduction in the number of structures between the epoch when this was maximal and the present epoch. In other words, it shows the merger fraction

where N_max is the maximum number of structures (attained at redshift z_max) and N₀ is the present number of structures in the simulation volume. This is only an estimated merger fraction because it is possible that new structures form at later epochs in these simulations. However, we expect this to be rather rare because the mass of the least massive structure rises with time (see Fig. 4).

The evolution of N_struc is qualitatively similar for all models. The total number of structures rises steeply after the onset of structure formation at about z_e = 4 until it reaches a maximum at z_e = 1. Only the highest and the lowest resolution model (b200l10, brown and b400l8, light blue) peak earlier and later, respectively. Afterwards, N_struc declines modestly to significantly, depending on the resolution: higher resolution models show a steeper decline. This indicates that mergers play a role in hydrodynamical MONDian models that lack star formation. Important to note is that the observed star formation rate density also peaks at a redshift of z_e ≈ 2 (Madau & Dickinson 2014), indicating that the number of structures may also peak at this time, which is compatible with our highest resolution model (b200l10). Interestingly, the models with the same box size end with nearly the same number of structures for l_max = 9 and 10, which seems to indicate that the last output (z_e = 0) is much closer to numerical convergence compared to the early universe in these simulations. Additionally, as expected, the simulations with an initially larger box (the b400 models shown in green, orange, and light blue) end up with a larger number of structures.

The right-hand panel of Fig. 8 shows the fraction of mergers of structures in each simulation (equation 14). There is a clear dependence on resolution, showing that increasing l_max and/or decreasing the box size increases f_merger. Also, the maximum of the number of structures mentioned above is more easily visible here, with most of the models peaking at z_e = 1. As shown in Fig. 4, the models computed here are not able to reach single galaxy scales at the onset of structure formation (here z_e ≈ 4). Sanders (1998) already estimated the time at which spheres of different masses should collapse in an idealized scenario, i.e. without the EFE and in spherical symmetry. These first analytic calculations estimate the onset of cluster formation at about z_e = 3, which is consistent with the onset of structure formation in the models shown here. Therefore, higher resolution simulations are needed that at least probe mass scales down to |$10^{10} \, \mathrm{M}_\odot$| in order to conclusively address if structures form early enough in the νHDM model compared to the most recent observations. However, it should be noted that the analytical work of Sanders (1998) was not based on the νHDM framework and in particular has a different background cosmology, i.e. an open universe with Ω_b = 0.02 and no dark matter or dark energy, potentially causing difficulties fitting the CMB and other cosmological data sets collected over the last quarter century. We note that their adopted cosmology gives almost twice as much time for structure formation by z_e = 10 because the different expansion history makes the universe older at that redshift. Moreover, increasing the resolution of the νHDM models causes only a modest evolution towards higher redshift for the onset of structure formation, so it seems unrealistic to form galaxy mass structures at z_e ≥ 10. This may already indicate a failure of the model.

While these first hydrodynamical cosmological MOND results suggest that the population of structures builds up until z_e ≈ 1−2 and then decreases again due to mergers, the results should be considered indicative rather than conclusive. This is because the merger rate is undoubtedly enhanced by virtue of all baryons being purely gaseous and thus able to dissipate energy rather efficiently during interactions. Future simulations will need to include star formation to obtain a more realistic assessment of the late-time merger rate and also a higher resolution to analyse the behaviour of individual galaxies.

Another aspect of the simulations which changes with time is the total baryonic mass inside structures. Its time evolution is shown in Fig. 9, where M_{bar, struc} is the total baryonic mass inside all structures and M_{bar, box} is the baryonic mass inside the box, which does not change with time. A distinct feature at z_e = 2 can also be seen here: at later times, the baryonic mass fraction increases significantly in the same way for all models up to a value of 0.3−0.5. This behaviour of all models indicates that the increase for z_e < 2 is driven by the accretion of baryons on to large structures rather than the formation of new structures or mergers thereof, which is further supported by the fact that f_merger is significantly lower for the low resolution models where small structures cannot form at high z_e. These results indicate that massive enough regions have collapsed and formed deep potential wells at z_e ≈ 2, such that significant amounts of baryonic mass move through filaments of the cosmic web towards these structures, where the filaments intersect and clusters form (this is compatible with the theoretical findings of Sanders 1998). What is more, a modest trend with box size is evident for the present epoch baryonic mass fraction in identifiable structures, which is closer to 0.4 for the b200 models and closer to 0.3 for the b400 models. Additionally, there is also a minor decrease with increasing resolution regardless of the box size.

3.3 The sterile neutrino mass fraction

These first ever hydrodynamical cosmological MOND simulations with both baryons (gas) and live sterile neutrino particles enable us to explore whether MOND might be consistent with the observed ratio of dynamical mass to baryonic mass, M_M/M_b, in massive galaxy clusters, where the ratio appears to be very close to the cosmic fraction (see section 7.1 of Banik & Zhao 2022, and references therein). Since MOND was originally proposed to explain observations with only baryonic mass, we compute the mass fraction of sterile neutrinos, f_ν ≡ M_ν/M_tot. Observationally, M_tot corresponds to M_M, the MOND dynamical mass. In our simulations, M_tot ≡ M_b + M_ν. This can readily be connected to the observed ratio of dynamical mass to baryonic mass: f_ν ≡ 1 − M_b/M_M.

Fig. 10 compares f_ν for every model (colours and panel placement are the same as in Fig. 4) with the findings from Pointecouteau & Silk (2005, hereafter PS05), who analysed M_M/M_b in a MOND context for 10 observed clusters out to a redshift of z_e = 0.15 in the mass range |$1.2\times 10^{14} \, \mathrm{M}_\odot$| to |$1.2\times 10^{15} \, \mathrm{M}_\odot$|⁠. They found that M_M/M_b = 4.94 ± 0.50 at distances from the cluster center of |$r\approx 0.5 \, R_\mathrm{vir}$|⁠, in contrast to the SMoC asymptotic (⁠|$r \ge 0.3 \, R_\mathrm{vir}$|⁠) value of |$7.7^{+1.4}_{-1.1}$| (in the SMoC most of the mass is CDM, with neutrinos playing a negligible role). This leads to an observational value for f_ν of 1 − M_b/M_M = 0.80 ± 0.02, which we show in Fig. 10. We bin our data in the same mass bins used in Fig. 5 (0.5 dex bin width) and show the mean value with its standard deviation.

It is quite striking that the models come out to be in good agreement with the observational data for f_ν at the high mass end, especially for the higher resolution models. The slightly higher asymptotic value of all models compared to the observed value is expected because f_ν for the simulations is calculated at R₁₈₀ but PS05 evaluated their findings at |$r\approx 0.5 \, R_\mathrm{200}$|⁠, which is closer to the cluster center. As a result, the simulated value should be closer to the cosmic mean value of Ω_ν/Ω_m ≈ 0.84, but since the dissipative baryons are expected to be more centrally concentrated than the sterile neutrinos, the neutrino mass fraction should be lower within a smaller radius. The highest resolution model (b200l10; panel e) shows a smooth decline of f_ν from the cosmic mean value for massive clusters down to a nearly vanishing sterile neutrino component for the least massive structures identified (see Fig. 5). This nicely shows the transition between clusters and galaxies, where the sterile neutrinos do not cluster by construction in the νHDM model. Indeed, our results show that the Local Group, with a baryonic mass of |$\approx 10^{11} \, \mathrm{M}_\odot$|⁠, does not have a significant sterile neutrino component in the νHDM paradigm, in agreement with the past dynamical history of the Milky Way-M31 galaxy pair (the timing argument; Banik, O’Ryan & Zhao 2018; Benisty & Guendelman 2020; McLeod & Lahav 2020; Banik et al. 2022). More generally, the almost purely baryonic nature of low-mass structures is very much in line with the aim of νHDM to preserve the empirical successes of MOND in galaxies. Important to note here is that the bins for the lowest masses may be populated by only a few structures, leading to a much larger statistical uncertainty than is shown with the mean value and its standard deviation. This may explain why f_ν sometimes unexpectedly decreases with increasing mass in the lowest mass bins, but this does not alter the general picture of negligible f_ν below |$10^{14} \, \mathrm{M}_\odot$| and f_ν ≈ 0.84 above |$10^{14} \, \mathrm{M}_\odot$|⁠. This unexpected behaviour noted here has to be re-evaluated with higher resolution simulations, investigating if this effect depends solely on resolution/mass of the sterile neutrino particles or has an entirely different cause.

3.4 Feasibility of the observed Local Group peculiar velocity

In addition to the mostly mass-related analysis presented here, we show in Fig. 11 the distribution and time evolution of the peculiar velocities of all structures. We use this to quantify the likelihood of finding a structure with a peculiar velocity as low as that of the Local Group (v_pec ≤ 630 km s⁻¹; Kogut et al. 1993).

Fig. 11 illustrates the peculiar velocity distribution of all structures at z_e = 0.0 and z_e = 1.0. The number of structures in each velocity bin of width 500 km s⁻¹ is normalized by the total number of structures, N_{tot, struc}, at the relevant redshift, with error bars showing the Poisson uncertainty. The black distributions correspond to z_e = 1.0 and the coloured distributions show the peculiar velocities today. The distributions in different simulations at z_e = 1.0 are very similar, showing a maximum between 1000 and 2000 km s⁻¹ followed by a smooth decline until the highest peculiar velocities are reached at about v_pec = 6000 km s⁻¹ for the b200 models and between 7000 and 8000 km s⁻¹ for the b400 models. By z_e = 0.0, the distributions have broadened significantly, shifting the mode and the maximum reached peculiar velocity to larger values. Also, the b200 models show more features in the extended high velocity tail of the distribution, which may be attributed to a smaller number of structures and higher density variations throughout the simulation box (see Fig. 3). Overall, the peculiar velocities seem to be nearly thermal, i.e. the distributions are comparable to Maxwell–Boltzmann distributions (Maxwell 1877), as is evident from the asymmetry due to the extended high velocity tail.

Earlier studies of particle-only νHDM simulations (Katz et al. 2013) show comparable behaviour for the evolution of the median peculiar velocity (see their fig. 4). Recent observations of large bulk flow velocities seem to indicate some tension with the SMoC (Watkins et al. 2023), which may be related to larger than expected peculiar velocities and therefore could be more natural in the νHDM framework (see also Vittorio, Juszkiewicz & Davis 1986). Additionally, compatible with the findings from the CMF (see Fig. 7), the peculiar velocities of the different models exceed those of local structures. In particular, the Local Group has a peculiar velocity of v_pec = 630 km s⁻¹ (Kogut et al. 1993). In the following, we quantify the fraction of structures in each model which have a smaller peculiar velocity. We express this fraction as a tension in terms of the equivalent number of standard deviations (results are presented in Table 3).

Haslbauer et al. (2020) point out in their section 4.2 that the peculiar velocity of the Local Group is much slower than typically expected in νHDM such that this observable causes the highest tension between their semi-analytical νHDM model and the observational constraints. The tension is still only 2.34σ, which is hardly impossible. We find a tension of 2.5σ−2.9σ, with the b200 models having a slightly smaller value of χ than the corresponding b400 models. The models computed for this work show a similar tension to that estimated by Haslbauer et al. (2020) based on the fraction of their simulated void where the peculiar velocity in the CMB frame is smaller than that observed for the Local Group. Therefore, we conclude that it is unlikely but not impossible to find structures with a comparably low peculiar velocity as the Local Group, which again suggests that the Local Volume is not an average volume of a νHDM universe − and might also not be for the observed Universe (see Section 3.1).

4 CONCLUSIONS AND FUTURE WORK

We report the first hydrodynamical MOND simulations of structure formation in a cosmological context. The calculations are performed in the νHDM cosmological framework proposed by Angus (2009) in which the CDM component in ΛCDM is replaced by a collisionless component in the form of sterile neutrinos with a mass of m_νs = 11 eV c⁻². This component is motivated from neutrino physics rather than cosmology but is equivalent to a hot ‘dark matter’ component. This conserves the MONDian behaviour on galactic scales but is compatible with galaxy cluster dynamics (Angus et al. 2010) and the CMB (Angus & Diaferio 2011) due to the strong gravitational fields in the early universe and a standard evolution of the cosmic scale factor a_e (for further details, see Haslbauer et al. 2020). We generate the initial conditions at z_e = 199 for this cosmological model by calculating the power spectrum with camb while including the necessary changes to only the namelist to replace the CDM component with 11 eV c⁻² sterile neutrinos. We then sample the power spectrum using the music program to generate the initial conditions of the live sterile neutrino particles and the gas distribution. As with other cosmological simulations (e.g. Candlish 2016), we solve the Poisson equation in supercomoving coordinates (equation (8); see also Appendix  A and Martel & Shapiro 1998; Teyssier 2002). Applying this to MOND implicitly assumes that the MOND Poisson equation should be applied only to the difference between the density and the cosmic mean value, which we argued to be reasonable in Section 2.3 (see also section 3.1.4 of Haslbauer et al. 2020). Recent analytic work in covariant MOND theories further supports this approach (Thomas et al. 2023).

Six models were computed in this work, exploring the impact of three different resolution settings and two box sizes (Section 2.4). We concentrate on small to intermediate sizes (200−400 cMpc h⁻¹) in order to resolve the smallest structures possible, while initializing the simulations with enough power to be able to form any structures at all. As the right-hand panel of Fig. 2 depicts, scales below 100 cMpc h⁻¹ initially have very little power in the matter power spectrum, so starting a simulation with an initial box size in this regime leads to vanishingly small density fluctuations. This would cause structures to arise from numerical noise only, if indeed any form at all. Structure formation must be top-down in the νHDM framework, so large box sizes are essential.

The first structures begin to form at z_e ≈ 4 for the higher resolution simulations. This result appears to be in tension with high redshift galaxies observed with the recently launched JWST, which also seem to be in some tension with ΛCDM (Haslbauer et al. 2022a). However, previous analytic estimates of structure formation in MOND indicate that galaxy scale structures (up to |$10^{11} \, \mathrm{M}_\odot$|⁠) should have formed by z_e = 5−10, while cluster-sized structures (⁠|$10^{14} \, \mathrm{M}_\odot$|⁠) should have reached maximum expansion and begun to recollapse at z_e = 3 (Sanders 1998). The latter is comparable to the evolution of the models shown here, but important to note is that the analytic estimates are calculated for an idealized spherical scenario without sterile neutrinos and dark energy, representing a rather different framework that may be in tension with the observed CMB. Simulations that adequately resolve galaxy mass scales (rather than just reaching this scale as we do here) are needed to conclusively test the onset of structure formation in the νHDM framework. We expect a higher resolution simulation to identify the first galaxies at earlier epochs than the structures shown here, with star formation in the first collapsing gas clouds occurring even earlier (see also Wittenburg et al. 2020; Eappen et al. 2022).

The most massive galaxy clusters (up to |$M_{180}=10^{17} \, \mathrm{M}_\odot$|⁠) form in the models with the largest box size (b = 400 cMpc h⁻¹) regardless of resolution (see Fig. 4). With improved resolution, it may be worthwhile to search for analogues to observed massive galaxies at high redshift (e.g. Furtak et al. 2022; Harikane et al. 2022; Labbe et al. 2022; Naidu et al. 2022; Adams et al. 2023; Atek et al. 2023; Yan et al. 2023). If these structures appear in improved simulations with a larger box size, it will become possible to estimate the number density of the most massive galaxy cluster collisions like El Gordo (Asencio et al. 2021) in a statistically meaningful manner.

We also show the differential mass function at z_e = 0 (Fig. 5) and compare the simulated CMF with weak lensing data (Wang et al. 2022b) up to z_e = 0.6 (Fig. 7). There is reasonable agreement in the intermediate and low mass regimes, but the slope of the simulated distribution is too shallow, leading to an overabundance of massive structures at the high mass end compared to the weak lensing catalogue. The overabundance is similar to that mentioned in Angus et al. (2013) and Katz et al. (2013), but we argue that, especially for the highest masses, the relatively local observations might not be representative of an average volume of the Universe. In particular, we are situated in a supervoid with a radius of r ≈ 300 Mpc (the KBC void; Keenan et al. 2013; Haslbauer et al. 2020), which inhibits the existence of very massive clusters as the most massive structures can only form in the most dense regions of the universe. Additionally, observations out to z_e ≈ 2 (Dainotti et al. 2021, 2022; Jia et al. 2022) show a smooth decrease of the inferred H₀ from the high local value (Riess et al. 2019) towards the smaller value inferred from the CMB (Planck Collaboration VI 2020). Assuming that the tension between the two values of H₀ vanishes when considering a large underdensity in MOND (Haslbauer et al. 2020), one interpretation of a smooth transition between the two values might be that the underdensity around us is larger than shown by the KBC void (the inference on the void radius in their fig. 4 does indeed show an extended tail towards large values). This would lead to the universe not being as homogeneous as predicted by the SMoC. There are indeed observed structures like El Gordo which are at relatively high redshift (z_e = 0.87) and are extremely massive (⁠|$M_{200}=2\times 10^{15} \, \mathrm{M}_\odot$|⁠; Kim et al. 2021), with El Gordo alone ruling out the ΛCDM cosmology at >6σ confidence (Asencio et al. 2021). We therefore argue that the νHDM model cannot yet be ruled out by the mismatch between observations and simulations at the high mass end of the cluster mass function.

To improve our understanding of this issue, a more detailed analysis would be required to check if the galaxy or galaxy cluster mass function for an observer in a low-density region matches the observations and if large volume (Gpc scale) simulations can reproduce El Gordo analogues with the observed frequency. As a first step towards this, we divide the b200l10 model into eight equally sized octants and compare their mass functions at z_e = 0 and z_e = 1. Fig. 6 shows that there is more than one dex difference in density between the sub-volumes and also that the most massive clusters only form in the densest regions. The figure also shows that our models form clusters of the same mass scale as the progenitors of El Gordo at z_e = 1, indicating that νHDM could be consistent with the observations of El Gordo (as argued in section 4.3 of Asencio et al. 2021).

We also examine the time evolution of the structures by showing how the number of structures and the baryonic mass inside all structures evolve with time (Figs 8 and 9, respectively). The former yields the first estimates for the prevalence of mergers in a cosmological MOND framework, while also showing that this seems to be resolution dependent. The latter, in combination with the former, suggests that the accretion of gas on to structures dominates the growth of mass in structures for z_e < 2 in these models.

Additionally, we are for the first time able to show the mass ratio of baryons to sterile neutrinos in bound structures in the form of the sterile neutrino mass fraction f_ν, which we compare to the observed ratio of MOND dynamical mass to baryonic mass (PS05). Fig. 10 shows that the higher resolution simulations are consistent with the observations in the observed mass range. In particular, f_ν becomes asymptotically flat at the cosmic value of Ω_ν/Ω_m ≈ 0.84 for structures more massive than |$10^{15} \, \mathrm{M}_\odot$|⁠, a result that is only weakly dependent on resolution and box size. The highest resolution simulation (panel e) confirms the expected νHDM behaviour that the fraction of mass in sterile neutrinos decreases from cluster mass scales down to galaxy mass scales, where it vanishes. This is required in any MOND cosmology because galaxies can be explained purely by their baryonic mass content in MOND (Famaey & McGaugh 2012; McGaugh 2020; Brouwer et al. 2021; Banik & Zhao 2022, and references therein).

Furthermore, we show the distribution of peculiar velocities of all structures and calculate the likelihood of finding structures with a peculiar velocity as small as that of the Local Group. We find this to cause between 2.5σ and 2.9σ tension depending on the model, with the b200 models having a slightly smaller tension compared to their b400 counterparts (see Table 3). Our result is compatible with the findings of Haslbauer et al. (2020), according to which a 2.34σ tension exists between observations and their semi-analytical νHDM model in this respect. This shows that a peculiar velocity as low as the observed 630 km s ⁻¹ is unlikely but not impossible in νHDM. Therefore, the Local Volume is rare in a νHDM universe and perhaps also in the real Universe, as discussed further in Section 3.1. We also note that the large observed bulk flow out to ≈200 Mpc is in some tension with ΛCDM (Watkins et al. 2023) and may be a sign of the typically much larger peculiar velocities in νHDM (see our Fig. 11 and also fig. 4 of Katz et al. 2013).

In this work, we present the first hydrodynamical cosmological MOND simulations. This is an important step towards a more detailed investigation of structure formation in a MOND cosmology. Further simulations similar to Katz et al. (2013) are needed to analyse the behaviour of this model on larger scales, where it appears to provide a promising explanation for the KBC void and Hubble tension (Haslbauer et al. 2020) and the massive high-redshift interacting galaxy cluster known as El Gordo (Asencio et al. 2021), phenomena which contradict ΛCDM at high significance. Moreover, the scales considered here need to be simulated at higher resolution with star formation and feedback prescriptions to investigate how single galaxies form and evolve in this model, which will be the task for future works. We are also conducting larger volume collisionless simulations to better explore the behaviour of our model on the largest scales (Russell et al., in preparation). We want to stress that improved simulations will be essential to check if the shallow mass function and the late onset of structure formation prevail in higher resolution simulations in the νHDM framework, which would falsify this model conclusively.

ACKNOWLEDGEMENTS

IB is supported by the Science and Technology Facilities Council grant ST/V000861/1. He acknowledges support from a ‘Pathways to Research’ fellowship from the University of Bonn in 2021 after an Alexander von Humboldt Foundation post-doctoral research fellowship (2018–2020). We thank the Deutscher Akademischer Austauschdienst-Eastern European exchange program at the University of Bonn for supporting the Bonn–Prague exchange and acknowledge financial support through the transdisciplinary research area TRA-Matter at the University of Bonn. The authors thank Benoit Famaey and Alfie Russell for useful discussions. They also thank the anonymous referee for helpful comments.

DATA AVAILABILITY

The algorithms used to prepare and run por simulations in a cosmological context and to extract their results into human-readable form are publicly available.¹¹ Nagesh et al. (2021) provides a user guide describing the operation of these algorithms.

Footnotes

References

APPENDIX A: THE POISSON EQUATION IN SUPERCOMOVING COORDINATES

Following the derivation of Martel & Shapiro (1998), we start with their equation (14), without the addition of the non-clumping X-component that they use:

where Φ is the Eulerian gravitational potential (equation (23) in Martel & Shapiro 1998),

with |$\bar{\rho }$| being the mean matter density, i.e. the left part of the right-hand side represents the background potential and ϕ the peculiar gravitational potential. We want ϕ as we need to apply the Jeans Swindle (Binney & Tremaine 1987, 2008), so we differentiate equation (A2) two times and additionally use equation (A1) to get the Poisson equation for ϕ:

This is the Poisson equation for cosmological calculations, which only applies to density contrasts rather than the total density field.

Now, we want to transform equation (A4) into supercomoving coordinates as these conserve the positions, the density, and the thermodynamic variables of an ideal gas in a cosmologically expanding volume in the absence of structure. To achieve this, several variables need to be transformed. We use the transformation from Martel & Shapiro (1998):¹²

where |$\Omega _{m,0}\equiv \bar{\rho }/\rho _{crit}$|⁠, the critical density |$\rho _{crit}\equiv 3H_0^2/8\mathrm{\pi } G$|⁠, ρ is the density, t is the time, ϕ is the gravitational potential, and |$\boldsymbol {r}$| is the spatial position. The tilde symbol denotes supercomoving coordinates. Every parameter with a star subscript is a necessary fixed parameter for the transformation, with r_⋆ being the only independent parameter − this is usually taken to be the simulation box size. Note that we use the t_⋆ definition from Teyssier (2002). We also need to convert the derivatives relative to r at fixed t to derivatives relative to |$\tilde {r}$| at fixed |$\tilde {t}$|⁠,

where |$\tilde {\nabla }$| is the gradient relative to |$\tilde {r}$| and f is an arbitrary function, which was not transformed here. The reader should keep in mind that if f depends on parameters that are transformed, then f also needs to be transformed accordingly.

We can start inserting the transformation now:

For QUMOND, ρ_tot is effectively ρ + ρ_P, with ρ_P in physical units:

Since |$\tilde {\rho }_{\it P}=a^3\rho _{\it P}/\rho _\star$|⁠, we get ρ_PDM in supercomoving units using again equations (A7), (A9), (A10)–(A12):

Important to note is that the argument of |$\tilde {\nu }$|⁠, i.e. |∇ϕ_N|/a₀, is dimensionless and therefore not affected by this coordinate transformation.

Finally, the Poisson equation for QUMOND in supercomoving coordinates is: