EDP Sciences
Free Access
Issue
A&A
Volume 552, April 2013
Article Number A36
Number of page(s) 16
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201321037
Published online 21 March 2013

© ESO, 2013

1. Introduction

To control the numerical integration accuracy of a finite self-gravitating system evolution in a static background, one may use total energy and angular momentum conservation (Aarseth 2003).In the cosmological case, instead of being conserved, the total (peculiar) energy decreases with time. Indeed, the cosmic energy equation, known as the Layzer-Irvine (LI) equation, describes the change of the total peculiar energy of a system due to the adiabatic expansion of the background (Peebles 1980). The quantification of the deviations from the LI equation is one of the very first tests made in all cosmological codes (see, e.g., Teyssier 2001; Miniati & Colella 2007). However, when typical initial conditions of standard models of structure formation are considered, the application of the LI test is limited by several reasons. First of all, high redshift initial conditions of typical cosmological N-body simulations are very uniform and cold, so that both the initial peculiar potential and kinetic energies are close to zero: this makes difficult the numerical estimation of the cosmic energy equation. In addition, at low redshifts, the density field is typically characterised by a great variety of non-linear structures of different sizes and by low density regions still in the linear regime. Structures of different sizes and in different physical regimes give very different contributions to total the peculiar potential and kinetic energies so that it is not straightforward to estimate which error in the cosmic energy equation maybe tolerated and/or which constraints on the accuracy of the simulation maybe placed (see also discussion in Joyce & Sylos Labini 2013, and references therein).

For this reason several authors have developed convergence tests to explore the role of the various physical and numerical parameters that characterise a given simulation with the aim of establishing the reliability of the results (see e.g., Power et al. 2003). Alternatively there are studies that compare different algorithms (see e.g., Knebe et al. 2003; Iannuzzi & Dolag 2011, and references therein). However, it is not a simple issue to control the numerical accuracy of the non-linear gravitational clustering in a cosmological simulation; for instance the question of determining the optimal spatial resolution, i.e. the best value for the gravitational force softening length, is still an open problem (see e.g., Splinter et al. 1988; Power et al. 2003; Romeo et al. 2008; Joyce et al. 2009a).

In this paper, following the work of Joyce & Sylos Labini (2013), we consider the evolution of a very simple system, namely an isolated over-density in an expanding background, as a toy model to control the accuracy of the numerical integration and to test the correctness of a cosmological simulations.

The evolution of this kind of structure in a static background, has been studied since the earliest numerical N-body experiments (Hénon 1964; van Albada 1982; Aarseth et al. 1988; David & Theuns 1989; Theuns & David 1990; Bertin 2000; Heggie & Hut 2003; Boily & Athanassoula 2006; Joyce et al. 2009b; Sylos Labini 2012). Briefly, driven by its own mean field, it collapses in a time scale , forming a virialized quasi-stationary state (QSS), i.e. a collision-less configuration that is not in thermodynamical equilibrium but that has a finite lifetime fixed by the collisional relaxation time-scale (1)where N is the system particles number and κ is a numerical factor that takes into account the “Coulomb logarithm” (Chandrasekhar, S. 1943; Theis & Spurzem 1999; Diemand et al. 2004; Gabrielli et al. 2010).

Joyce & Sylos Labini (2013) have shown that an isolated over-density in a cosmological simulation should in principle evolve, in physical coordinates, just as if it were in a non-expanding space with open boundary conditions. This corresponds to the fact that such a structure is in the stable clustering regime. Namely, it evolves as it were truly isolated from the rest of universe. While the stable clustering regime is assumed to approximately describe the evolution of a structure in the strongly non-linear regime (Peebles 1980), in the case of a single over-density with periodic boundary conditions and negligible finite size effects it must be precisely followed: the stable clustering regime thus corresponds to a QSS for a structure in a static background.

Indeed, in a cosmological simulation, where periodic boundary conditions are used, a single over-density can be considered isolated if finite size effects are negligible. Such mass distribution is isolated but for the interaction with its “copies” included in the infinite system over which the force is summed. It is easy to show that the copies contribute to the force with terms that represent perturbations suppressed by positive powers of R0/L, where R0 is the over-density size and L is the side of the periodic box (Gabrielli et al. 2006; Joyce & Sylos Labini 2013). When the latter condition is satisfied, the only difference in the evolution between the expanding and the static case can arise from the smoothing of the gravitational force.

Beyond the comparison of the expanding case with the template given by the static background evolution, one may consider the analysis of the LI equation. In fact, Joyce & Sylos Labini (2013) found that is a useful diagnostic to control the accuracy of cosmological simulations when the evolution of a single over-density is considered. In this case, the breakdown of the LI test occurs when the structure leaves the stable clustering regime: this is due to the difficulty of numerically integrating the system evolution with sufficient accuracy because hard two-body collisions take place.

We consider in this paper the evolution ofan initially spherical and uniform cloud of particles with a velocity dispersion such that the initial virial ratio is in the range b0 ∈  [−1,0]. When the system is warm enough, i.e. b0 ≈  −1, it approximately maintains its original size and shape during and after the collapse. Instead, when the system is cold enough, i.e. b0 ≈ 0, it undergoes to a stronger contraction so that its size is considerably reduced during the collapse forming thereafter a virialized structure with a density profile decaying as r-4 (Joyce et al. 2009b; Sylos Labini 2012). Our goal is to study the collapse of such a single cloud of particles in an expanding universe: by monitoring the LI equation and the stability of the clustering (in real space) we may understand to which accuracy the cosmic energy equation must be satisfied, in the strongly non-linear regime, so that the simulation still provides with physically meaningful results. In particular we will study the accuracy of the numerical integration by changing the numerical and physical parameters that characterise a simulation.

The paper is organised as follows. In Sect. 2 we recall the main results of the evolution of an isolated over-density in a static background. Then we summarise the results of Joyce & Sylos Labini (2013) concerning the approximate equivalence of an isolated cloud of particles in a static background with open boundary conditions and the same cloud when it is placed alone in a box of an infinite periodic system with space expansion – i.e., when the simulation is run in cosmological conditions. In Sect. 3 we discuss the preparation of the initial conditions, the choice of the parameters that control the numerical integration and we briefly recall the LI test. The comparison between the evolution of an isolated over-density in a static and expanding background is presented in Sect. 4. A series of systematic tests to determine the role of the integration parameters such as the gravitational force softening length and the accuracy of the time-stepping is presented in Sect. 5. Then in Sect. 6 we discuss several tests performed by changing the physical parameters of the initial conditions, such as the number of particles, the system size and the velocity dispersion. Finally in Sect. 7 we discuss the principal results and then we draw our main conclusions.

2. Isolated over-density in a cosmological and a static background

Firstly we briefly recall the main features of the collapse of an isolated over-density in a static background with open boundary conditions (hereafter STAT case). Then we discuss the equivalence of the evolution of an isolated structure in STAT case and in an expanding background with periodic boundary conditions (EXP case).

2.1. Isolated over-densities in a static background

In the STAT case, the global collapse of an isolated and spherical cloud of self-gravitating particles, of density ρ, is determined by its mean gravitational field and it is characterised by the time scale (2)When the initial velocity dispersion is such that the virial ratio is b0 = 2K0/W0 ≈  − 1 (where K0/W0 is the system initial kinetic/potential energy) the collapse consists in a few oscillations that drive the system to a virialized QSS: this maintains approximately the size R0 and shape of the initial configuration (Sylos Labini 2012).

On the other hand, when the initial velocity dispersion is small enough (or zero), the system undergoes to a stronger contraction reducing its size by a large factor. In this case the features of the virialized QSS depend explicitly on the number of system particles (at constant density): for instance, when particles are initially Poisson distributed, the intrinsic characteristic size in Eq. (3) scales as rc ≈ N− 1/3 (Aarseth et al. 1988; Joyce et al. 2009b). Other noticeable features are: (i) during the collapse a fraction of particles gain enough kinetic energy to escape from the system; (ii) the density profile after the collapse is characterised by the power law decay described by (3)(that we denominate the quasi equilibrium profile – QEP) where rc and nc are two constants.

