## Abstract

In this paper, a series of high-resolution *N*-body simulations is presented in which the equations of motion have been changed to account for modified Newtonian dynamics (MOND). It is shown that a low-Ω_{0} MONDian model with an appropriate choice for the normalization σ_{8} can lead to clustering properties at redshift *z* = 0 similar to the commonly accepted (standard) Λ cold dark matter (ΛCDM) model. However, such a model shows no significant structures at high redshift, with only very few objects present beyond *z* > 3 that can be readily ascribed to the low Ω_{0} value adopted. The agreement with ΛCDM at redshift *z* = 0 is driven by the more rapid structure evolution in MOND. Moreover, galaxy formation appears to be more strongly biased in MONDian cosmologies. Within the current implementation of MOND, density profiles of gravitationally bound objects at *z* = 0 can still be fitted by the universal Navarro, Frenk & White (NFW) profile, but MOND haloes are less clumpy.

## Introduction

Although the currently favoured ΛCDM model has proven to be remarkably successful on large scales [cf. *Wilkinson Microwave Anisotropy Probe* (*WMAP*) results, Spergel et al. 2003], recent high-resolution *N*-body simulations seem to be in contradiction with observation on subgalactic scales: the cold dark matter (CDM) ‘crisis’ on small scales is far from being over. The problem with the steep central densities of galactic haloes, for instance, is still unsolved as the highest-resolution simulations favour a cusp with a logarithmic inner slope for the density profile of approximately −1.2 (Power et al. 2003), whereas high-resolution observations of low surface brightness galaxies are best fitted by haloes with a core of constant density (de Block & Bosma 2002; Swaters et al. 2003). Suggested solutions to this include the introduction of self-interactions into collisionless *N*-body simulations (e.g. Bento et al. 2000; Spergel & Steinhardt 2000), replacing CDM with warm dark matter (e.g. Colin, Avila-Reese & Valenzuela 2000; Bode, Ostriker & Turok 2001; Knebe et al. 2002), or by non-standard modifications to an otherwise unperturbed CDM power spectrum (e.g. bumpy power spectra: Little, Knebe & Islam 2003; tilted CDM: Bullock 2001). Some of the problems – for instance, the overabundance of satellites – can be resolved with such modifications, but none of the proposed solutions have been able to rectify *all* shortcomings of ΛCDM simultaneously.

Therefore, there might be alternative solutions worthy of exploration, one of which is to abandon dark matter completely and to adopt the equations of modified newtonian dynamics (MOND; Milgrom 1983; Bekenstein & Milgrom 1984). It has already been shown by other authors that this simple idea might explain many properties of galaxies without the need of non-baryonic matter (e.g. Begeman, Broeils & Sanders 1991; Milgrom 1994; Sanders 1996; McGaugh & de Blok 1998; Scarpa, Marconi & Gilmozzi 2003). MOND is also successful in describing the dynamics of galaxy groups and clusters (Milgrom 1998; Sanders 1999), globular clusters (Scarpa et al. 2003) and, to a limited extent, gravitational lensing (Qin & Zou 1995; Mortlock & Turner 2001). A recent review of MOND is given by Sanders & McGaugh (2002), which also summarizes (most of) the successes alluded to above.

However, there has yet to come a detailed study of the implications of MOND in cosmological simulations of structure and galaxy formation, which is the aim of the current study. Nusser (2002) already investigated MOND of the large-scale structure using the *N*-body approach. His simulations, however, are of lower resolution, both in terms of spatial and mass resolution, making a study of individual objects difficult. Moreover, his implementation of the MOND equations is slightly different to our treatment.

The outline of the paper is as follows. In Section 2, we present the way we modified our *N*-body code mlapm to account for MOND. Section 3 introduces the cosmological models under investigation, whereas Section 4 summarizes the numerical details. The analysis can be found in Section 5. We finish with a summary and our conclusions in Section 6.

## The Mondian Equations of Motion

In an *N*-body code, one integrates the (comoving) equations of motion

In these equations, ** x** =

