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  

(1)
formula
which are completed by Poisson's equation  
(2)
formula

In these equations, x =r/a is the comoving coordinate, p is the canonical momentum, ρ(x) is the local comoving density, graphic the mean comoving density, ∇x· is the divergence operator (Δx is the nabla operator) with respect to x and Fpec(x) = −∇Φ(x) is the 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  

3
formula
is replaced by  
4
formula
where r =ax is the proper coordinate, g0 is the fundamental acceleration of MOND and gM is the MONDian acceleration field. Note that equations (3) and (4) are given in proper coordinates, in contrast to equation (2), where the solution Fpec describes only the peculiar acceleration.

If we now compare equations (3) and (4), we find the relation between MONDian acceleration gM and Newtonian acceleration g to be  

(5)
formula
The field equation (4) is non-linear and is difficult to solve in general. However, it has been shown by Bekenstein & Milgrom (1984) that ∇×h decreases like 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 ∇×h 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).

Using Milgrom's (1983) suggested interpolation function  

(6)
formula
one now only needs to solve  
(7)
formula
to get gM as a function of g. The relevant solution of equation (7) is  
(8)
formula

Equation (8) allows us to obtain the MONDian acceleration gM for a given Newtonian acceleration g.

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

(9)
formula
whereas combing equation (1) leads to  
(10)
formula
Using the second Friedmann equation, we can then rewrite equation (9) as  
(11)
formula
This equation shows that the peculiar acceleration in proper coordinatesgpec should be defined as  
(12)
formula
where Fpec 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 FM,pec is as follows:

  • (i)

    solve equation (2) using mlapm, which gives Fpec;

  • (ii)

    use equation (12) to transfer Fpec to gpec =Fpec/a2;

  • (iii)

    transfer gM,pec back to FM,pec =a2gM,pec;

  • (iv)

    use equation (8) to calculate gM,pec from gpec;

  • (v)

    use FM,pec rather than Fpec 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.

Table 1.

The best-fitting model in each magnitude bin using all the data.

Table 1.

The best-fitting model in each magnitude bin using all the data.

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 σnorm8 of the power spectrum. Moreover, we differentiate between σnorm8 and the actual measure of σz = 08 in the simulations at redshift z = 0 because the OCBMond model requires a lower P(k)-normalization σnorm8 to arrive at a comparable σz = 08 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 1283 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 2563 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 mp = 1.30 × 109 and 0.17 × 109h−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(rvir) = Δvirρb, where ρb is the background density. This defines the outer radius rvir 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  

(13)
formula
and derived values for C1 and C2 in all three models. In ΛCDM and OCBM, those constants C1,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 C1 and C2 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 σnorm8 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 = 08 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 δ∝a2, as opposed to Newtonian theory, where δ∝a. This explains why the OCBMond model with the (initially) low σnorm8 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 σnorm8 = 0.4 needed to bring the MOND simulation into agreement with the standard ΛCDM model is closer to the COBE normalization σCOBE8≈ 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.

Figure 1.

Comparison of the large-scale density field of the three models under investigation at redshifts z = 0 (left-hand panels) and 5 (right-hand panels). The upper panels show the ΛCDM simulation, the middle panels show the OCBMond model and the lower panels show the OCBM model.

Figure 1.

Comparison of the large-scale density field of the three models under investigation at redshifts z = 0 (left-hand panels) and 5 (right-hand panels). The upper panels show the ΛCDM simulation, the middle panels show the OCBMond model and the lower panels show the OCBM model.

Figure 5.

Two-point correlation function at redshift z = 0 for the 500 most massive objects.

Figure 5.

Two-point correlation function at redshift z = 0 for the 500 most massive objects.

Figure 2.

Matter power spectra for all three models at redshifts z = 0 (thick lines) and 5 (thin lines).

Figure 2.

Matter power spectra for all three models at redshifts z = 0 (thick lines) and 5 (thin lines).

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 zstop≃ 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 > 1011h−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 σnorm8 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).

Figure 3.

Cumulative mass function for all three models at redshifts z = 5, 3, 1 and 0.

Figure 3.

Cumulative mass function for all three models at redshifts z = 5, 3, 1 and 0.

Figure 4.

Redshift evolution of the abundance of haloes with mass M > 1011h−1 M.

Figure 4.

Redshift evolution of the abundance of haloes with mass M > 1011h−1 M.

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 Mmin 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 < 1011h−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 (zstop≃ 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).

Figure 6.

Most massive galactic halo in ΛCDM (upper panel), OCBMond (middle panel) and OCBM (lower panel). The line in the lower right corner of each panel is a scale indicating 100h−1 kpc.

Figure 6.

Most massive galactic halo in ΛCDM (upper panel), OCBMond (middle panel) and OCBM (lower panel). The line in the lower right corner of each panel is a scale indicating 100h−1 kpc.

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),  