The key to understand the formation of this density profile is represented by the physics of the ejection mechanism. The main characteristics of the virialized QSS that can be framed in a simple physical model introduced by Sylos Labini (2012)are where vr is the radial velocity and Φ(r) is the so-called pseudo phase-space density. Note that, once the QSS is formed, its potential and kinetic energy reach a constant value: this fact represents one of the diagnostics to numerically detect the quasi stationarity of the particle configuration (Joyce et al. 2009b; Sylos Labini 2012; Joyce & Sylos Labini 2013).

2.2. Equivalence of the evolution in a static and in an expanding background

Joyce & Sylos Labini (2013) have shown that when a cosmological code is used to simulate an isolated structure (of size R0) in an expanding universe, it should reproduce the same results, in physical coordinates, as that obtained for the same structure in open boundary conditions without expansion. Let us briefly recall in which limit this equivalence is verified.

Dissipation-less cosmological N-body simulations solve numerically the equations (Peebles 1980) (7)where xi are the comoving particle coordinates of the i = 1,...,N particles of equal mass m, enclosed in a cubic box of side L, and subject to periodic boundary conditions, a(t) is the appropriate scale factor for the cosmology considered, and H(t) = ȧ/a is the Hubble constant. Note that the force on a particle is that due to the N − 1 other particles and to all their copies. The leading non-zero term of the force due to the copies is a repulsive “dipole”; by neglecting the quadrupole and higher multipole terms, that are convergent and suppressed by positive powers of (R0/L) compared to the dipole term (Gabrielli et al. 2006; Joyce & Sylos Labini 2013), we can rewrite Eq. (7) as (8)where the sum is now only over the N − 1 particles in the simulation box. Assuming for simplicity an Einstein-de Sitter (EdS) cosmology, for which (9)Eq. (7) may be written, in physical coordinates ri ≡ a(t)xi, simply as (10)These are the equations of motion of N purely self-gravitating particles. Therefore, up to finite size corrections that vanish in the limit R0/L → 0, the equations of motion (in physical coordinates) in an expanding universe are the same of those describing the evolution of the same finite structure in a static background with open boundary conditions.

We stress that Eq. (10) represents an important approximation that justifies the analogy between the same cloud of self-gravitating particles in STAT and EXP conditions. This is not a trivial approximation as it consists in replacing the sum in Eq. (7) on a infinite periodic system, i.e. with infinite terms, with the sum on a finite system in Eq. (10). The periodic sum can be written as a local sum on N particles and a non-local term which takes into account the contribution of the infinite copies: this latter term can be neglected for the system considered, i.e. for R0 < L. Only in this situation the evolution of a single over-density in a periodic system, i.e., an infinite number of over-densities one in each of the infinite boxes of the periodic system, in an expanding background is equivalent to the evolution of a truly isolated over-density in a static background with open boundary conditions.

Note that the simulations are run, as usual, in comoving coordinates. In what follows we present analyses both in comoving and in physical coordinates because the stable clustering regime should manifest itself in physical coordinates, i.e. the structure remains almost invariant only in these coordinates (Joyce & Sylos Labini 2013).

3. The simulations

3.1. Initial conditions

We have considered a system with total mass M in a spherical volume : this can be discretized by using N randomly distributed particles of mass m such that (11)We have used 103 ≤ N < 104 (series S), N = 104 (series M) and 104 < N ≤ 105 (series L). Details of the simulations are reported in Table 1: all length scales are given in units of the box side L. Note that the results we report require only choice of units for length and energy: for the former we will take units defined by L = 1, and for the latter units in which W0 = 3   GM2/(5   R0) =  −1, i.e. in which the initial continuum limit gravitational potential energy is minus unity.

Table 1

Details of the N-body simulations.

Note that the system radius R0 is chosen to be R0 < L/4 so to minimise, as discussed above, the effect of the periodic boundary conditions. In particular, the distance between any pair of structure particle i,j is less than that separating i from any particle of the copies in the infinite periodic system.

The spherical cloud corresponds, in a cosmological simulation, to an over-density with initial mass density contrast (12)where is the cloud density and ρ0 = M/L3 is the mean mass density of the universe (see Table 2). Particles mass is such that ρ0 is equal to the critical density.

Table 2

Parameters of the spherical clouds with different initial radius R0 and with over-density δ (Eq. (12)).

The scale factor in an EdS universe obeys to (Peebles 1980) (13)where . Assuming a0 = 1 from Eqs. (12), (13) we find (Joyce & Sylos Labini 2013) (14)where τc is given by Eq. (2). Note that from Eq. (14) we find that only for very high values of the over-density, i.e. δ > 104, at moderate expansion factors a ≈ 10, (t − t0)/τc is of the order of a thousand i.e. of the order of τ2 for N = 103 (see Eq. (2)).

3.2. Numerical integration parameters

We have used the publicly available tree-code GADGET (Springel et al. 2001; Springel 2005). In this code one may choose the values of a number of parameters and we have followed the prescriptions of the user guide for GADGET-21 to fix the basic ones.

The gravitational potential used by GADGET has the exact Newtonian potential above a separation of 2.8ε while for smaller separations the potential first derivative vanishes2. The different values of ε used are reported in Table 1. Note that we have chosen a softening length that is fixed in comoving coordinates and thus it varies in physical ones.

Further we have set, following Power et al. (2003); Springel (2005) a Barnes-Hut opening criterion with an opening angle of the tree algorithm θ = 0.7 for the first force computation and a dynamical updating criterion thereafter3. The accuracy of the relative cell-opening criterion (i.e., the parameter ErrTolForceAcc)4 is set to φ = 5 × 10-3. Finally, we have chosen the time-stepping criterion 0 of GADGET where the value of the parameter controlling its accuracy η (i.e., the parameter ErrTolIntAccuracy) is reported in Table 1.

3.3. Basic definitions

Let us introduce some basic definitions that will be useful in what follows. The kinetic energy in physical coordinates is (15)and the potential energy in physical coordinates is (16)where (17)and where g(r) is the exact GADGET pair potential (Springel 2005). The virial ratio is defined as (18)Note that the initial virial ratio is b0 = 2K0/W0, where K0,W0 are respectively the initial kinetic and potential energy. We define the potential energy in comoving coordinates (19)and the GADGET gravitational potential energy (20)i.e., the two-body potential is calculated by using periodic boundary conditions. As mentioned above, in the limit R0 ≪ L, we have Wg ≈ Wc.

The peculiar kinetic energy is (21)where (22)If the system is virialized in a region of size R0 ≪ L then (23)In this approximation (24)so that the kinetic energy in physical coordinates (Eq. (15)) is well approximated by the peculiar kinetic energy (Eq. (21)), i.e. Kp ≈ Kc.

3.4. The Layzer-Irvine equation

Let us briefly recall the expression of the LI equation. By defining the total peculiar energy as the LI equation can be written as (25)Integrating Eq. (25) we find (26)To test whether the LI equation is numerically satisfied during the time integration, one may measure the deviation of A(a) from 1. In particular, we will determine the expansion factor aLI such that (27)and larger than 0.1 for a > aLI. This allows us to quantitatively determine the range of redshifts in which the LI equation is well approximated in the simulations, i.e. in which the LI test is satisfied. Note that at very early times, when the initial velocities are close (or equal) to zero the denominator of Eq. (26) approaches zero and thus A(a) is difficult to be precisely determined. Otherwise for long enough times there is not any particular numerical problem in the determination of Eq. (26).

4. Comparison between the evolution in a static and expanding background

We first compare the evolution of an isolated and uniform over-density in an expanding background with periodic boundary conditions (EXP) with the evolution of the same system in a static background and with open boundary conditions (STAT). We consider here the simulation M4a (see Table 1), which we have chosen as a reference for the reasons illustrated below. Results for the other simulations are discussed in the next sections.

Note that the difficulty of integrating the cold collapse occurs only around the strongest collapse phase of the system, i.e. at t ≈ τc, a time scale which is much smaller than any relevant value of the expansion factor considered in the what follows.

4.1. Short time evolution

