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.
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 these equations, x =r/a is the comoving coordinate, p is the canonical momentum, ρ(x) is the local comoving density, 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) and (4) are given in proper coordinates, in contrast to equation (2), where the solution Fpec describes only the peculiar acceleration. (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(7) is
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 , the proper peculiar acceleration gpec and the peculiar acceleration in comoving coordinates Fpec. The second derivative with respect to time of r =ax gives(1) leads to (9) as 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:
solve equation (2) using mlapm, which gives Fpec;
use equation (12) to transfer Fpec to gpec =Fpec/a2;
transfer gM,pec back to FM,pec =a2gM,pec;
use equation (8) to calculate gM,pec from gpec;
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:
a standard ΛCDM model;
an open, low-Ω0 model with the same σ8 as ΛCDM;
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 σ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 sumsKnebe et al. 2001) that adaptive softening in general does not guarantee precise momentum conservation and our values for C1 and C2 are as expected.
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.
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).
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).
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).
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),
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 byTable 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 vcircmax∝M1/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 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 vvir∝M1/3 (which simply follows from v2∝M/r and r∝M1/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.
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 . 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.
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),Frenk et al. 1988; Warren et al. 1992; Cole & Lacey 1996; Gardner 2001; Maller, Dekel & Somerville 2002) 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)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.
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 dispersionFig. 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.
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).
The OCBMond model shows an extremely fast evolution of the number density of haloes for redshifts z < 5.
Virialized objects in a MONDian universe are slightly more correlated.
The density profile of MOND objects still follows the universal NFW density profile, but they are
far less concentrated;
show less substructure; and
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.
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.