{'fraction_non_alphanumeric': 0.03417301308338215, 'fraction_numerical': 0.016559265768404607, 'mean_word_length': 4.640088105726872, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modified Baryonic Dynamics: two-component cosmological simulations with light sterile neutrinos

In this article we continue to test cosmological models centred on Modified Newtonian Dynamics (MOND) with light sterile neutrinos, which could in principle be a way to solve the fine-tuning problems of the standard model on galaxy scales while preserving successful predictions on larger scales. Due to previous failures of the simple MOND cosmological model, here we test a speculative model where the modified gravitational field is produced only by the baryons and the sterile neutrinos produce a purely Newtonian field (hence Modified Baryonic Dynamics). We use two component cosmological simulations to separate the baryonic N-body particles from the sterile neutrino ones. The premise is to attenuate the over-production of massive galaxy cluster halos which were prevalent in the original MOND plus light sterile neutrinos scenario. Theoretical issues with such a formulation notwithstanding, the Modified Baryonic Dynamics model fails to produce the correct amplitude for the galaxy cluster mass function for any reasonable value of the primordial power spectrum normalisation.


Introduction
Modified Newtonian Dynamics (MOND; [1], and see [2] for a recent review) is a modification to gravity in the ultra weak-field regime. It has a remarkable ability to reproduce the dynamics of galaxies of all types, shapes and sizes ( [3,4,5,6]) -with only a few notable caveats ( [7,8]). On these scales, the standard cosmological model still struggles to reproduce the observed galaxy regularities and has a number of other problems (see e.g. [9,10]. On large scales, like clusters of galaxies, the MOND approach does not work ( [11,12,13,14]). There is a missing mass problem for MOND in clusters of galaxies. There is also no clear way to describe cosmological phenomena like the anisotropies in the cosmic microwave background (CMB; [15,16,17,18], but see [19]).
Approaching the missing mass problem from the other direction, cosmologists have found a satisfactory representation of the Universe at large scales using standard gravity (i.e. General Relativity). This comes through the combination of cold dark matter and dark energy (ΛCDM), with baryons making up a mere 5%. Invoking this combination allows the ΛCDM (or concordance) model to successfully match the acoustic peaks in the CMB ( [20]). Furthermore, the observed distribution of matter on large scales ( [21,22,23]) is well reproduced by N-body simulations of structure formation ( [24,25,26,27]). Something that would confirm the ΛCDM model would be cosmological N-body simulations that produce galaxies that replicate the properties of observed galaxies ( [9] and references therein). These observed properties of galaxies include the galaxy luminosity function ( [28,29,30,31]), the relative frequencies of the various Hubble types ( [32]), the colour bimodality ( [33,34,35]), and -perhaps most importantly -dynamical properties elegantly encapsulated by MOND. Two of the most significant of these dynamical properties are (i) the baryonic Tully-Fisher relation ( [36,37]) and (ii) the scaling between the observed density profile of baryons and the inferred gravitational field (from a dynamical measure - [38]). This second property includes the requirement for dark matter (DM) to have a centrally cored distribution, rather than cusped ( [39,40]).
Although there are promising efforts to satisfy the aforementioned criteria ( [41,42,40,43,44,45]), these last two requirements remain illusive. The most common belief is that more sophisticated modelling of the complex hydrodynamical feedback processes, from central black holes and massive stars, will allow the observations and predictions to be reconciled ( [46,47,48]).
Another approach is to assume MOND describes the dynamics of galaxies well and that we must blend a combination of MOND and DM together. The first attempt of this type came from [12] who added MOND and 2 eV active neutrinos together. This was not able to produce the measured acoustic peaks in the CMB, nor satisfy the requirement of DM in MOND clusters of galaxies ( [49,13]) due to phase space constraints. [50] suggested combining MOND with an 11 eV sterile neutrino, which was shown to be consistent with the CMB (under the ansatz that MOND is irrelevant at z > 10 3 ) and the phase-space constraints from clusters of galaxies ( [14]). The added attraction of this idea (which was reviewed by [51]) was that the 11 eV neutrinos would be fully thermalised, i.e. half of all quantum states would be occupied (just like the active neutrinos) -removing any fine tuning of their abundance. The other reason MOND plus sterile neutrinos is so attractive is that the free streaming properties of neutrinos on small scales means they would not influence galaxies. Thus, they would leave the impressive results of MOND in galaxies unblemished.
Using a purpose built MOND cosmological N-body code, [52] showed this model was inconsistent with the observed cluster mass function -it produced too many high mass clusters and too few low mass clusters. [53] further demonstrated this behaviour for MOND cosmological simulations to produce too many superclusters was not limited to 11 eV sterile neutrinos. In fact, all masses of sterile neutrinos were ruled out, regardless of (i) the value of the MOND acceleration constant, (ii) the redshift at which MOND is "switched on", (iii) the normalisation of the initial conditions and (iv) the interpolating function. Other factors like the equation of state of dark energy were also shown to be ineffective.
In summary it does not appear to be possible to form the correct halo mass function in standard MOND from any sterile neutrino initial conditions that grew from an initially Harrison-Zel'dovich power spectrum under GR until z ∼ 200. So if MOND is the correct description of gravitational dynamics on galaxy scales, then either the initial conditions are not as described above and yet conspire to produce the correct CMB angular power spectrum, or perhaps MOND does not affect the sterile neutrinos. This is important because although galaxies require MOND to form (and stably exist) without CDM, the clusters clearly do not require MOND at all and one should not ignore how well Newtonian gravity reproduces the cluster mass function.
At minimum, the ΛCDM model gives the correct cluster scale halo mass function at z = 0, whether some additional boost to gravity is required to form the clusters early enough has been discussed in the literature ( [54,55,56,57,58,59,60,61]). MOND has a double negative effect on the cluster mass function if it influences the sterile neutrinos. Not only does it facilitate more rapid growth and the formation of much larger and denser structures than in Newtonian gravity, but these more massive halos now have MOND gravity meaning their dynamical masses are further enhanced, causing poorer agreement with the data. This result might suggest that if the MOND gravitational field is not produced by the sterile neutrinos (meaning only a Newtonian gravitational field is produced by them), but is only produced by the baryons, that it will have a positive influence on the halo mass function. Below we address how significant this influence is.
There is another factor to consider here which rules out the idea of MOND not being activated until some low redshift. That is that galaxies in MOND must form without the aid of a dark matter halo (cold, warm or hot) and galaxy formation without dark matter (if it is possible at all) is only possible with the added benefit of stronger than Newtonian gravitational attraction between the baryons. Thus, if MOND was not in effect until z = 1, then galaxies would not begin to form until then and galaxies are clearly formed long before this.
The simple options to blend MOND with DM are therefore ruled out. There are also more involved ideas, such as the tantalising theory of dipolar dark matter [62,63] and bimetric MOND ( [64,65]), which require further investigation. However, a simpler framework may still exist with sterile neutrinos and MOND.
In the traditional framework of MOND, there exists no dark matter in galaxies. Therefore, it is assumed that only the baryons produce a modified gravitational field. In the extended frameworks, where MOND is blended with some species of neutrinos, it has always been assumed that the baryons and the neutrinos contribute to the modified gravitational field. On the other hand, there do exist modified theories of gravity (completely different to MOND) where the dark matter has a mutual fifth-force interaction, to which the baryons do not participate (e.g. [66,67]). Here we go the other way, and suggest the baryons produce a modified gravitational field, but the sterile neutrinos do not. This can in principle retain all the benefits of adding sterile neutrinos to MOND -such as addressing the acoustic peaks of the CMB and solving the missing mass problem in MOND clusters. It also might evade the problems engendered when the sterile neutrinos produce a modified gravitational field -like overproducing massive clusters. Furthermore, in the original MOND + sterile neutrino framework, the effectiveness of MOND had to be "switched off" (i.e. the acceleration constant of MOND had to be substantially decreased) prior to recombination in an ad hoc way. This was to avoid MOND altering the shape of the CMB acoustic power spectrum. In a model where the sterile neutrinos do not produce a modified gravitational field, this may not be necessary.
In this article we run MOND cosmological simulations in a framework where the modified gravitational field is produced only by the baryons, not the sterile neutrinos. We present the framework and describe the simulation setup in §2, show the results of our simulations in §3 and in §4 we give our conclusions.

Background and Method
The Quasi-Linear formulation of MOND (QUMOND; [68]) requires solution of a modified version of the Poisson equation. Specifically, the ordinary Poisson equation for cosmological simulations is solved to give the Newtonian potential, Φ N , at scale factor a, from the ordinary matter density ρ that includes baryons and neutrinos. This would also include cold dark matter if there was any in our model. The QUMOND potential, Φ, is found from the Newtonian potential as follows and y = ∇Φ N /a o a, with a o being the MOND acceleration constant chosen here to be 3.6 ( km s −1 ) 2 pc −1 . The interpolating (ν) function is parametrised as per [2] Eqs. 51 and 53 where and α = 1 is the so-called simple ν-function and α = 2 is the standard ν-function.
The specifics of how to solve Eqs 2.1 & 2.2 are also explained in AD11, but we review the main points here.
The code we wrote is particle-mesh based, with a grid-mesh that has a limit of 257 cells in each dimension. In all our simulations we use 256 particles per dimension. Particlemesh solvers are required to handle the non-linear MOND Poisson equation. Direct or tree-code methods fail because in MOND we cannot co-add particle gravities. The full MOND Poisson equation must be solved because the trivial MOND equation of spherical symmetry does not satisfy the conservation laws.
Our one-component simulations are described by the following procedure • The particle positions and velocities are read in and the density of the particles is assigned to the various cells with the cubic cloud-in-cell method.
• We use multi-grid methods (see Numerical Recipes §19.6) to solve the Poisson equation to find the Newtonian potential (Eq 2.1). A 3D black-red sweep to update the cells with the new approximation of the potential in that cell and we iterate until we have fractional accuracy of 10 −10 .
• We take the divergence of the vector in the square brackets of Eq 2.2 which gives us the source of the QUMOND potential.
• We then repeat the Poisson solving step with the new source density to give the QUMOND potential, Φ, which we take the gradient of to find the gravity at each cell.
• We then interpolate to each particle's position to find the appropriate gravity and move each particle with a second order leapfrog.
• The procedure repeats from the second stage until the simulation reaches z = 0.

Modified Baryonic Dynamics
The modification of gravity we propose is the following: the Newtonian potentials of the sterile neutrinos and baryons are found by solving the Newtonian Poisson equation The gravitational field produced by the sterile neutrinos is not modified and is thus −∇Φ N,νs . On the other hand, the gravitational field produced by the baryons is found from the analogue of the QUMOND equation (Eq 2.2), which is where crucially The total modulus of the gravitational field experienced by any particle, sterile neutrino or baryon, is Given that a modified gravitational field (relative to Newtonian gravity) is only produced by the baryons, we refer to this model as Modified Baryonic Dynamics (MBD). Note that, at this point, it is a phenomenological approach at the classical level: whether a Lagrangian producing these equations of motion can be found (both at the classical and covariant level) could be the subject of further work, if necessary. It is plausible that such a Lagrangian will require couplings of the baryons and sterile neutrinos to two different metrics, thereby leading to a violation of the weak-equivalence principle, which would have additional interesting consequences in their own right (see, e.g., [69]).

Initial Conditions and Simulation Setup
In [52] and [53] we made use of the original COSMICS/GRAFICS package of [70] to generate our initial conditions. We chose to input our own transfer functions using the massive neutrino parametrisation of [71] (their equations 10-12) and the resulting linear matter power spectra were plotted for a sample of neutrino masses in [53] Fig 1. We exploit the new CMB results provided by the PLANCK mission ( [20]) which identifies Ω b , Ω νs , Ω Λ , h, n s )=(0.049, 0.267, 0.683, 0.671, 0.962).
We use the quadrupole temperature, Q rms−P S , as a free parameter to fit the amplitude of the cluster mass function in MBD.
The CMB quadrupole, Q rms−P S , is used to normalise the initial power spectrum of perturbations in the same way as σ 8 typically is for CDM simulations, because one cannot use linear theory in MOND to estimate σ 8 at z = 0.
In order to exploit the MBD formulation we use the GRAFICS-2 package of [72] which generates two sets of particles: baryons and sterile neutrinos (or any other DM particle one chooses), giving 256 3 for each species. The masses of the sterile neutrino and baryonic particles are found from m = 1.4 × 10 11 where Ω i should be replaced with either Ω νs or Ω b . It is then perfectly straight-forward to code the above expressions.
We then proceed as follows: • We read in the positions of the sterile neutrinos, compute their density (ρ νs ) on the grid and then solve for their Newtonian potential (Φ N,νs , as per Eq 2.4) using finite differencing and multi-grid methods as described in AD11.
• We read in the positions of the baryons, compute their density on the grid (ρ b ) and then their Newtonian potential (Φ N,b , as per Eq 2.4).
• We derive ∇ · [ν(y)∇Φ N,b ], the QUMOND source density for the baryons, which uses an argument for the interpolating function as per Eq 2.6.
• We find the QUMOND potential for the baryons Φ b , using Eq 2.5.
• We add the QUMOND potential for the baryons together with the Newtonian potential for the sterile neutrinos (as in Eq 2.7) to give the total potential, from which the gravitational field is derived to move both sets of particles.
• We then interpolate to each particle's position to find the appropriate gravity and move each particle with a second order leapfrog.
• The procedure repeats from the second stage until the simulation reaches z = 0.

Results
Using two component simulations with sterile neutrinos, the baryons almost exactly trace the sterile neutrino distribution, which is not surprising. As with the original single component MOND simulations ( [53]), when computing the mass function, we use the Newtonian equivalent mass of halos -that is the dynamical mass our MOND halos would appear to have if probed using Newtonian dynamics. With MBD it is merely a case of separately running the AMIGA halo finder ( [73,74]) on both the baryons and sterile neutrinos (which gives identical halo numbers and halo particle mass distributions for our cluster sized halos) and calculating the Newtonian equivalent mass, M m (r) = M b (r)ν GM ν+b (r) r 2 ao + M ν (r) since the modification to gravity only affects the baryonic gravitational field. M b and M ν represent the particle mass profiles of baryons and sterile neutrinos respectively. From here it is merely a case of finding the radius where M m (r) encloses an average density of matter 200 times the critical density and building the mass function. We only consider the mass bins of our mass function to be relevant if more than 8 halos are found within the bins of 0.233 dex per decade of mass.
We ran 12 two component simulations with different initial normalisations through the CMB quadrupole and different sterile neutrino masses. For each set of parameters, we ran two simulations: one with a 128 Mpc/h box and the other with a 256 Mpc/h box. We did this because we do not have adaptive mesh refinement built into our code. Therefore, in the simulations with the larger box, we have insufficient spatial resolution which leads to a suppression of the mass function. On the other hand, with the smaller boxes we naturally form fewer halos and have poor statistics. This issue is described in more detail in [53]. Since we have twice the spatial resolution in the 128 Mpc/h simulations, we use the amplitude of their mass functions to re-normalise the 256 Mpc/h simulations by matching their amplitude to the 128 Mpc/h boxes at a single mass bin. From comparison with theoretical mass functions ( [53]), we know the 128 Mpc/h box simulations produce roughly the correct mass function normalisation.
In Fig 1 we plot the mass functions from simulations with the same neutrino mass of 100 eV, but various normalisations through the CMB quadrupole. The red lines represent the 128 Mpc/h boxes and the black lines represent the renormalised mass functions of the 256 Mpc/h boxes. We only plot the curves for mass bins that have eight or more halos (as mentioned previously). The four normalisations used, in increasing amplitude, represent Q rms−P S = 4, 5, 9 and 17 µK. Clearly Q rms−P S < 5 µK is required to have a mass function amplitude which is comparable to the data (for a 100 eV sterile neutrino).
In general, the mass functions found using different normalisations, through the CMB quadrupole (Q rms−P S ), in pure MOND showed no variation (Fig 7, [53]) i.e. the mass functions had the same amplitude at z = 0 regardless of initial normalisation. The mass functions found with MBD clearly vary with normalisation.
In Fig 2 we plot the mass functions from simulations with the same normalisation through the CMB quadrupole, of Q rms−P S = 4.5 µK, but various neutrino masses 50, 100, 150 and 300 eV. We find that the lower the sterile neutrino mass, the lower the amplitude of the mass function and better the agreement with the observed cluster mass function. We only tested sterile neutrino masses of 50, 100, 150 and 300 eV, but using a mass of 30 eV or lower would not allow a fit to the cluster mass function. However, regardless of mass, a very low normalisation is always required.
To test if there was any sensitivity to the ν function (Eq 2.3) employed, we ran two further MBD simulations using α = 2 (the standard ν function), where the other MBD simulations used α = 1 (the simple ν function). We chose Q rms−P S =4.5 and 12.5 µK and the result was that there is very little difference between simulations with α = 1 and α = 2. This means using any reasonable ν function in MBD will still require a low normalisation for the initial spectrum of perturbations, through the CMB quadrupole.

The significance of the CMB quadrupole
The measured CMB quadrupole value is related to the multipole moment l = 2 coef- . [75] and [76] plot likelihood versus quadrupole value in their Figs 3 and 38 respectively. The allowed range at one sigma has shifted towards lower values in the latest paper and currently the lower one sigma confidence level is C T T 2 ≈ 50µK 2 (the maximum likelihood is around 160µK 2 ). This leads to Q rms−P S ≈ 4.5µK at one sigma and roughly 8 µK at maximum likelihood.
The fitted CMB quadrupole we require to be consistent with the galaxy cluster mass function is Q rms−P S < 5µK, which is roughly one sigma from the observed quadrupole.
The measured quadrupole is obviously separate from the fitted quadrupole. A quadrupole value is fitted to the CMB angular power spectrum.
In the ΛCDM model, a fitted quadrupole of around 21.4 µK ( [76]) is required so that the amplitude of the theoretical CMB power spectrum at mulitpoles l > 2 matches the observed spectrum. There should also be consistency between the theoretical value of the quadrupole required to fit the CMB and other cosmological probes such as baryonic acoustic oscillations, galaxy clustering and the galaxy cluster mass function i.e. σ 8 should be the same for each probe. The theoretical quadrupole is between one and two sigma larger than the measured one. However, consistency with the observed one can be achieved by models that suppress the theoretical one (due to inflationary theories, alternative gravity theories, etc) or increase the observed quadrupole (due to foreground modelling).
In the MBD case, we also expect the fitted quadrupole used to normalise the CMB (21.4 µK as per the ΛCDM value) and found from matching the galaxy cluster mass function to be the same. Clearly, the galaxy cluster mass function requires Q rms−P S < 5µK and the CMB requires Q rms−P S ≈ 21.4µK, so they are hugely discrepant. This is far more serious than a single, consistent fitted quadrupole overestimating the measured one -as in the ΛCDM case.
Therefore, our basic assumption that the CMB power spectrum would be the same in MBD as in ΛCDM must be revised. Instead, to check the validity of the MBD model, a full treatment of CMB anisotropies is required and if a model is found which matches the CMB power spectrum, its normalising CMB quadrupole amplitude should be consistent with the one used to fit the galaxy cluster mass function.

Conclusion
In previous works ( [52,53]) we tested the hypothesis that combining MOND with sterile neutrinos could produce the observed mass function of clusters of galaxies. We found it could not. This meant that, if MOND is the correct description of weak-field gravitational dynamics on galaxy scales, then either the whole cosmology and/or the initial conditions are not as described above (see, e.g., Sect. 9.2 of [2] for a discussion in the context of covariant MOND theories), and yet would conspire to produce the correct CMB angular power spectrum, or perhaps MOND does not affect the sterile neutrinos in the same way as the baryons. For this reason, we propose that the modified gravitational field of MOND is only produced by the baryons. In this scenario, the sterile neutrinos produce a Newtonian gravitational field. However, the baryons -subject to the total Newtonian gravitational field (both from the sterile neutrinos and baryons) -produce a MOND-like gravitational field.
There are similarities between our MBD model and many of the fifth force models referenced in [77]. For instance [67] employed a long range scalar field, generated only by the dark matter (not by the baryons), to increase the mutual acceleration of the two clusters comprising the colliding bullet cluster. In our MBD model, only the baryons generate a stronger than Newtonian gravitational attraction, whilst the sterile neutrinos produce a purely Newtonian gravitational field.
Our MBD model is far more conducive to producing a halo mass function that satisfies observational constraints than our previous attempts. Nevertheless, it requires fine tuning the normalisation of the primordial power spectrum, as is often done, through the CMB quadrupole, Q rms−P S . It was found that a normalisation of Q rms−P S < 5 µK was required to have an adequate match and a lower sterile neutrino mass was preferred over a higher one.
The CMB quadrupole required is around 1 σ lower than the measured value from CMB observations ( [76]). More significantly, this is much lower than the normalisation required to fit the CMB in ΛCDM (21.4 µK). This is crucial because we previously adopted the ΛCDM fits to the CMB as a basis for any MOND plus sterile neutrino theory, under the assumption the two models would give an identical CMB power spectrum.
We conclude that the simple model considered here is not viable, and that more involved MOND models should be considered if the MOND phenomenology is indeed derived from a fundamental modification of the Lagrangian of Nature rather than an emergent phenomenon at galaxy scales. These theories must include at least one new degree of freedom corresponding to a non-trivial dark matter fluid with a non-trivial coupling to baryons. Interesting possibilities in this vein include dipolar dark matter [62,63,78] or BIMOND with twin matter [64,65], as well as other possible frameworks still to be conceived.

Acknowledgements
GWA's research is supported by the FWO -Vlaanderen.