© 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 R₀/L, where R₀ 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 b₀ ∈  [−1,0]. When the system is warm enough, i.e. b₀ ≈  −1, it approximately maintains its original size and shape during and after the collapse. Instead, when the system is cold enough, i.e. b₀ ≈ 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 b₀ = 2K₀/W₀ ≈  − 1 (where K₀/W₀ 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 R₀ 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 r_c ≈ 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 r_c and n_c 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 v_r 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 R₀) 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 x_i 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 (R₀/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 r_i ≡ a(t)x_i, 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 R₀/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 R₀ < 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 10³ ≤ N < 10⁴ (series S), N = 10⁴ (series M) and 10⁴ < N ≤ 10⁵ (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 W₀ = 3   GM²/(5   R₀) =  −1, i.e. in which the initial continuum limit gravitational potential energy is minus unity.

Note that the system radius R₀ is chosen to be R₀ < 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 ρ₀ = M/L³ is the mean mass density of the universe (see Table 2). Particles mass is such that ρ₀ is equal to the critical density.

The scale factor in an EdS universe obeys to (Peebles 1980) (13)where . Assuming a₀ = 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. δ > 10⁴, at moderate expansion factors a ≈ 10, (t − t₀)/τ_c is of the order of a thousand i.e. of the order of τ₂ for N = 10³ (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-2¹ 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 vanishes². 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 thereafter³. The accuracy of the relative cell-opening criterion (i.e., the parameter ErrTolForceAcc)⁴ 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 b₀ = 2K₀/W₀, where K₀,W₀ 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 R₀ ≪ L, we have W_g ≈ W_c.

The peculiar kinetic energy is (21)where (22)If the system is virialized in a region of size R₀ ≪ 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. K_p ≈ K_c.

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 a_LI such that (27)and larger than 0.1 for a > a_LI. 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 W₀). 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. W_p ≈ const., and it satisfies the virial equilibrium (Fig. 1 – bottom panel). A fraction of about f_p ≈ 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 < a_LI (see Eq. (27)).

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.

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 W_p ≈ const., i.e. the structure is in a QSS and in the stable clustering regime, only for a < a_LI ≈ 6 (see Fig. 3). Instead, for a > a_LI we find that | W_p |  ~ 1/a (and thus W_c ≈ const). Thus the breakdown of the LI test corresponds to a departure from the stable clustering regime.

Indeed, the density profile in comoving coordinates for a > a_LI 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 = aε. When the structure is in the stable clustering regime, its size in physical coordinates R_p is constant, and thus R_p/ε_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 > a_LI 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 > a_LI, 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 < a_LI, 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 e_p ≲ 0 may go asymptotically far from the structure core (Sylos Labini 2012).

On the other hand, for a > a_LI, 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 > a_LI (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.

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 < a_LI the different curves collapse on each other, implying that the virialized structure is in the stable clustering regime, for a > a_LI 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 W_p increases with the scale factor for a > a_LI (as W_c ≈ const.): as time passes and thus W_p increases proportionally to − 1/a.

Thus for a > a_LI 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 =  −ΔK_p. 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 < a_LI we find that f(v) remains almost stable, for a > a_LI 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 > a_LI. 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.

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 < a_LIω(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 > a_LIω(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 < a_LI and for r > r_c, scales as Φ(r) ~ r^− 5/2, i.e., the behaviour expected for the QEP (see Eq. (5)). On the other hand, for a > a_LI, at large radii, we find Φ(r) ~ r^-3 × r^3/2 ~ r^-1.5.

4.3. Summary

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

• 1.

  For a < a_LI the structure is in a QSS (i.e., | W_p |  ≈ 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 > a_LI the structure is (almost) frozen in comoving coordinates and therefore its size in physical coordinates increases and its potential energy tends to zero as W_p ∝  − 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 > a_LI 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 > a_LI, as it should if stable clustering were satisfied. We thus conclude that the inaccurate integration of hard collisions induces the change of regime at a_LI.

Note that a_LI ≈ 6 corresponds to a time (see Eq. (14)) (t_break − t₀)/τ_c ≈ 50 that is much smaller than the collisional relaxation time scale (see Eq. (1)) τ₂ ≈ 10⁴τ_c for a system with N = 10⁴ 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.1. Softening length

The softening length⁶ has been varied in a range such that: (i) ε < R₀, 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 W_p and the IL test (Fig. 9) for simulations of the same system with N = 10⁴ particles and with different values of ε (see Table 1) – changed by a factor 10⁴ – show behaviours coherent with those discussed in Sect. 4.2. Namely the two regimes found above, corresponding to W_p ≈ const. for a < a_LI and | W_p |  ~ a^-1 for a > a_LI, still characterise the temporal evolution. However, the transition between them depends on the softening length, i.e., a_LI = a_LI(ε).

In particular, we find that a_LI ≈ 6 for ε ∈  [3.7 × 10^-5,3.7 × 10^-4], while outside this range a_LI ≈ 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.

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

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 < a_LI, the density profile obeys to Eq. (3).

Then, for a > a_LI, 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,R₀,b₀) and numerical parameters (ε, θ, φ) are fixed and the time-stepping accuracy η is varied by a factor 10² (see Table 1 for details).

From the LI test (see Fig. 12) we find a_LI = a_LI(η). 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.

The behaviour of a_LI 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 a_LI. In agreement with the results discussed above, we find that the NFW profile provides a very good fit of the density profile for a > a_LI (see Fig. 14).

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 > a_LI, 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 W_p ≈ 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 a_LI = a_LI(ε,η). 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 < a_LI 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 > a_LI, 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 a_LI. 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 < a_LI this is well approximated by the QEP shape, while for a > a_LI 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.

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 a_LI 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 R₀ 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 b₀.

6.1. Scaling with particles number

When the number of particles is small enough, i.e. N < N^∗ ≈ 10⁴ 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^∗ ≈ 10⁴ 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. a_LI = a_LI(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 a_LI.

To interpret this behaviour we recall that the density of the core after the collapse, n_c in Eq. (3), scales with the number of particles as n_c ∝ N²: at fixed mass density, by increasing the number of particles, the density of the core increases as ρ_c = n_c × 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 τ₂ is almost independent on N, at least for the limited range of N we considered, i.e. N = 10⁴ ÷ 10⁵.

In Fig. 16 it is shown the density profile in physical coordinates for simulations with N = 10³,10⁵ 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 r_c ≈ N^− 1/3 and n_c ≈ N². 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 = 10⁵ the simulation is still in the QSS regime, for the N = 10³ case the density profile is already well fitted by the NFW shape (see Fig. 16).

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 R₀ we can conclude that a_LI = a_LI(R₀). In particular, the smaller is R₀ and the smaller is a_LI.

This fact can be simply interpreted by considering that when the system size decreases, the amplitude of the over-density δ = δ(R₀) increases and so the value of the characteristic time scale τ_c = τ_c(R₀) (see Table 2). Thus, to compare the behaviours for different R₀ we plot (see Fig. 17 – bottom panels) W_p(a) and A(a) as a function of time in units, for each value of R₀, of the corresponding value of τ_c = τ_c(R₀). The residual difference in the observed behaviours is due to the effect of the different ratio R₀/L: when R₀ > 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 R₀ (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 R₀ and in a coherent way with the behaviour of A(a).

Finally we note that for long enough values of the scale factor the density profile for different R₀ 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 R₀, for long enough times, all curves collapse on each other.

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 b₀ →  − 1 the halo is already almost virialized while when b₀ → 0 the halo is far from such a configuration. However, for b₀ ∈  [ − 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 b₀ ∈  [ − 1,0], i.e. with different initial velocity dispersion. We find that a_LI = a_LI(b₀): in particular, the smaller is b₀ the larger is a_LI. 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 b₀ at the higher is the density of the core reached after the collapse.

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 b₀ >  − 1/2. In this situation the system evolution is very similar to the b₀ = 0 case: at the beginning the density profile is well approximated by the QEP and then, for a > a_LI it develops a r^-3 tail approaching the NFW shape also at small radii. Thus, although the b₀ = 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 b₀ <  − 1/2 the virialized structure after the collapse does not present a density profile with a 1/r⁴ tail. However, for long enough times the density profile approaches the 1/r³ behaviour for any value of b₀. This shows that the 1/r³ 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 b₀ and the slower is the evolution. In particular, the smaller b₀ the longer it takes to form the 1/r³ tail, a behaviour that is approached at long enough times for any value of b₀.

6.4. Summary

In summary we found that a_LI = a_LI(R₀,N,b₀), 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 a_LI on N, we have noticed that: (i) for relatively small particles number, i.e. N < N^∗ ≈ 10⁴, 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 ∈  [10⁴,10⁵], 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 τ₂ (see Eq. (1)) is very weakly dependenton N.

By decreasing the system size R₀ we find that a_LI decreases as well: this is coherent with the fact that the smaller is R₀ 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(R₀) 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 R₀ 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 a_LI is found to increase. This behaviour can be simply explained by considering that when b₀ → 0 the collapse is much more violent than for b₀ →  − 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 τ₂ ≫ τ_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 b₀ −  ∈  [ − 1,0]. When b₀ →  − 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. b₀ → 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. b₀ = 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 b₀ → 0, once the system has reached a virialized configuration, its evolution can be characterised by two different regimes:

• (i)

  for scale factors a < a_LI the potential energy in physical coordinates W_p 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 > a_LI we find that | W_p |  ~ 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 τ₂ that diverges with the number N of system particles because of collisional effects. However, we found that the time scale corresponding to a_LI is (i) t_LI < τ₂ 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 a_LI 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 a_LI 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 > a_LI, 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 R₀ and the initial velocity dispersion (and thus the initial virial ratio b₀). As result we found that When decreasing R₀ we find that a_LI decreases: the initial inter-particle distance gets smaller, as ℓ ∝ R₀, and thus hard scatterings are more probable. The dependence of a_LI on the particle number clearly shows the effect of two-body collisions. For small particle number, i.e. N < 10⁴, the system never reaches a QSS and the LI test fails at all times. We then find that, at fixed ε and η, and for N ∈  [10⁴,10⁵] the breakdown expansion factor a_LI 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 τ₂ in Eq. (1) becomes almost independent on N.

Finally by increasing the initial velocity dispersion, i.e. decreasing b₀, we find that a_LI increases: indeed, as shown in Sylos Labini (2012), when the initial configuration is closer to the virial equilibrium, i.e. b₀ =  −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 b₀ = 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 b₀ ∈  [ − 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/r³ behaviour for any value of b₀ 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 ∈  [10³,10⁵], 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

---

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

All Figures