**/**

*r**a*is the comoving coordinate,

**is the canonical momentum, ρ(**

*p**x*) is the local comoving density, the mean comoving density, ∇

_{x}· is the divergence operator (Δ

_{x}is the nabla operator) with respect to

**and**

*x*

*F*_{pec}(

**) = −∇Φ(**

*x***) is the**

*x**peculiar acceleration field in comoving coordinates*. We now need to modify these (comoving) equations to account for MOND.

Despite MOND being a modification to Newton's second law rather than to gravity, one has the option to actually interpret MOND as an alteration of the law of gravity (cf. Sanders & McGaugh 2002). In that case, Poisson's equation

is replaced by where**=**

*r**a*

**is the proper coordinate,**

*x**g*

_{0}is the fundamental acceleration of MOND and

*g*_{M}is the MONDian acceleration field. Note that equations (3) and (4) are given in proper coordinates, in contrast to equation (2), where the solution

*F*_{pec}describes only the

*peculiar*acceleration.

If we now compare equations (3) and (4), we find the relation between MONDian acceleration *g*_{M} and Newtonian acceleration ** g** to be

**decreases like**

*h***(**

*O**r*

^{−3}) with increasing scale. Moreover, in cases of spherical, planar or cylindrical symmetry, the curl term vanishes. Under the assumption that objects forming in the Universe show at least one of these symmetries, we are able to neglect the curl term completely. One might argue that this kind of symmetry is probably very weak at high redshifts. However, later in this Section, we also make the assumption that MOND only affects peculiar accelerations (in proper coordinates) which are well above the MOND acceleration at early times (cf. equation 12). Therefore, we presume that ∇×

**is unimportant for the growth of large-scale structures as well as for the internal properties of virialized objects (under the assumption that they are at least symmetric along one of their axes).**

*h*Using Milgrom's (1983) suggested interpolation function

one now only needs to solve to get*g*

_{M}as a function of

*g*. The relevant solution of equation (7) is

Equation (8) allows us to obtain the MONDian acceleration *g*_{M} for a given Newtonian acceleration ** g**.

However, we are actually solving equation (2) in our *N*-body code mlapm, which gives *F*_{pec}, the Newtonian peculiar acceleration in *comoving coordinates*. Therefore, we also need to derive a relation between the proper acceleration , the proper peculiar acceleration *g*_{pec} and the peculiar acceleration in comoving coordinates *F*_{pec}. The second derivative with respect to time of ** r** =

*a*

**gives**

*x**peculiar acceleration in proper coordinates*

*g*_{pec}should be defined as where

*F*_{pec}is the solution of Poisson's equation (equation 2) as obtained by mlapm.

If we now make the assumption that MOND *only* affects the peculiar acceleration in proper coordinates but leaves the Hubble acceleration unchanged, the recipe for getting the MONDian peculiar accelerations in comoving coordinates *F*_{M,pec} is as follows:

- (i)
solve equation (2) using mlapm, which gives

*F*_{pec}; - (ii)
use equation (12) to transfer

*F*_{pec}to*g*_{pec}=*F*_{pec}/*a*^{2}; - (iii)
transfer

*g*_{M,pec}back to*F*_{M,pec}=*a*^{2}*g*_{M,pec}; - (iv)
use equation (8) to calculate

*g*_{M,pec}from*g*_{pec}; - (v)
use

*F*_{M,pec}rather than*F*_{pec}in equation (1).

This scheme has been employed for the simulations carried out with mlapm described below.^{1} We would like to point out that our treatment of MOND agrees with the one advocated by Sanders (2001) in the limit β= 0.

## The Cosmological Models

The reason for introducing MOND by Milgrom (1983) was to explain the flat rotation curves of galaxies *without* the need for dark matter. Having this in mind, we decided to use an Ω_{0} value for the MOND simulation that is close to the upper bound allowed by big bang nucleosynthesis and agrees with the recent measurements of cosmological parameters by the *WMAP* experiment (Spergel et al. 2003). Our data base of simulations is made up of the following three runs:

- (i)
a standard ΛCDM model;

- (ii)
an open, low-Ω

_{0}model with the same σ_{8}as ΛCDM; - (iii)
an open, low-Ω

_{0}model with MOND and adjusted σ_{8}.

These runs are labelled ΛCDM, OCBM and OCBMond, respectively, and their cosmological parameters are summarized in Table 1. The OCBM model is only to be understood as a gauge for the MOND model run, rather than as an alternative to ΛCDM.

We can see from equation (12) that the peculiar accelerations (which are subject to MOND) are large at early times, and therefore a ‘Newtonian treatment’ in the early Universe is justified. Thus, the input power spectra to our initial conditions generator were calculated using the cmbfast code (Seljak & Zaldarriaga 1996) with Ω_{0} = Ω_{b} = 0.04 for OCBM and OCBMond, respectively. This explains the choice for using the expression CBM rather than CDM. Such a relatively high Ω_{b} value (compared to Ω_{0}) actually introduces ‘baryon wiggles’ into the primordial power spectrum (oscillations frozen into the plasma at the epoch of recombination which are suppressed in dark-matter-dominated models). However, fluctuations on scales of the box size employed (*B* = 32 *h*^{−1} Mpc, see below) and smaller are not affected by it other than an overall ‘damping’ (e.g. Silk 1968; Eisenstein & Hu 1998). This damping, however, is compensated by the choice of normalization σ^{norm}_{8} of the power spectrum. Moreover, we differentiate between σ^{norm}_{8} and the actual measure of σ^{z = 0}_{8} in the simulations at redshift *z* = 0 because the OCBMond model requires a lower *P*(*k*)-normalization σ^{norm}_{8} to arrive at a comparable σ^{z = 0}_{8} value. This is due to a much faster growth of structures when using MOND, which will be emphasized in more detail in Section 5.

## The *N*-Body Simulations

Using the input power spectra derived with the cmbfast code, we displace 128^{3} particles from their initial positions on a regular lattice using the Zel'dovich approximation (Efstathiou, Frenk & White 1985). The box size was chosen to be 32 *h*^{−1} Mpc on a side. This choice guarantees proper treatment of the fundamental mode, which will still be in the linear regime at *z* = 0 (the scale that turns, non-linear at *z* = 0 is roughly 20 *h*^{−1} Mpc for the models under investigation). The particles were evolved from redshift *z* = 50 until *z* = 0 with the open-source adaptive mesh refinement code mlapm (Knebe, Green & Binney 2001). We employed 500 steps on the domain grid built of 256^{3} cells, and in all three runs, a force resolution of 11 *h*^{−1} kpc was reached in the highest density regions. The mass resolution of the runs is *m*_{p} = 1.30 × 10^{9} and 0.17 × 10^{9}*h*^{−1} M_{⊙} for the ΛCDM model and for the two low-Ω_{0} models, respectively. We output the particle positions and velocities at redshifts *z* = 5, 3, 1, 0.5 and 0. These ‘snapshots’ are then analysed with respect to the large-scale clustering as well as the properties of individual objects. Even though in OCBM(ond) we made the assumption of the non-existence of dark matter, we still refer to these objects as ‘haloes’; due to the absence of dissipation, they show similar shapes and density distributions when being compared to their ΛCDM counterparts, as can be seen in the analysis below.

Gravitationally bound objects are identified using the bound density maxima method (BDM, Klypin & Holtzman 1997). The BDM code identifies local overdensity peaks by smoothing the density field on a particular scale of the order of the force resolution. These peaks are prospective halo centres. For each of these halo centres, we step out in (logarithmically spaced) radial bins until the density reaches ρ_{halo}(*r*_{vir}) = Δ_{vir}ρ_{b}, where ρ_{b} is the background density. This defines the outer radius *r*_{vir} of the halo. However, one needs to choose the correct virial overdensity Δ_{vir} carefully, which is much higher for the OCBM and OCBMond model due to the low Ω_{0} value. The parameters used are Δ_{vir} = 340 for ΛCDM and Δ_{vir} = 2200 for OCBM/OCBMond (see Gross 1997, Appendix C, and references therein).

Brada & Milgrom (1999) raised the issue that care should be taken when numerically integrating the equations of motion using MOND. Even more care should be taken in our case, where we made the assumption that the curl term in equation (5) vanishes, which gives rise to forces that are not strictly conservative. Therefore, we performed a simple test on the final output to ensure that our implementation of MOND does not violate momentum conservation: the net momentum of all particles in the box as well as the net force should vanish. We calculated the sums

and derived values for*C*

_{1}and

*C*

_{2}in all three models. In ΛCDM and OCBM, those constants

*C*

_{1,2}are of the order 10

^{−4}; even though in OCBMond they are about an order of magnitude larger, we still believe that they are sufficiently small. It has been shown elsewhere (Knebe et al. 2001) that adaptive softening in general does not guarantee precise momentum conservation and our values for

*C*

_{1}and

*C*

_{2}are as expected.

## Analysis

### Large-scale clustering properties

We start with inspecting the large-scale density field in all three runs. Fig. 1 shows a projection of the whole simulation, with each individual particle grey-scaled according to the local density at redshifts *z* = 0 and 5. This figure indicates that the MOND simulation looks fairly similar to the other two models in terms of the locations of high-density peaks (dark areas), filaments and the large-scale structure, respectively. However, one should bear in mind that the OCBMond simulation was started with a much lower σ^{norm}_{8} normalization than the other two runs. In fact, this is reflected in the right panel showing the density field at redshift *z* = 5; the MOND simulation is less evolved. Moreover, without any further analysis, one might even be inclined to conclude from Fig. 1 that the MOND model is more strongly clustered at *z* = 0, which is affirmed by the slightly higher σ^{z = 0}_{8} value given in Table 1. As we will see later on, this is not necessarily true; we can confirm a higher amplitude of the two-point correlation function for galaxies in OCBMond (cf. Fig. 5 later) while at the same time showing a comparable amplitude in the matter power spectrum on small scales. The latter can be viewed in Fig. 2, where we plot the matter power spectrum for all three models at redshift *z* = 0. There we observe that the OCBMond model shows a slightly larger amplitude for *k* < 0.8 *h* Mpc^{−1} (scales close to the fundamental mode). Nusser (2002), who treated the MOND equation similarly to (but differently from) us,^{2} already pointed out that the linear evolution of the growing mode solution for the density contrast δ scales like δ∝*a*^{2}, as opposed to Newtonian theory, where δ∝*a*. This explains why the OCBMond model with the (initially) low σ^{norm}_{8} normalization outruns the evolution of the Newtonian OCBM simulation. In other words, the MOND model had to be started with a lower σ_{8} normalization to provide competitive results at redshift *z* = 0. Sanders (2001) also noted that, due to a much faster growth of structures in MOND universes, the amplitude of *P*(*k*) for purely baryonic models matches the standard ΛCDM cosmology. We also like to point out that the required value σ^{norm}_{8} = 0.4 needed to bring the MOND simulation into agreement with the standard ΛCDM model is closer to the *COBE* normalization σ^{COBE}_{8}≈ 0.1 for that particular cosmology. Another feature worth noting is that the OCBMond power spectrum does not show a distinctive ‘break’ due to the transfer of power from large to small scales, as seen in both the ΛCDM and the OCBM model.

In Fig. 3, we are plotting the cumulative mass function of objects identified using the BDM technique (Klypin & Holtzman 1997). This figure highlights again that the hierarchical formation of gravitationally bound objects is driven much faster in the MOND simulation, but is being initiated at later times. At a redshift of *z* = 5, we can see that the abundance of objects on all mass-scales is nearly identical in the OCBM and ΛCDM model, with a much lower amplitude for OCBMond. Whereas the evolution for the Newtonian OCBM run already ceases at a redshift of around *z*_{stop}≃ 1/Ω_{0}− 1 ≃ 24, we still see a very strong increase (by more than one order of magnitude) in the number density of objects of all masses in the OCBMond simulation. To emphasize this, Fig. 4 shows the (integral) abundance evolution of objects with mass *M* > 10^{11}*h*^{−1} M_{⊙}. The OCBM model undoubtedly experiences very little evolution from *z*∼ 5 to today, whereas both other models show a very steep evolution. The discrepancy between the OCBMond and the other two models at redshift *z* = 5 can, however, be ascribed to the lower initial σ^{norm}_{8} value; OCBMond was set up with much smaller initial density perturbation which only grew to a comparable level of clustering via the effects of MOND. Again, this is in agreement with the findings of Sanders (2001), who showed that the collapse of spherically symmetric overdensities becomes MOND dominated for redshifts *z*≲ 5 and hence starts to outrun Newtonian models (cf. fig. 5 in Sanders 2001).

The question now arises as to what degree the (formation) sites of haloes in ΛCDM and the two OCBM models are correlated. To this extent, we calculated the two-point correlation function of the 500 most massive objects; the result can be found in Fig. 5. We chose to use a fixed number of haloes rather than a mass cut for the calculation of ξ_{gal}(*r*) in order not to introduce an artificial bias. We do have far more objects of a given mass in the ΛCDM model, and therefore a mass limit *M*_{min} used with the estimator for ξ_{gal}(*r*) would lead to different correlation amplitudes. The agreement between the Newtonian OCBM and the ΛCDM model is not surprising. As already pointed out by other authors, the correlation function is expected to be (nearly) identical in cases of equal σ_{8} normalizations, irrespective of the cosmological model (Martel & Matzner 2000). Moreover, if the model is fixed and only the σ_{8} normalization is varied, it should leave no imprint on ξ_{gal}(*r*) (Croft & Efstathiou 1994). However, the OCBMond model stands out again, having a much higher amplitude on all scales. This clearly attributes for the differences already mentioned in the discussion of Fig. 1 and indicates that the OCBMond model is more evolved at redshift *z* = 0 than the other two models, even though the matter power spectrum has a lower amplitude on corresponding scales; structure formation in MONDian cosmologies is even more biased than in Newtonian models. This is in agreement with findings that small scales enter the MOND regime before large-scale fluctuations (cf. Nusser 2002; Sanders 2001).

### Galactic haloes

Having analysed the large-scale clustering properties, we now turn to the investigation of the internal properties of galactic haloes. To this extent, we use the halo catalogues based on the BDM code (Klypin & Holtzman 1997) and neglect objects less massive than *M* < 10^{11}*h*^{−1} M_{⊙} again. This sets the minimum number of particles per halo to 77 for ΛCDM and 488 for OCBM(ond).

A visualization of the density fields throughout the most massive BDM halo is given in Fig. 6. It is quite striking that neither of the low-Ω_{0} haloes shows substantial substructure. However, this is easily understood for OCBM, because in a low-density universe, structure formation ceases at early times (*z*_{stop}≃ 1/Ω_{0}− 1). This means that clusters in such cosmologies should show fewer substructures because they had more time to virialize (cf. Knebe & Müller 2000). However, this explanation obviously does not hold for OCBMond, as nearly all haloes in that particular model formed exceptionally late (cf. Fig. 4).

The most interesting question, however, is probably the shape of the density profile and the rotation curve for haloes in MONDian cosmologies. Fig. 7 now shows the matter profile of the most massive halo in all three models along with fits (thin solid lines) to a Navarro, Frenk & White (NFW) profile (Navarro, Frenk & White 1997),

The scale radius*r*

_{s}is being used to define the concentration of the halo

*c*=

*r*

_{vir}/

*r*

_{s}, where

*r*

_{vir}is the radius at which the density reaches the virial overdensities Δ

_{vir}≈ 340 and ≈2200 for ΛCDM and OCBM(ond), respectively.

We observe that – even for the OCBMond model – the data is equally well described by the functional form of a NFW profile (out to the virial radius, at least). However, the central density of that halo in the OCBMond model is lower than in ΛCDM and especially than in OCBM. The high central density for OCBM is readily explained by the fact that haloes in that particular cosmology form at a time when the Universe is still very dense. This result is also supported by the values of the concentration parameter presented in Table 2: the most massive object in the MOND model shows the lowest concentration *c*, mostly due to the late onset of formation, as observed in Fig. 4. When inverting the density profile into a (Newtonian) circular velocity curve simply by using *v*^{2}_{circ}(*r*) =*GM*(<*r*)/*r*, this also entails a shift of the maximum of the rotation curve to higher radii, as can be seen in Fig. 8 (*r*^{circ}_{max}≈ 2*r*_{s}, Navarro et al. 1997). The functional shape of the rotation curves for an NFW density profile is given by

*v*

_{vir}being the (Newtonian) circular velocity at the virial radius

*r*

_{vir}and where

*x*=

*r*/