Figure 1 (upper panel) shows the potential and kinetic energy in physical coordinates (normalised to the initial potential energy W0). As for the STAT case the potential (kinetic) energy decreases (increases) to reach its minimum (maximum) at t ≈ τc, i.e. at the strongest phase of the collapse. For t > τc the bounded structure is stationary, i.e. Wp ≈ const., and it satisfies the virial equilibrium (Fig. 1 – bottom panel). A fraction of about fp ≈ 25% of the particles gain kinetic energy so that, after the collapse, they are ejected from the system moving as free particles thereafter.

Energy conservation for the STAT case can easily be controlled to a precision better than ≈ 0.01% (see, e.g., Aarseth (2003) and references therein). From the LI test (Fig. 1 – bottom panel) in the EXP case we find A(t) ≈ 1 with relatively large fluctuations up to the maximum contraction phase. However, these fluctuations are damped after the collapse: we thus conclude that in this range of time the system is still in the regime a < aLI (see Eq. (27)).

thumbnail Fig. 1

Upper panel: energy in physical coordinates (Kp Eq. (15) and Wp Eq. (16)) normalised to the initial total energy as a function of time (computed by using Eq. (14) for the EXP case), together with the corresponding behaviour of the same system in the STAT case. Bottom panel: virial ratio | b(t) |, together with the LI test (Eq. (26)) for the EXP case.

Open with DEXTER

The density profile (Fig. 2 – upper panel) is well fitted by Eq. (3) both in the EXP and STAT case. The radial velocity dispersion (Fig. 2 – bottom panel) shows the Keplerian behaviour described by Eq. (4) Therefore the pseudo phase-space density obeys to Eq. (5), i.e. it decays as Φ(r) ~ r− 5/2.

thumbnail Fig. 2

Comparison of the density (upper panel) and velocity (bottom panel) profiles of the same system evolved in an expanding background (EXP at a = 1.36 – black circles) and in a static background with open boundary conditions (STAT at t = 5τc – red squares). The best fit with Eq. (3) (QEP) and a line with slope − 1 are shown as references (see Eq. (4)).

Open with DEXTER

4.2. Long time evolution

The situation discussed in Sect. 4.1 drastically changes at longer times in the EXP simulation. Indeed, already for moderate values of the expansion factor, the evolution of the system rapidly deviates from the short-time behaviour, and correspondingly the properties of the QSS become different from Eqs. (3)–(5).

We find that Wp ≈ const., i.e. the structure is in a QSS and in the stable clustering regime, only for a < aLI ≈ 6 (see Fig. 3). Instead, for a > aLI we find that | Wp |  ~ 1/a (and thus Wc ≈ const). Thus the breakdown of the LI test corresponds to a departure from the stable clustering regime.

thumbnail Fig. 3

Upper panel: absolute value of the potential energy in physical coordinates Wp(a) (Eq. (16)) as a function of the scale factor. The approximate value of aLI, separating the regime where Wp ≈ const. from the one where | Wp |  ~ a-1, is reported as reference. Bottom panel: LI test (Eq. (26)) and virial ratio (absolute value) for particles with negative energy.

Open with DEXTER

thumbnail Fig. 4

Evolution of the density profile in comoving coordinates. Two reference lines with slope –4 and –3 are reported as a reference. The arrow shows the direction of the density profile evolution. Inset panel: behaviours in function of the scale factor of the length scale R(0.4) (R(0.8)), including the 40% (80%)of the mass of the virialized structure, normalised to 2.8ε.

Open with DEXTER

Indeed, the density profile in comoving coordinates for a > aLI stabilises to an almost “asymptotic” shape which does not change anymore at small radii and slowly changes at large ones (see Fig. 4). This is characterised by an approximate x-3 decay extending over about three decades, i.e. from ~ 10-3   L to L.

Note that the softening length is fixed in comoving coordinates and thus in physical coordinates it increases as εp = . When the structure is in the stable clustering regime, its size in physical coordinates Rp is constant, and thus Rp/εp ∝ 1/a: if the stable clustering regime persists long enough then the size of the structure will become of the order of, or even smaller than, the softening length. However, we noticed that beyond a certain time the structure becomes stable in comoving coordinate and thus its size in physical coordinates scales as , where is the constant size in comoving coordinates. Therefore for a > aLI one approximately finds that const.

The length scale can be estimated, for instance, by measuring the radius R(0.4) (R(0.8)) including the 40% (80%) of the mass of the virialized structure. As shown in inset panel of Fig. 4, for a > aLI, R(0.4) and R(0.8) remain respectively larger than 10 and 100 times the scale beyond which the force is exactly Newtonian, i.e. 2.8ε. In this condition the virial ratio fluctuates around –1 at all times (see Fig. 3). This implies that the Newtonian value of the mean field potential is well approximated by the smoothed potential and thus the softening length does not perturb the mean field dynamics. (Note that this result is confirmed by the behaviours measured by using other values of the softening length, as well as by an analysis of the gravitational potential energy in function of ε – see discussion in Sect. 5.1.)

Coherently with the behaviour of the energy, the density profile in physical coordinates (see Fig. 5) does not change for a < aLI, and it is well fitted by Eq. (3): this situation corresponds to the QSS in the stable clustering regime. Note that, as time passes, the spatial extension of the halo, i.e. the range of radii where n(r) ~ r-4, increases. Indeed, as time passes particles with energy ep ≲ 0 may go asymptotically far from the structure core (Sylos Labini 2012).

On the other hand, for a > aLI, we find that n(r) continues to evolve: at large radii it decays as ~r-3 and at small radii the density does not flatten as for the STAT case. A Navarro et al. (1997) (NFW) profile (28)gives a reasonable fit for a > aLI (bottom panels of Fig. 5): in particular the large scale tail shows a well-defined r-3 decay. Therefore only when stability in physical coordinates is observed, the density profile is consistent with the expected one, i.e. Eq. (3); when stability in physical coordinates is broken then the density profile changes its functional behaviour approaching the NFW one.

thumbnail Fig. 5

Density profile in physical coordinates for the simulation M4a at different scale factors a (see labels) together with the density profile for the STAT case. The arrow shows the direction of the density profile evolution in the cosmological case. Further it is shown the best fit given by Eq. (3) (QEP) and a reference line with slope –3.

Open with DEXTER

Similarly, the radial velocity dispersion profile in physical coordinates (Fig. 6) shows an approximated Keplerian decay at large scale as in the STAT case. However, while for a < aLI the different curves collapse on each other, implying that the virialized structure is in the stable clustering regime, for a > aLI the velocity profile marks as wella breaking of the stable clustering regime. In particular, the slope remains − 1, i.e. halo particles are still orbiting in a central gravitational potential but with an amplitude that decreases with time. This is consistent with the fact that the potential energy in physical coordinates Wp increases with the scale factor for a > aLI (as Wc ≈ const.): as time passes and thus Wp increases proportionally to − 1/a.

thumbnail Fig. 6

Radial velocity dispersion profile in physical coordinates at different scale factors a (see labels) together with the profile for the STAT case. The arrow shows the direction of the density profile evolution in the cosmological case. A line with slope r-1 (see Eq. (4)) is reported as reference.

Open with DEXTER

Thus for a > aLI there is an energy injection into the system which causes an increase of the total energy: this occurs slowly enough that the system can rearrange its phase space configuration to satisfy the virial theorem (see the bottom panel of Fig. 3). Therefore the energy increase corresponds to decrease of the kinetic energy: as b ≈  − 1 we have ΔE =  −ΔKp. Note the physical size of the structure increases almost simultaneously with the violation of the LI equation: this implies that the energy injection is spurious as otherwise the LI test should be verified.

The decrease of the kinetic energy corresponds to the fact that high velocity particles are slowed down. This is shown in Fig. 7 (upper panel) where it is reported the evolution of the probability distribution function (PDF) f(v) of the absolute value of the velocity. While for a < aLI we find that f(v) remains almost stable, for a > aLI high velocity particles are slowed down and correspondingly the high velocity tail of f(v) is significantly reduced. However, the shape of f(v) remains almost the same for a > aLI. The failure of the LI test and consequent breakdown of the stable clustering regime is thus induced by the difficulty of integrating high velocity particles.