(14)
formula
The scale radius rs is being used to define the concentration of the halo c =rvir/rs, where rvir is the radius at which the density reaches the virial overdensities Δvir≈ 340 and ≈2200 for ΛCDM and OCBM(ond), respectively.

Figure 7.

Density profile for the most massive halo in all three models. The vertical thin lines indicate the virial radii.

Figure 7.

Density profile for the most massive halo in all three models. The vertical thin lines indicate the virial radii.

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 v2circ(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 (rcircmax≈ 2rs, Navarro et al. 1997). The functional shape of the rotation curves for an NFW density profile is given by  

(15)
formula
with vvir being the (Newtonian) circular velocity at the virial radius rvir and where x =r/rvir. 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 vcircmax, as the scaling between those two quantities is roughly vcircmaxM1/3vir. This scaling relation can be derived when using xcircmax≈ 2/c (cf. Bullock et al. 2001b; Navarro et al. 1997) together with equation (15), giving  
(16)
formula
Furthermore, if we assume cM−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 vvirM1/3 (which simply follows from v2M/r and rM1/3 for a spherical overdensity), the same scaling then holds for vmaxcirc, explaining the observed drop of the maximum of the rotation curve for the OCBMond model.

Table 2.

Properties of the most massive gravitationally bound object. The mass M is measured in h−1 M?, velocities are measured in km s−1 and radii are measured in h−1 kpc.2 is actually the reduced ?2 value as returned by the curvefit routine of IDL.a

Table 2.

Properties of the most massive gravitationally bound object. The mass M is measured in h−1 M?, velocities are measured in km s−1 and radii are measured in h−1 kpc.2 is actually the reduced ?2 value as returned by the curvefit routine of IDL.a

Figure 8.

Rotation curve of the most massive halo.

Figure 8.

Rotation curve of the most massive halo.

However, for the OCBMond model, we should take into account that accelerations are not Newtonian and hence v2circ(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 =v2/r is transfered to the MONDian acceleration gM according to equation (8). In turn, gM is used to calculate graphic. The resulting vcirc(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 Mvir≈ 2.8 × 1013h−1 M when using Newtonian dynamics. MOND, on the contrary, would give the real value of Mvir = 0.3 × 1013h−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 g0 indicated by the thin solid line.

Figure 9.

Acceleration curve of the most massive halo. The solid vertical line indicates the MONDian accleration parameter g0 = 1.2 × 108 cm s−2.

Figure 9.

Acceleration curve of the most massive halo. The solid vertical line indicates the MONDian accleration parameter g0 = 1.2 × 108 cm s−2.

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 rvir, the triaxiality T, ellipticities e1 and e2, 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), 

(17)
formula
which proved to be a more stable measurement than the usual graphic definition. For the OCBMond model, vvir 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)  
(18)
formula
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).

Figure 10.

Spin parameter distribution in all three models. Lines show fits obtained using the lognormal distribution given by equation (18).

Figure 10.

Spin parameter distribution in all three models. Lines show fits obtained using the lognormal distribution given by equation (18).

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.

Figure 11.

Radial distribution of spin parameter λ. Again, the vertical lines show the virial radius.

Figure 11.

Radial distribution of spin parameter λ. Again, the vertical lines show the virial radius.

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)  

(19)
formula
The triaxialities T are accompanied by the ellipticities  
(20)
formula
and the values for the most massive object are summarized in Table 2. The mean values 〈T〉, 〈e1〉, and 〈e2〉, when averaging over all haloes more massive than 1011h−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.

Table 3.

Mean triaxiality parameter T and mean ellipticities e1, e2 when averaging over haloes with M > 1011h−1 Mc.

Table 3.

Mean triaxiality parameter T and mean ellipticities e1, e2 when averaging over haloes with M > 1011h−1 Mc.

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  

(21)
formula
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.

Figure 12.

Variance σδ(r) for the most massive halo. The vertical lines indicate the virial radius of the respective halo.

Figure 12.

Variance σδ(r) for the most massive halo. The vertical lines indicate the virial radius of the respective halo.

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 σnorm8 = 0.4, the simulation finished with σz = 08 = 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.

References

Bahcall
N.
et al
,
2003
,
ApJ
 ,
585
,
182
Begeman
K. G.
Broeils
A. H.
Sanders
R. H.
,
1991
,
MNRAS
 ,
249
,
523
Bekenstein
J.
Milgrom
M.
,
1984
,
ApJ
 ,
286
,
7
Bento
M. C.
Bertolami
O.
Rosenfeld
R.
Teodoro
L.
,
2000
,
Phys. Rev. D
 ,
62
Bode
P.
Ostriker
J. P.
Turok
N.
,
2001
,
ApJ
 ,
556
,
93
Brada
R.
Milgrom
M.
,
1999
,
ApJ
 ,
519
,
590
Bullock
J. A.
,
2001
, ()
Bullock
J.
Dekel
A.
Kolatt
T. S.
Kravtsov
A. V.
Klypin
A. A.
Porciani
C.
Primack
J. R.
,
2001
a
ApJ
 ,