*r*

_{vir}. The drop of the maximum of the rotation curve by about a factor of 3.2 is merely a reflection of the scatter in mass for the most massive halo throughout the three models. As can be seen in Table 2, the halo is more than ten times as massive in ΛCDM than in OCBMond. This should give an approximately 2.7 times higher

*v*

^{circ}

_{max}, as the scaling between those two quantities is roughly

*v*

^{circ}

_{max}∝

*M*

^{1/3}

_{vir}. This scaling relation can be derived when using

*x*

^{circ}

_{max}≈ 2/

*c*(cf. Bullock et al. 2001b; Navarro et al. 1997) together with equation (15), giving Furthermore, if we assume

*c*∝

*M*

^{−0.13}, as shown by Bullock et al. (2001b, cf. equation 18), we find that the factor under the square root in equation (16) is roughly constant for the mass range under consideration. Also, as

*v*

_{vir}∝

*M*

^{1/3}(which simply follows from

*v*

^{2}∝

*M*/

*r*and

*r*∝

*M*

^{1/3}for a spherical overdensity), the same scaling then holds for

*v*

^{max}

_{circ}, explaining the observed drop of the maximum of the rotation curve for the OCBMond model.

However, for the OCBMond model, we should take into account that accelerations are *not* Newtonian and hence *v*^{2}_{circ}(*r*) =*GM* (<*r*)/*r* does not hold. Therefore, we also show in Fig. 8 the actual MONDian rotation curve calculated as follows. The Newtonian acceleration *g* =*v*^{2}/*r* is transfered to the MONDian acceleration *g*_{M} according to equation (8). In turn, *g*_{M} is used to calculate . The resulting *v*_{circ}(*r*) is labelled ‘MOND’ in Fig. 8. We note that the MONDian velocities are actually larger than the Newtonian ones, bringing them closer to the ΛCDM model. This has implications for dynamical mass estimates of galaxy clusters as presented in Sanders (1999). Sanders showed that the dynamically estimated mass of galaxy clusters is too large compared with the observed mass when using Newtonian physics. However, MOND rectifies dynamical masses, bringing them into better agreement with observed masses. A similar phenomenon can be observed in Fig. 8: the MONDian curve for OCBMond would be measured observationally and hence be translated into too high a cluster mass of *M*_{vir}≈ 2.8 × 10^{13}*h*^{−1} M_{⊙} when using Newtonian dynamics. MOND, on the contrary, would give the real value of *M*_{vir} = 0.3 × 10^{13}*h*^{−1} M_{⊙}. Later on, we will see that this leaves an imprint on the (radial) distribution of the spin parameter, too.

In addition to the rotation curve, we also show the acceleration as a function of radius in Fig. 9. The object that formed in OCBMond has MONDian accelerations throughout the whole halo, whereas the central parts of the objects in the other two models are above the universal acceleration *g*_{0} indicated by the thin solid line.

We would like to stress that in both Figs 8 and 9, the Newtonian curves for the OCBMond model are only plotted for completeness; they do *not* carry observable information, as the halo actually follows MONDian rather than Newtonian physics. However, they are educational in the sense of gauging the importance MOND has on the internal structure of the halo.

Table 2 now lists some internal properties in addition to the ones already mentioned, namely the velocity dispersion σ_{v}, the virial radius *r*_{vir}, the triaxiality *T*, ellipticities *e*_{1} and *e*_{2}, the spin parameter λ as well as the best-fit parameters λ_{0} and σ_{λ} when fitting the probability distribution *P*(λ) to a lognormal distribution. The spin parameter λ was calculated using the definition given in Bullock et al. (2001a),

*v*

_{vir}is again the MONDian circular velocity at the virial radius. The distribution of λ,

*P*(λ), has been fitted to a lognormal distribution (e.g. Frenk et al. 1988; Warren et al. 1992; Cole & Lacey 1996; Gardner 2001; Maller, Dekel & Somerville 2002) The results are presented in Fig. 10, showing that the OCBMond distribution

*P*(λ) is nearly indistinguishable from the other two models. However, the reduced χ

^{2}value for OCBMond is about a factor of two larger (cf. Table. 2). As has been noted by Maller et al., the lognormal distribution given by equation (18) is not as good a fit to models with haloes that are recent merger remnants. This is definitely one of the effects that has an influence on the spin parameter distribution for the OCBMond model, as we expect a high level of recent merger activity (cf. Fig. 4).

We also calculated the radial distribution of λ(<*r*) throughout the haloes and present the result for the most massive one in Fig. 11. We see that λ(<*r*) is roughly constant for the Newtonian models of structure formation, whereas there is a sharp increase of λ(<*r*) in the MOND halo towards the virial radius. This implies that the material in that particular halo moves on more circular orbits and that the object itself is closer to solid body rotation, respectively. These results are in agreement with our previous finding, where we showed that the MONDian rotation curve is not given by simply ‘inverting the density profile’, as in the Newtonian case (cf. Fig. 9); the velocity in the outskirts of the halo is much larger and hence we expect a rise in λ(*r*), too.

We are now going to focus on the shape of the haloes. First, we show measurements of the overall shape of haloes; to this extent, we calculated the eigenvalues *a* > *b* > *c* of the inertia tensor. They were used in turn to construct the triaxiality parameter (e.g. Franx, Illingworth & de Zeeuw 1991)

*T*are accompanied by the ellipticities and the values for the most massive object are summarized in Table 2. The mean values 〈

*T*〉, 〈

*e*

_{1}〉, and 〈

*e*

_{2}〉, when averaging over all haloes more massive than 10

^{11}

*h*

^{−1}M

_{⊙}, can be found in Table 3. We observe a trend for the MOND haloes to be more triaxial with higher ellipticities at the same time.

Secondly, we would like to quantify subclumps within the virial radius of the halo itself. One possibility to measure the substructure content of a halo is to calculate the radial profile of the density dispersion

where ρ_{i}(

*r*) is the local density at a particle position which was estimated using the nearest 20 neighbours and 〈ρ(

*r*)〉 is the average taken over all

*N*(

*r*) particles within a spherical shell [

*r*,

*r*+

*dr*]. We used the same binning as already applied to the density profile and the rotation curve, respectively. The result for the most massive halo in all three runs is presented in Fig. 12. This figure shows that the dispersion is always smallest for the MOND halo, at least out to the virial radius; the density at each individual particle position is always close to the mean density. If there were subclumps present, one would expect to find peaks in the σ

^{2}

_{δ}(

*r*) curve due to local deviations from the mean density profile.

## Summary and Discussion

In this paper, we presented three cosmological simulations: a fiducial standard ΛCDM model, a very low-Ω_{0} model, and the same Ω_{0} model but with MONDian equations of motions. Putting aside the classical arguments against MOND, we showed that structures found in a MONDian low-Ω_{0} universe are not significantly different from the standard ΛCDM model. However, we derived several differences which might easily validate or rule out MONDian cosmologies observationally. Our main results can be summarized as follows.

- (i)
Even though the OCBMond run was set up using a low value for σ

^{norm}_{8}= 0.4, the simulation finished with σ^{z = 0}_{8}= 0.92 (which is close to the cluster normalization of the ΛCDM model). - (ii)
The OCBMond model shows an extremely fast evolution of the number density of haloes for redshifts

*z*< 5. - (iii)
Virialized objects in a MONDian universe are slightly more correlated.

- (iv)
The density profile of MOND objects still follows the universal NFW density profile, but they are

- (a)
far less concentrated;

- (b)
show less substructure; and

- (c)
are closer to solid body rotation.

We conclude that the most distinctive feature of a MONDian universe is the late epoch of galaxy formation; in a MONDian universe that is normalized to match a ΛCDM model at redshift *z* = 0, we expect *not* to observe any galaxies until recently (*z* < 5). This actually holds for any low-σ_{8} universe, either MONDian or Newtonian, but only the assumption of MOND can bring such a model into agreement with observations of the local Universe again. Another interesting finding is that the outer parts of MOND haloes are closer to solid body rotation than their standard ΛCDM counterparts, even though the overall distribution *P*(λ) is nearly indistinguishable from the ΛCDM model.

However, there have still been many assumptions made during the course of this study which are hard to justify, and hence all results are to be understood simply as preliminary until our understanding of MOND actually improves. Not only did we assume that MOND only affects peculiar accelerations in proper coordinates, but we also neglected the curl term in equation (5). The effect of this rotational component is not clear, but it is guaranteed to decrease rapidly on large scales and vanish completely for objects that have spherical, planar or cylindrical symmetry. Another disclaimer is that MOND is based on the non-existence of dark matter, but our simulations only model gravity; if the Universe consists only of baryons, one definitely would need to include gas physics to be able to make credible predictions for internal properties of galaxies. However, this is beyond the scope of this study. None the less, we have shown that cosmology with MOND does not necessarily lead to completely odd results. Our treatment of MOND, though, gives clustering properties comparable to the favourite concordance ΛCDM model.

## Acknowledgments

The simulations presented in this paper were carried out on the Beowulf cluster at the Centre for Astrophysics & Supercomputing, Swinburne University. We acknowledge the support of the Australian Research Council through its Discovery Project program (DP0343508). We would like to thank the two anonymous referees for valuable comments that helped to improve the scientific content of the paper.