thumbnail Fig. 7

Evolution of the probability distribution function (PDF) of the absolute value of the velocity at different scale factors (see labels). The arrow shows the direction of the evolution.

Open with DEXTER

High velocities are typically associated with close scatterings and the limited ability to resolve hard scattering between close neighbours, through which particles can increase their velocity, can also be traced by looking at the evolution of the nearest neighbour distribution ω(x) in comoving coordinates (see Fig. 7): indeed, one may note that while for a < aLIω(x) shows a peak that is shifted at smaller and smaller scale, corresponding to the fact that the structure is stable in physical coordinates, for a > aLIω(x) does not change anymore reaching an “asymptotic” shape. This is a clear effect of the limited small scale resolution of the numerical code.

Finally we note that from the behaviours of the density profile and of the radial velocity dispersion, we can easily find that the pseudo phase-space density, for a < aLI and for r > rc, scales as Φ(r) ~ r− 5/2, i.e., the behaviour expected for the QEP (see Eq. (5)). On the other hand, for a > aLI, at large radii, we find Φ(r) ~ r-3 × r3/2 ~ r-1.5.

4.3. Summary

Let us summarise the situation in the two regimes that we have identified, i.e. a < aLI and a > aLI.

  • 1.

    For a < aLI the structure is in a QSS (i.e., | Wp |  ≈ const.) and the density and velocity profiles are described by the collision-less equilibrium behaviours (Eqs. (3)–(5)). This situation corresponds to the stable clustering regime: the physical size and shape of the structure remains invariant during time evolution, modulo a marginal spatial extension of the power-law tail due to high velocity, but bounded, particles.

  • 2.

    For a > aLI the structure is (almost) frozen in comoving coordinates and therefore its size in physical coordinates increases and its potential energy tends to zero as Wp ∝  − 1/a. The breakdown of the IL test marks a problem in the numerical integration of the particles equations of motion in this regime: this occurs because there is an injection of energy which causes the expansion of the structure in physical space. Correspondingly the density profile develops a r-3 tail,approaching a NFW profile. The velocity dispersion profile still decays in a Keplerian way but with an amplitude that decreases with time equilibrating the smaller potential energy (in absolute value) of the structure, so that the virial ratio remains b ≈  − 1.

The decrease of the system kinetic energy for a > aLI corresponds to the fact that high velocity particles are slowed down and the PDF of the particles velocity f(v) shifts towards smaller velocity values. This is indeed the consequence of the inaccurate integration of two-body collisions: it becomes more difficult to numerically resolve high velocity particles. As shown by the behaviour of the nearest neighbour distribution ω(x) these high velocity particles are also close neighbours. Indeed, ω(x) does not evolve anymore for a > aLI, as it should if stable clustering were satisfied. We thus conclude that the inaccurate integration of hard collisions induces the change of regime at aLI.

thumbnail Fig. 8

Evolution of the nearest neighbour distribution in comoving coordinates at different scale factors (see labels). The arrow shows the direction of the evolution. The initial distribution (IC) is also shown for comparison.

Open with DEXTER

Note that aLI ≈ 6 corresponds to a time (see Eq. (14)) (tbreak − t0)/τc ≈ 50 that is much smaller than the collisional relaxation time scale (see Eq. (1)) τ2 ≈ 104τc for a system with N = 104 particles. This implies that the break-down of the QSS is not related to the expected evaporation due to two-body collisions for t > τc.

Finally we note that the breakdown of the QSS, if physically originated, should not be associated with the breakdown of the LI test as we find. This shows that such a breakdown is induced by spurious reasons that we identify in the role of the inaccurate integration of high velocities particles and thus of hard scatterings. We will return on this important point in Sect. 5.1.

5. Role of numerical integration parameters

In this section we study the effect of varying some of the most relevant numerical parameters that characterise the numerical integration: the softening length ε, the parameter that controls the time-stepping accuracy η, the opening angle θ of the tree algorithm and the accuracy of the relative cell-opening criterion φ5.

5.1. Softening length

The softening length6 has been varied in a range such that: (i) ε < R0, i.e. it is smaller than the initial size of the structure so that the gravitational mean field force is well approximated at least at early times and (ii) the lower limit to ε is due to numerical constraints: ε > 10-6 otherwise the integration time rapidly becomes too long and the simulation becomes unfeasible.

The potential energy in physical coordinates Wp and the IL test (Fig. 9) for simulations of the same system with N = 104 particles and with different values of ε (see Table 1) – changed by a factor 104 – show behaviours coherent with those discussed in Sect. 4.2. Namely the two regimes found above, corresponding to Wp ≈ const. for a < aLI and | Wp |  ~ a-1 for a > aLI, still characterise the temporal evolution. However, the transition between them depends on the softening length, i.e., aLI = aLI(ε).

In particular, we find that aLI ≈ 6 for ε ∈  [3.7 × 10-5,3.7 × 10-4], while outside this range aLI ≈ 1. This result is in agreement with Joyce & Sylos Labini (2013) who concluded that at most in that narrow range of ε can the behaviour required, a QSS in the stable clustering regime, be reproduced by the N-body method.

thumbnail Fig. 9

Behaviour of A(a) (Eq. (26)) for simulations with different ε (see Table 1 for details).

Open with DEXTER

To understand the effect of a softening length which is too large with respect to the size of the structure, we have computed the gravitational potential energy of the structure with density profile described by Eq. (3), in the case in which the pair potential is the GADGET one with different values ε. We find that the gravitational potential energy is very well approximated by the Newtonian value when R(0.4)/(2ε) ≈ 10 and/or R(0.8)/(2ε) ≈ 100, where R(0.4) and R(0.8) are the radii in comoving coordinates including respectively the 40% and 80% of the mass of the virialized structure.

The behaviours of R(0.8)/(2ε) as a function of the expansion factor and for different values of the softening length are reported in Fig. 10 (see the bottom panel). For ε ≥ 3.7 × 10-4 we have that R(0.4)/(2.8ε) < 10: in this situation the effective mean field force due to all particles is different from Newtonian, i.e. the smoothening length is too large with respect to the size of the system. The effect of a too large softening thus corresponds to a deviation of the virial ratio from –1 (see the upper panel of Fig. 10).

thumbnail Fig. 10

Upper panel: virial ratio as a function of the scale factor in simulations with different ε. Bottom panel: length scale R(0.8), including the 80% of the mass of the virialized structure, normalised to 2.8ε.

Open with DEXTER

Instead, for ε ≤ 3.7 × 10-5 we find that, because the system remains stable in comoving coordinates R(0.8)/(2.8ε) is larger than 100 for all values of the scale factor considered, and thus the virial ratio fluctuates around –1 at all times. Therefore we conclude that the relevant characteristic size of the structure remains large enough with respect to the force softening and that in this situation the Newtonian value of the mean field potential is well approximated by the smoothed potential.

The behaviour of the density profile is also coherent with what we have already found in Sect. 4: for a < aLI, the density profile obeys to Eq. (3).

thumbnail Fig. 11

Density profile in physical coordinates at a = 4.6 for the simulations M4a and M1a. The best fit with the NFW and QEP profiles are also reported.

Open with DEXTER

Then, for a > aLI, we observe a breaking of the stable clustering regime, and the density profile rapidly changes shape. In particular, it firstly deviates at small radii, displaying a power-law behaviour, and then at large ones where the slope tends to − 3. The whole density profile can be fitted by a NFW behaviour (see Fig. 11).

The evolution of the velocity PDF f(v) and the nearest neighbours distribution ω(x) for different values of ε are coherent with the behaviours previously discussed. Namely, for large enough scale factors, i.e. a ≥ 10, the simulation which behaves better according to the LI test (i.e., ε = 3.7 × 10-5) shows a more extended tail of f(v) at high velocities than the others. This means a better ability of the code to follow high velocity particles. Correspondingly also ω(x) is peaked at smaller separations i.e. close particles can be better resolved for that value of the softening length. Both these facts point towards a better spatial (i.e., close neighbours) and temporal (i.e., high velocity particles) resolution when ε = 3.7 × 10-5.