240
,
555
Bullock
J.
Kolatt
T. S.
Sigad
Y.
Sommervile
R. S.
Kravtsov
A. V.
Klypin
A. A.
Primack
J. R.
Dekel
A.
,
2001
b
MNRAS
 ,
559
,
321
Cole
S.
Lacey
C.
,
1996
,
MNRAS
 ,
281
,
716
Colin
P.
Avila-Reese
V.
Valenzuela
O.
,
2000
,
ApJ
 ,
542
,
622
Croft
R. A. C.
Efstathiou
G.
,
1994
,
MNRAS
 ,
267
,
390
Davis
M.
Efstathiou
G.
Frenk
C. S.
White
S. D. M.
,
1985
,
ApJ
 ,
292
,
371
de Blok
W. J. G.
Bosma
A.
,
2002
,
A&A
 ,
385
,
816
Dressler
A.
Shectman
S. A.
,
1988
,
Astron. J.
 ,
95
,
985
Eisenstein
D. J.
Hu
W.
,
1998
,
ApJ
 ,
496
,
605
Efstathiou
G.
Frenk
C. S.
White
S. D. M.
,
1985
,
ApJS
 ,
57
,
241
Franx
M.
Illingworth
G.
de Zeeuw
T.
,
1991
,
ApJ
 ,
383
,
112
Frenk
C. S.
White
S. D. M.
Davis
M.
Efstathiou
G.
,
1988
,
ApJ
 ,
327
,
507
Gardner
J.
,
2001
,
ApJ
 ,
557
,
616
Gross
M. A. K.
,
1997
,
PhD thesis
 ,
Univ. California
,
Santa Cruz
Klypin
A. A.
Holtzman
J.
,
1997
, ()
Knebe
A.
Müller
V.
,
2000
,
A&A
 ,
354
,
761
Knebe
A.
Green
A.
Binney
J.
,
2001
,
MNRAS
 ,
325
,
845
Knebe
A.
Devriendt
J.
Mahmood
A.
Silk
J.
,
2002
,
MNRAS
 ,
329
,
813
Little
B.
Knebe
A.
Islam
R. R.
,
2003
,
MNRAS
 ,
341
,
617
McGaugh
S. S.
de Blok
W. J. G.
,
1998
,
ApJ
 ,
499
,
66
Maller
A. H.
Dekel
A.
Somerville
R. S.
,
2002
,
MNRAS
 ,
329
,
423
Martel
H.
Matzner
R.
,
2000
,
ApJ
 ,
530
,
525
Milgrom
M.
,
1983
,
ApJ
 ,
270
,
365
Milgrom
M.
,
1994
,
ApJ
 ,
270
,
371
Milgrom
M.
,
1998
,
ApJ
 ,
496
,
L89
Mortlock
D. J.
Turner
E.
,
2001
,
MNRAS
 ,
327
,
557
Navarro
J.
Frenk
C. S.
White
S. D. M.
,
1997
,
ApJ
 ,
490
,
493
Nusser
A.
,
2002
,
MNRAS
 ,
331
,
909
Pinkney
J.
Roettiger
K.
Burns
J. O.
Bird
C. M.
,
1996
,
ApJ
 ,
104
,
1
Power
C.
Navarro
J. F.
Jenkins
A.
Frenk
C. S.
White
S. D. M.
Springel
V.
Stadel
J.
Quinn
T.
,
2003
,
MNRAS
 ,
338
,
14
Qin
B.
Zou
Z. L.
,
1995
,
A&A
 ,
296
,
264
Sanders
R. H.
,
1996
,
ApJ
 ,
473
,
117
Sanders
R. H.
,
1999
,
ApJ
 ,
512
,
L23
Sanders
R. H.
,
2001
,
ApJ
 ,
560
,
1
Sanders
R. H.
McGaugh
S. S.
,
2002
,
A&A
 ,
40
,
263
Scarpa
R.
Marconi
G.
Gilmozzi
R.
,
2003
,
A&A
 ,
405
,
15
Scarpa
R.
,
2003
, ()
Seljak
U.
Zaldarriaga
M.
,
1996
,
ApJ
 ,
469
,
437
Silk
J.
,
1968
,
ApJ
 ,
151
,
459
Spergel
D. N.
Steinhardt
P. J.
,
2000
,
Phys. Rev. Lett.
 ,
84
,
3760
Spergel
D. N.
et al
,
2003
,
ApJS
 ,
148
,
175
Swaters
R. A.
Madore
B. F.
van den Bosch
F. C.
Balcells
M.
,
2003
,
ApJ
 ,
583
,
732
Warren
M. S.
Quinn
P. J.
Salmon
J. K.
Zurek
W. H.
,
1992
,
ApJ
 ,
399
,
405
1
The latest release version of mlapm includes the MOND implementation, which can be activated using-dmond on compile time.
2
Nusser (2002) did not use Milgrom's interpolation function, as given by equation (6), but rather a spontaneous transition from Newtonian to MOND accelerations.