Similarly to the behaviour of the density profile, whose small scale properties depend on ε the nearest neighbour distribution and its peak shown an ε dependence implying that the small scale properties of the system are determined by the spatial resolution used in the simulation.

5.2. Time-stepping accuracy

The accuracy of the time-stepping in the GADGET code is controlled by the parameter η: in particular, the time-stepping is fixed by (29)where is the particle acceleration. It thus depends both on the softening length ε and on the parameter η, i.e. the spatial and temporal resolution are strongly entangled. To determine the role of η, we consider a set of simulations in which all physical (N,R0,b0) and numerical parameters (ε, θ, φ) are fixed and the time-stepping accuracy η is varied by a factor 102 (see Table 1 for details).

From the LI test (see Fig. 12) we find aLI = aLI(η). Note that, as previously found, the LI test is verified in about the same range of scale factors where the potential energy is roughly constant and the structure is in the stable clustering regime.

thumbnail Fig. 12

Behaviour of A(a) (Eq. (26)) for simulations with different η and ε (see labels).

Open with DEXTER

thumbnail Fig. 13

Behaviour of aLI as a function of η for different values of the softening length ε.

Open with DEXTER

The behaviour of aLI for different values of ε and η is reported in Fig. 13: the larger is η and the sooner occurs (i) the breakdown of the LI test; (ii) the departure from the stable clustering regime; (iii) the deviation from the QSS with the profile given by Eq. (3). That is, at fixed ε, the larger is η the smaller is aLI. In agreement with the results discussed above, we find that the NFW profile provides a very good fit of the density profile for a > aLI (see Fig. 14).

thumbnail Fig. 14

Density profile in physical coordinates at a = 2.5 respectively for η = 0.0025 and η = 0.25. In both the cases the softening length is ε = 3.7 × 10-5.

Open with DEXTER

The behaviours of the velocity f(v) and of the nearest neighbours ω(x) distributions show that the smaller is η and the better is the resolution of close neighbours and of high velocity particles. In addition, the high velocity tail of f(v) is substantially reduced in the regime a > aLI, in agreement with the results obtained by changing the softening length.

5.3. Other numerical parameters

Results obtained by varying the opening angle θ of the tree algorithm from 0.7 to 0.1 and the accuracy of the relative cell-opening criterion φ, i.e. the parameter denominated ErrTolForceAcc, from 5 × 10-3 to 5 × 10-4 do not show sensible differences. This fact gives us reasonable confidence that the results we have presented are quite well converged with respect to these parameters.

5.4. Summary

The two main parameters that control the spatial and temporal resolution of the numerical integration, namely the softening length ε and the time-stepping accuracy η, determine the range of time in which a simulation is well performing, i.e., the structure is a QSS and Wp ≈ const., the evolution is in agreement with stable clustering and the LI test is reasonably well satisfied. Thus the scale factor at which the breakdown of the LI test occurs is such that aLI = aLI(ε,η). We note that most of the discussion in the literature has focused the attention on the role of the softening length (see e.g. Splinter et al. 1988; Heitmann et al. 2005; Romeo et al. 2008; Joyce et al. 2009a; Power et al. 2003): here we have found that a crucial parameter is also represented by the time-stepping accuracy.

Coherently with these results we found that the velocity PDF f(v) and nearest neighbour distribution ω(x) depends on both the two parameters ε and η. In particular the high velocity tail of f(v) remains stable for a < aLI while there is a sensible deficit of high velocity particles for larger scale factors, due to the limited temporal resolution. Similarly ω(x) does not change anymore for a > aLI, implying the departure from the stable clustering regime. In this situation the peak of ω(x) is determined by the values of both parameters ε,η.

In summary, tests performed by varying the softening parameter show that at fixed ε <  the smaller is η and the larger is aLI. Our interpretation of this behaviour is that at a given level of numerical accuracy, determined by the choice of ε,η, hard scattering between particles can be resolved only up to a certain scale factor. Finally we note that the results of the evolution of the density profile are coherent with the ones found in Sect. 4, namely for a < aLI this is well approximated by the QEP shape, while for a > aLI the profile shows the NFW behaviour.

We have discussed that there is an artificial injection of energy into the system during time evolution. We now conjecture that this is due to an inaccurate integration of two-body scatterings. When an hard scattering occurs the acceleration can be huge for a very short time interval δτ: however, if the time stepping is larger than δτ then the particle velocity will be overestimated by a large factor: this is a known numerical problem discussed in detail by, e.g., Knebe et al. (2000). Indeed, when η decreases the amplitude of the time stepping (see Eq. (29)) gets shorter allowing a better sampling of the change of velocity during hard scatterings.

thumbnail Fig. 15

Behaviour of A(a) (Eq. (26)) for simulations with different N (see Table 1). Insert panel: behaviour of the virial ratio for the same simulations.

Open with DEXTER

Instead, the change of particle velocity is so that for hard scatterings we have : the over-estimation of the velocity change increases when ε gets smaller. Therefore we may conclude that the measured dependence of aLI on ε,η and N is coherent with the interpretation that there is a numerical problem in the integration of hard scatterings. Collisionality for the smaller softenings we use manifests itself most clearly in the breakdown of the LI test that represents a useful tool to study these effects.

It is interesting to note that collisional effects in cosmological simulations have been documented in previous studies by considering the comparison of different codes. In particular, Knebe et al. (2000) very explicitly demonstrates the effects of poor integration of hard two body collisions at small softenings. They showed that two-body scattering have effects which can be greatly amplified by inaccuracies of time integration and that the effects of such poorly integrated hard collisions is to lead to an artificial injection of energy into the system which causes an increasing its size. These results are indeed fully compatible with our results and interpretation (see also the discussion in Joyce & Sylos Labini 2013).

6. Role of the physical parameters

In this section we first study the evolution of systems with same physical and numerical parameters but different particles number. That is, we discretize the system with a different number N of particles taking fixed the total mass M = m × N, with m ∝ 1/N. Then we consider systems of different size R0 but with the number of particles, their mass and the numerical parameters [ε,η,θ,φ] fixed: this corresponds to change only the value of δ in Eq. (12). Finally we change only the initial velocity dispersion and thus the initial virial ratio b0.

6.1. Scaling with particles number

When the number of particles is small enough, i.e. N < N ≈ 104 we observe that during the collapse phase the LI test always fails as A(a) presents relatively large deviations from unity since a ≳ 1 which persist at all times (see Fig. 15). This can be interpreted as due to the fact that a few hard collisions occurring in the collapse phase substantially perturb the integration of the equation of motions. Note that according to the results obtained in Sect. 5 the value of N should in general depend on ε,η.

In simulations with N > N ≈ 104 we find instead that A(a) ≈ 1 soon after the collapse phase, with smaller fluctuations for larger N values. The scale factor signing the break of the LI test weakly depends on N, i.e. aLI = aLI(N): indeed it is found to be almost independent on N. This latter result seems surprising as from Eq. (1) the collisional relaxation time should increase with N so that one would naively expect that larger is N and larger is aLI.

To interpret this behaviour we recall that the density of the core after the collapse, nc in Eq. (3), scales with the number of particles as nc ∝ N2: at fixed mass density, by increasing the number of particles, the density of the core increases as ρc = nc × m ∝ N, given that particles mass scales as m ∝ N-1 (Joyce et al. 2009b). From Eq. (1) we thus find that , as . Given that κ is weakly dependent on N, i.e. κ ∝ 1/log    (N), we thus find that τ2 is almost independent on N, at least for the limited range of N we considered, i.e. N = 104 ÷ 105.

In Fig. 16 it is shown the density profile in physical coordinates for simulations with N = 103,105 particles: note that n(r) has been normalised by recalling the results of Joyce et al. (2009b) where it was shown that the two parameters in Eq. (3) scale as rc ≈ N− 1/3 and nc ≈ N2. As we found in the previous sections, the breakdown of the stable clustering regime corresponds to a deformation of the profile. In the case shown in Fig. 16, while for N = 105 the simulation is still in the QSS regime, for the N = 103 case the density profile is already well fitted by the NFW shape (see Fig. 16).

thumbnail Fig. 16

Density profile in physical coordinates for simulations with 103 and 106 particles, at a = 4.9, renormalised as explained in the text at different values of the expansion factor. The best fit with Eq. (28) (NFW) and Eq. (3) (QEP), respectively for N = 103 and N = 105 are reported as reference.

Open with DEXTER

6.2. Finite size effects

By considering the behaviours of the gravitational potential energy and of the LI test for simulations with different initial size R0 we can conclude that aLI = aLI(R0). In particular, the smaller is R0 and the smaller is aLI.

thumbnail Fig. 17

Behaviour of A(a) (Eq. (26)) for simulations with different initial size R0 as a function of time (normalised to τc(R0)).

Open with DEXTER

This fact can be simply interpreted by considering that when the system size decreases, the amplitude of the over-density δ = δ(R0) increases and so the value of the characteristic time scale τc = τc(R0) (see Table 2). Thus, to compare the behaviours for different R0 we plot (see Fig. 17 – bottom panels) Wp(a) and A(a) as a function of time in units, for each value of R0, of the corresponding value of τc = τc(R0). The residual difference in the observed behaviours is due to the effect of the different ratio R0/L: when R0 > 1/10   L finite size corrections, i.e. the terms neglected in approximating Eq. (7) with Eq. (10) become important (Joyce & Sylos Labini 2013).

The evolution of the density profile in physical coordinates for the systems with different R0 (see Fig. 18) shows that they firstly reach the QEP behaviour; then the density profile approaches the r-3 behaviour, at different values of the scale factor for different R0 and in a coherent way with the behaviour of A(a).

thumbnail Fig. 18

Density profile in physical coordinates at a = 2.5 for different values of R0 (see labels). The vertical dashed line correspond to a × ε. A reference line with slope − 4 is plotted as reference.

Open with DEXTER

Finally we note that for long enough values of the scale factor the density profile for different R0 converges precisely to the same behaviour characterised by a x-3 decay. As discussed in Sect. 5, at fixed time-stepping accuracy, the small scale properties of the distributions are determined by the value of the softening length. Indeed, independently on the initial value of R0, for long enough times, all curves collapse on each other.

thumbnail Fig. 19

Density profile in comoving coordinates at a = 70 for the simulations with ε = 3.7 × 10-4 and a different value of the initial radius R0 (see labels). A reference line with slope − 4 is plotted as reference.

Open with DEXTER

6.3. Role of the initial velocity dispersion

The initial conditions we have considered can be thought to represent an halo which virializes in a relatively short time and which then evolve in isolation from the rest of the mass of the universe. Indeed, when b0 →  − 1 the halo is already almost virialized while when b0 → 0 the halo is far from such a configuration. However, for b0 ∈  [ − 1,0] the structure undergoes to a collapse, which can be less or more violent, that drives the system, in about the same time-scale τc (see Eq. (2)), into a virialized QSS that can be well described as a collision-less equilibrium configuration (Sylos Labini 2012).

Figure 20 shows the behaviour of A(a) for simulations with different initial virial ratio b0 ∈  [ − 1,0], i.e. with different initial velocity dispersion. We find that aLI = aLI(b0): in particular, the smaller is b0 the larger is aLI. This behaviour is easily explained by considering that the closer the system is initially to the virial configuration and the less violent is the dynamics that drives the system into the virialized QSS (Sylos Labini 2012): the larger is b0 at the higher is the density of the core reached after the collapse.

thumbnail Fig. 20

Behaviour of A(a) (Eq. (26)) for for simulations with different initial size b0.

Open with DEXTER

We recall that, in the STAT case, the formation of the r-4 is observed only when the collapse is violent enough that a fraction of the initial mass and energy is ejected from the system at the formation of the QSS: in Sylos Labini (2012) it was found that this occurs for b0 >  − 1/2. In this situation the system evolution is very similar to the b0 = 0 case: at the beginning the density profile is well approximated by the QEP and then, for a > aLI it develops a r-3 tail approaching the NFW shape also at small radii. Thus, although the b0 = 0 is peculiar because the collapse is instantaneous in the continuous limit, it does not represent a pathological case and the r-4 tail is formed as long as the collapse is violent enough that part of the mass and energy is ejected from the virialized structure.

For values of the initial virial ratio b0 <  − 1/2 the virialized structure after the collapse does not present a density profile with a 1/r4 tail. However, for long enough times the density profile approaches the 1/r3 behaviour for any value of b0. This shows that the 1/r3 is formed through the dynamical process and that it is thus independent on the initial conditions.

In Fig. 21 it is shown the comparison of the density profile, in comoving coordinates, for simulations with different initial virial ratio and at different expansion factors: this confirms that the smaller is b0 and the slower is the evolution. In particular, the smaller b0 the longer it takes to form the 1/r3 tail, a behaviour that is approached at long enough times for any value of b0.

thumbnail Fig. 21

Comparison of the density profile in comoving coordinates for simulation with different initial virial ratio a a = 70.

Open with DEXTER

6.4. Summary

In summary we found that aLI = aLI(R0,N,b0), i.e. the regime in which the system evolution is correctly numerically integrated depends on the physical parameters of the initial conditions.

For what concerns the dependence of aLI on N, we have noticed that: (i) for relatively small particles number, i.e. N < N ≈ 104, the system never reaches the QSS in the stable clustering regime. We have interpreted this behaviour as due the fact that a few hard scatterings are enough to perturb the macroscopic dynamics. (ii) On the other hand for large enough particles number, i.e., N ∈  [104,105], we find that the deviation from the QSS slowly depends on N. This occurs because when initial conditions are cold, the system contracts more for larger N values, so that the virialized state has a core mass density ρc ∝ N and thus the typical time scale for crossing is : in this situation the collisional time scale τ2 (see Eq. (1)) is very weakly dependenton N.

By decreasing the system size R0 we find that aLI decreases as well: this is coherent with the fact that the smaller is R0 and the smaller is the average distance between nearest neighbour particles in the core, and thus the higher the probability of hard collisions. However, when we renormalize the behaviours to the proper characteristic time τc = τc(R0) we find that the differences in the departure from the stable clustering regime become smaller: these residual differences are due to finite size effects.

In addition we have found that for long enough times all systems with initially different R0 converge to the same density profile characterised by a r-3 decay. This convergence is obtained because the numerical integration is badly performed, as shown by the LI test.

Finally, by increasing the initial velocity dispersion, the breakdown scale factor aLI is found to increase. This behaviour can be simply explained by considering that when b0 → 0 the collapse is much more violent than for b0 →  − 1; thus in the former case the core may reach higher densities sooner, then having a similar evolution to the cold case. In addition, when the initial velocity dispersion is high enough, the system remains in the stable clustering regime longer and until the softening length, that increases in physical coordinates, becomes large enough to perturb the evolution. At later times the system shows an evolution that is again similar to the cold case.

7. Discussion and conclusion

In this paper we have tested a cosmological code (GADGET) by choosing very simple and idealised initial conditions. While we do not claim that this choice can be justified by any physical argument, i.e. the toy model does not represent any realistic physical system, we stress that in this way we have been ableto make a series of controlled tests on the accuracy of a cosmological code. Given what we have learnt from this case, we discuss here some general ideas for the testing of a cosmological code in the case of more realistic initial conditions, as those predicted by standard models of galaxy formation. To consider separately these two different but, as we argue below, probably related issues, we first discuss the results of the simulations we have presented in this paper and then we briefly comment on the more complicated problem that requires more tests and work to be properly explored.

7.1. The control of the accuracy of a cosmological code: a toy model

The evolution of an isolated over-density in a cosmological N-body simulation is approximately equivalent, in physical coordinates, to the evolution of the same structure in an open system without expansion (Joyce & Sylos Labini 2013). This equivalence is strictly verified in the limit of vanishing force smoothing and when the size of the over-density is much smaller than the size of the simulating box. When these conditions are satisfied, the evolution of the open system can be used as a template to compare the evolution in an expanding background with the aim of determining the accuracy of a cosmological simulation. The open system evolution is such that it relaxes to a virialized QSS in a time τc (see Eq. (2)) (van Albada 1982; Aarseth et al. 1988; Joyce et al. 2009b; Sylos Labini 2012) and it evaporates, due to collisional effects, in a time τ2 ≫ τc (see Eq. (1)). The QSS corresponds, in a cosmological N-body simulation, to an isolated structure in the stable clustering regime.

Joyce & Sylos Labini (2013) considered the evolution of an isolated over-density initially close to virial equilibrium: for this reason during the collapse its size and density vary a little. In this paper we considered a more general case, namely an isolated over-density with initial virial ratio in the range b0 −  ∈  [ − 1,0]. When b0 →  − 1 the initial conditions are warm enough that the collapse consists in a relatively small modification of system size. Instead when initial velocities are cold, i.e. b0 → 0, the system largely contracts during the collapse, reducing considerably its size and shape and thus rapidly reaching, in its central core, a very high density (Sylos Labini 2012).

Note that an initially cold, i.e. b0 = 0, and spherical over-density does not represent a pathological case from the dynamical point of view: indeed, as firstly noticed by Stiavelli & Bertin (1987), a density profile with a r-4 tail is a generic property of systems formed by a violent relaxation process. For instance, recently Sylos Labini (2012, 2013) found that this profile is formed also in the collapse of an initially spherical and isolated cloud of particles with a small enough velocity dispersion, with a uniform or power-law density profiles.

We found that, for b0 → 0, once the system has reached a virialized configuration, its evolution can be characterised by two different regimes:

  • (i)

    for scale factors a < aLI the potential energy in physical coordinates Wp isapproximately constant, the Layzer Irvine (LI) equation,describing the energy change in an expanding background, iswell satisfied, the density profile is characterised by a r-4 tail,characteristic of the QEP in the open case (seeEq. (3)), and the system is in the stable clusteringregime, which means that it remains stable in physicalcoordinates.

  • (ii)

    On the other hand for a > aLI we find that | Wp |  ~ a-1, there is a breakdown of the LI equation, the density profile becomes stable in comoving coordinates, which corresponds to the breakdown of the stable clustering regime, it develops a r-3 tail and it can reasonably be fitted by the NFW profile.

In principle, the QSS has a lifetime τ2 that diverges with the number N of system particles because of collisional effects. However, we found that the time scale corresponding to aLI is (i) tLI < τ2 and (ii) it depends on the numerical parameters of the simulation. Therefore this time scale is a numerical artifact: the relevant question is what determines it.

To answer to this question, we have performed various series of simulations by varying the integration parameters, such as the softening length of the gravitational force ε (which, in the simulations we have performed, is fixed in comoving coordinates during the evolution) and the accuracy of the time-stepping η. We found that In particular, we found that there is a narrow range of ε (at fixed η) for which the simulation behaves “well” (i.e. regime (i) above). When ε is much smaller than the initial inter-particle distance the integration is more difficult and aLI decreases: this fact can be interpreted as due the effect of hard collisions which require a high numerical accuracy to be correctly integrated. On the other hand, when ε is of the order of the characteristic scale of the virialized system the mean field potential is no longer Newtonian so that the system deviates from the virialized configuration and from the QEP.

That there is a problem in the accuracy of the numerical integration is confirmed by varying η at fixed ε (with ε < ): in particular when η decreases aLI increases. This fact can be again interpreted in terms of a better accuracy of the numerical integration and in particular a more accurate integration of high velocity particles, which require a smaller time-stepping to be properly followed. Both tests by varying ε and η are thus consistent with the interpretation outlined above, that confirms the results by Joyce & Sylos Labini (2013): two-body scatterings perturb the accuracy of the simulation evolution. Thus the temporal and spatial resolution of the simulation are entangled and both must be considered for the determination of the optimal values of the numerical parameters.

Note that codes with adaptive spatial resolution (see e.g., Knebe et al. 2000; Teyssier 2001; Miniati & Colella 2007) may perform much better than the fixed resolution case which we have considered in this paper. Indeed, it is clear such codes will, by construction, avoid respect. The tests we have discussed in this paper can be easily applied to other codes and a detailed study of codes with an adaptive mesh refinement is postponed to a forthcoming work.

Because of the violation of the LI, for a > aLI, the system expands in physical coordinates, thus increasing its potential energy: this is due to an injection of energy into the system. From the tests performed by varying the softening length and the time-stepping parameter we conclude it is the effects of poor integrated hard two body collisions that leads to an artificial injection of energy into the system. This same conclusion was drawn by Knebe et al. (2000) by considering a comparison of different N-body codes. In summary we find, in agreement with Joyce & Sylos Labini (2013), that given a certain choice of numerical parameters, one may follow structure formation in the non-linear regime, for a limited time which depends on several physical properties of the system. The results obtained here clearly show the role of the numerical (in)accuracy in the determination of the shape of the density profile structure.

In order to further test the result that the numerical integration becomes inaccurate because of collisional effects, we have, at fixed numerical parameters, changed the physical parameters of the initial conditions. Namely we have varied the number of particles N (at fixed density), the size of the structure R0 and the initial velocity dispersion (and thus the initial virial ratio b0). As result we found that When decreasing R0 we find that aLI decreases: the initial inter-particle distance gets smaller, as  ∝ R0, and thus hard scatterings are more probable. The dependence of aLI on the particle number clearly shows the effect of two-body collisions. For small particle number, i.e. N < 104, the system never reaches a QSS and the LI test fails at all times. We then find that, at fixed ε and η, and for N ∈  [104,105] the breakdown expansion factor aLI is almost independent with N. This behaviour can be interpreted as due to the fact that the larger is N, the larger the core mass density and thus the characteristic time τc: in this situation the collisional time scale τ2 in Eq. (1) becomes almost independent on N.

Finally by increasing the initial velocity dispersion, i.e. decreasing b0, we find that aLI increases: indeed, as shown in Sylos Labini (2012), when the initial configuration is closer to the virial equilibrium, i.e. b0 =  −1, the collapse is less violent so that the size and density of the structure after the collapse are not largely changed as it occurs when b0 = 0. In this latter case the core of the structure after the collapse may reach very high densities, i.e. a more favourable condition for the occurrence of two-body scatterings. However, we have found that for b0 ∈  [ − 1,0] the density profile, in the long time evolution, approaches the same NFW shape. Thus, the density profile, for long enough times, approaches the 1/r3 behaviour for any value of b0 a fact that shows that this tail is formed by the same dynamics process, independent on the initial conditions, in all cases.

While the origin of the r-4 can be understood by considering that the energy of particles orbiting around the core is conserved Sylos Labini (2012), the effect of inaccurately integrated two-body scatterings is indeed that of perturbing such energy conservation. The energy injection, is especially efficient for the fastest particles, i.e. those that form the tail of the density profile. A quantitative study of this effect is postponed to a forthcoming work, but qualitatively we may conclude that the deformation of the density profile is related to the inability to resolve numerically high velocity particles and thus to conserve their energy properly.

7.2. Some comments on more complex initial conditions

The simplicity of the system considered as initial conditions has let us to develop systematic tests to control the accuracy of a cosmological code. The principal result that we have obtained is that the range of redshifts where the numerical integration of a certain structure is reliable, i.e. it follows the stable clustering regime and it passes the LI test, is in general limited to narrow one that is dependent on the various physical (number of particles, velocity dispersion, density and amplitude of the over-density) and numerical (softening length, time-stepping) parameters of a given simulation. Therefore, our results suggest that the limited numerical resolution of a given simulation can be not sufficient to integrate the formation of structures with a wide range in size and mass as those formed in the currently favoured cosmological models structures through a hierarchical aggregation. A more detailed study of the more complex cosmological case is thus required.

Here we note that, if one is not able to integrate correctly the evolution of this simplified system, one will encounter more difficult problems if one wants to properly simulate the formation of structures of different sizes, with different particle numbers and in a complex environment as that occurring in the case of typical initial conditions of structure formation models. On the other hand the choice of the physical and numerical parameters adopted in this paper could be completely unrealistic and thus irrelevant for the case of a cosmological simulations that uses initial conditions as those predicted by standard models of galaxy formation.

In this respect, we note that the numerical parameters have been varied in a relatively large range, i.e. ε ∈  [10-6,10-3] and η ∈ 2.5 ×  [10-3,10-1]. In addition, the number of particles used, ranging in N ∈  [103,105], is comparable to the number of particles which form a single halo in a typical cosmological simulation (Springel et al. 2005). Thus the system we considered, once it reaches a virialized state, can be though to represent a simplified toy model approximating a virialized cosmological halo at its formation time. However, this structure then evolves as if it were isolated from the rest of the universe. Thus, such a self-gravitating system presents several differences with respect to the haloes formed in typical cosmological simulations: it is spherical and isolated so that it has initially a simple geometry, it is not subjected neither to mass accretion nor to tidal interactions with the neighbouring structures, it does not fragment into (macroscopic) substructures during its evolution and it does not collide with structures of comparable size.

A possible method to control both the role of the physical and the numerical parameters in a more complex cosmological simulations was presented in Joyce & Sylos Labini (2013), where it was outlined a sort of recipe that can be simply implemented and that can furnish a general manner for testing the accuracy of a simulation of haloes containing N particles. This is simply based on the testing of whether an halo formed in a complex environment, once it is isolated, evolves in agreement with stable clustering and represents a quasi-stationary configuration, given the same numerical parameter of the numerical integration (softening length, time-stepping, etc.) and the same physical parameters of the structure (number of particles, size and velocity dispersion). The implementation of this test is postponed to a forthcoming work.

As a final remark, we stress that we have observed that the density profile of the virialized structure formed from the simple initial conditions we have considered converges, for a wide range the numerical and physical parameters of the simulation, to the NFW profile, i.e. the same profile observed in cosmological N-body simulations that start from the typical initial conditions predicted by standard models of galaxy formation. Whether this is a coincidence or whether also standard cosmological haloes are seriously affected by the numerical (in)accuracy of the numerical integration is an open problem that will explored in full details in forthcoming works


2

The value of ε is the value of the parameter with this name in the code.

3

But θ = 0.1 in the simulation M4b.

4

But φ = 5 × 10-4 in the simulation M4c.

5

I.e., the parameter denominated ErrTolForceAcc in the code.

6

Recall that ε is given in units of the box side L.

Acknowledgments

I thank Michael Joyce for useful discussions and comments, Roberto Ammendola and Nazario Tantalo for the valuable assistance in the use of the Fermi supercomputer where the simulations have been performed. I also would like to thank Alessandro Romeo for a number of very useful comments and suggestions that have allowed to improve the presentation of the paper.

References

All Tables

Table 1

Details of the N-body simulations.

Table 2

Parameters of the spherical clouds with different initial radius R0 and with over-density δ (Eq. (12)).

All Figures

thumbnail Fig. 1

Upper panel: energy in physical coordinates (Kp Eq. (15) and Wp Eq. (16)) normalised to the initial total energy as a function of time (computed by using Eq. (14) for the EXP case), together with the corresponding behaviour of the same system in the STAT case. Bottom panel: virial ratio | b(t) |, together with the LI test (Eq. (26)) for the EXP case.

Open with DEXTER
In the text
thumbnail Fig. 2

Comparison of the density (upper panel) and velocity (bottom panel) profiles of the same system evolved in an expanding background (EXP at a = 1.36 – black circles) and in a static background with open boundary conditions (STAT at t = 5τc – red squares). The best fit with Eq. (3) (QEP) and a line with slope − 1 are shown as references (see Eq. (4)).

Open with DEXTER
In the text
thumbnail Fig. 3

Upper panel: absolute value of the potential energy in physical coordinates Wp(a) (Eq. (16)) as a function of the scale factor. The approximate value of aLI, separating the regime where Wp ≈ const. from the one where | Wp |  ~ a-1, is reported as reference. Bottom panel: LI test (Eq. (26)) and virial ratio (absolute value) for particles with negative energy.

Open with DEXTER
In the text
thumbnail Fig. 4

Evolution of the density profile in comoving coordinates. Two reference lines with slope –4 and –3 are reported as a reference. The arrow shows the direction of the density profile evolution. Inset panel: behaviours in function of the scale factor of the length scale R(0.4) (R(0.8)), including the 40% (80%)of the mass of the virialized structure, normalised to 2.8ε.

Open with DEXTER
In the text
thumbnail Fig. 5

Density profile in physical coordinates for the simulation M4a at different scale factors a (see labels) together with the density profile for the STAT case. The arrow shows the direction of the density profile evolution in the cosmological case. Further it is shown the best fit given by Eq. (3) (QEP) and a reference line with slope –3.

Open with DEXTER
In the text
thumbnail Fig. 6

Radial velocity dispersion profile in physical coordinates at different scale factors a (see labels) together with the profile for the STAT case. The arrow shows the direction of the density profile evolution in the cosmological case. A line with slope r-1 (see Eq. (4)) is reported as reference.

Open with DEXTER
In the text
thumbnail Fig. 7

Evolution of the probability distribution function (PDF) of the absolute value of the velocity at different scale factors (see labels). The arrow shows the direction of the evolution.

Open with DEXTER
In the text
thumbnail Fig. 8

Evolution of the nearest neighbour distribution in comoving coordinates at different scale factors (see labels). The arrow shows the direction of the evolution. The initial distribution (IC) is also shown for comparison.

Open with DEXTER
In the text
thumbnail Fig. 9

Behaviour of A(a) (Eq. (26)) for simulations with different ε (see Table 1 for details).

Open with DEXTER
In the text
thumbnail Fig. 10

Upper panel: virial ratio as a function of the scale factor in simulations with different ε. Bottom panel: length scale R(0.8), including the 80% of the mass of the virialized structure, normalised to 2.8ε.

Open with DEXTER
In the text
thumbnail Fig. 11

Density profile in physical coordinates at a = 4.6 for the simulations M4a and M1a. The best fit with the NFW and QEP profiles are also reported.

Open with DEXTER
In the text
thumbnail Fig. 12

Behaviour of A(a) (Eq. (26)) for simulations with different η and ε (see labels).

Open with DEXTER
In the text
thumbnail Fig. 13

Behaviour of aLI as a function of η for different values of the softening length ε.

Open with DEXTER
In the text
thumbnail Fig. 14

Density profile in physical coordinates at a = 2.5 respectively for η = 0.0025 and η = 0.25. In both the cases the softening length is ε = 3.7 × 10-5.

Open with DEXTER
In the text
thumbnail Fig. 15

Behaviour of A(a) (Eq. (26)) for simulations with different N (see Table 1). Insert panel: behaviour of the virial ratio for the same simulations.

Open with DEXTER
In the text
thumbnail Fig. 16

Density profile in physical coordinates for simulations with 103 and 106 particles, at a = 4.9, renormalised as explained in the text at different values of the expansion factor. The best fit with Eq. (28) (NFW) and Eq. (3) (QEP), respectively for N = 103 and N = 105 are reported as reference.

Open with DEXTER
In the text
thumbnail Fig. 17

Behaviour of A(a) (Eq. (26)) for simulations with different initial size R0 as a function of time (normalised to τc(R0)).

Open with DEXTER
In the text
thumbnail Fig. 18

Density profile in physical coordinates at a = 2.5 for different values of R0 (see labels). The vertical dashed line correspond to a × ε. A reference line with slope − 4 is plotted as reference.

Open with DEXTER
In the text
thumbnail Fig. 19

Density profile in comoving coordinates at a = 70 for the simulations with ε = 3.7 × 10-4 and a different value of the initial radius R0 (see labels). A reference line with slope − 4 is plotted as reference.

Open with DEXTER
In the text
thumbnail Fig. 20

Behaviour of A(a) (Eq. (26)) for for simulations with different initial size b0.

Open with DEXTER
In the text
thumbnail Fig. 21

Comparison of the density profile in comoving coordinates for simulation with different initial virial ratio a a = 70.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.