A toy model to test the accuracy of cosmological Nbody simulations
^{1} Centro Studi e Ricerche Enrico Fermi, via Panisperna 89 A, Compendio del Viminale, 00184 Rome, Italy
email: Francesco.SylosLabini@roma1.infn.it
^{2} Istituto dei Sistemi Complessi CNR, via dei Taurini 19, 00185 Rome, Italy
Received: 3 January 2013
Accepted: 31 January 2013
The evolution of an isolated overdensity represents a useful toy model to test the accuracy of a cosmological Nbody code in the nonlinear regime as it is approximately equivalent to that of a truly isolated cloud of particles, with same density profile and velocity distribution, in a nonexpanding background. This is the case as long as the system size is smaller than the simulation box side, so that its interaction with the infinite copies can be neglected. In such a situation, the overdensity rapidly undergoes to a global collapse forming a quasi stationary state in virial equilibrium. However, by evolving the system with a cosmological code (GADGET) for a sufficiently long time, a clear deviation from such quasiequilibrium configuration is observed. This occurs in a time t_{LI} that depends on the values of the simulation numerical parameters such as the softening length and the timestepping accuracy, i.e. it is a numerical artifact related to the limited spatial and temporal resolutions. The analysis of the LayzerIrvine cosmic energy equation confirms that this deviation corresponds to an unphysical dynamical regime. By varying the numerical parameters of the simulation and the physical parameters of the system we show that the unphysical behaviour originates from badly integrated close scatterings of highvelocity particles. We find that, while the structure may remain virialized in the unphysical regime, its density and velocity profiles are modified with respect to the quasiequilibrium configurations, converging, however, to well defined shapes, the former characterised by a Navarro Frenk Whitetype behaviour.
Key words: methods: numerical / galaxies: halos / galaxies: formation / largescale structure of Universe / dark matter
© ESO, 2013
1. Introduction
To control the numerical integration accuracy of a finite selfgravitating 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 LayzerIrvine (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 Nbody 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 nonlinear 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 nonlinear 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 overdensity 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 Nbody 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 quasistationary state (QSS), i.e. a collisionless configuration that is not in thermodynamical equilibrium but that has a finite lifetime fixed by the collisional relaxation timescale (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 overdensity in a cosmological simulation should in principle evolve, in physical coordinates, just as if it were in a nonexpanding 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 nonlinear regime (Peebles 1980), in the case of a single overdensity 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 overdensity 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_{0}/L, where R_{0} is the overdensity 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 overdensity 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 twobody 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_{0} ∈ [−1,0]. When the system is warm enough, i.e. b_{0} ≈ −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} ≈ 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 nonlinear 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 overdensity 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 overdensity 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 timestepping 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 overdensity in a cosmological and a static background
Firstly we briefly recall the main features of the collapse of an isolated overdensity 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 overdensities in a static background
In the STAT case, the global collapse of an isolated and spherical cloud of selfgravitating 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_{0} = 2K_{0}/W_{0} ≈ − 1 (where K_{0}/W_{0} 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_{0} 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 socalled pseudo phasespace 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_{0}) 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.
Dissipationless cosmological Nbody 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 nonzero 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_{0}/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 Einsteinde 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 selfgravitating particles. Therefore, up to finite size corrections that vanish in the limit R_{0}/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 selfgravitating 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 nonlocal 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_{0} < L. Only in this situation the evolution of a single overdensity in a periodic system, i.e., an infinite number of overdensities one in each of the infinite boxes of the periodic system, in an expanding background is equivalent to the evolution of a truly isolated overdensity 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^{3} ≤ N < 10^{4} (series S), N = 10^{4} (series M) and 10^{4} < N ≤ 10^{5} (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_{0} = 3 GM^{2}/(5 R_{0}) = −1, i.e. in which the initial continuum limit gravitational potential energy is minus unity.
Details of the Nbody simulations.
Note that the system radius R_{0} is chosen to be R_{0} < 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 overdensity with initial mass density contrast (12)where is the cloud density and ρ_{0} = M/L^{3} is the mean mass density of the universe (see Table 2). Particles mass is such that ρ_{0} is equal to the critical density.
The scale factor in an EdS universe obeys to (Peebles 1980) (13)where . Assuming a_{0} = 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 overdensity, i.e. δ > 10^{4}, at moderate expansion factors a ≈ 10, (t − t_{0})/τ_{c} is of the order of a thousand i.e. of the order of τ_{2} for N = 10^{3} (see Eq. (2)).
3.2. Numerical integration parameters
We have used the publicly available treecode 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 GADGET2^{1} 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^{2}. 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 BarnesHut opening criterion with an opening angle of the tree algorithm θ = 0.7 for the first force computation and a dynamical updating criterion thereafter^{3}. The accuracy of the relative cellopening criterion (i.e., the parameter ErrTolForceAcc)^{4} is set to φ = 5 × 10^{3}. Finally, we have chosen the timestepping 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_{0} = 2K_{0}/W_{0}, where K_{0},W_{0} 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 twobody potential is calculated by using periodic boundary conditions. As mentioned above, in the limit R_{0} ≪ 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_{0} ≪ 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 LayzerIrvine 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 overdensity 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_{0}). 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)).
Fig. 1 Upper panel: energy in physical coordinates (K_{p} Eq. (15) and W_{p} 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 phasespace density obeys to Eq. (5), i.e. it decays as Φ(r) ~ r^{− 5/2}.
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 shorttime 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.
Fig. 3 Upper panel: absolute value of the potential energy in physical coordinates W_{p}(a) (Eq. (16)) as a function of the scale factor. The approximate value of a_{LI}, separating the regime where W_{p} ≈ const. from the one where  W_{p}  ~ 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 
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 > 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 welldefined 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.
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 < 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.
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 > 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.
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 < 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 phasespace 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 collisionless 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 powerlaw 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 twobody 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}.
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 a_{LI} ≈ 6 corresponds to a time (see Eq. (14)) (t_{break} − t_{0})/τ_{c} ≈ 50 that is much smaller than the collisional relaxation time scale (see Eq. (1)) τ_{2} ≈ 10^{4}τ_{c} for a system with N = 10^{4} particles. This implies that the breakdown of the QSS is not related to the expected evaporation due to twobody 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 timestepping accuracy η, the opening angle θ of the tree algorithm and the accuracy of the relative cellopening criterion φ^{5}.
5.1. Softening length
The softening length^{6} has been varied in a range such that: (i) ε < R_{0}, 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^{4} particles and with different values of ε (see Table 1) – changed by a factor 10^{4} – 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 Nbody method.
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).
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 < a_{LI}, the density profile obeys to Eq. (3).
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 > 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 powerlaw 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. Timestepping accuracy
The accuracy of the timestepping in the GADGET code is controlled by the parameter η: in particular, the timestepping 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_{0},b_{0}) and numerical parameters (ε, θ, φ) are fixed and the timestepping accuracy η is varied by a factor 10^{2} (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.
Fig. 12 Behaviour of A(a) (Eq. (26)) for simulations with different η and ε (see labels). 

Open with DEXTER 
Fig. 13 Behaviour of a_{LI} as a function of η for different values of the softening length ε. 

Open with DEXTER 
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).
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 > 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 cellopening 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 timestepping 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 timestepping 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 twobody 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.
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 overestimation 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 twobody 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_{0} 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_{0}.
6.1. Scaling with particles number
When the number of particles is small enough, i.e. N < N^{∗} ≈ 10^{4} 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^{4} 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^{2}: 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 τ_{2} is almost independent on N, at least for the limited range of N we considered, i.e. N = 10^{4} ÷ 10^{5}.
In Fig. 16 it is shown the density profile in physical coordinates for simulations with N = 10^{3},10^{5} 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^{2}. 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^{5} the simulation is still in the QSS regime, for the N = 10^{3} case the density profile is already well fitted by the NFW shape (see Fig. 16).
Fig. 16 Density profile in physical coordinates for simulations with 10^{3} and 10^{6} 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 = 10^{3} and N = 10^{5} 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 R_{0} we can conclude that a_{LI} = a_{LI}(R_{0}). In particular, the smaller is R_{0} and the smaller is a_{LI}.
Fig. 17 Behaviour of A(a) (Eq. (26)) for simulations with different initial size R_{0} as a function of time (normalised to τ_{c}(R_{0})). 

Open with DEXTER 
This fact can be simply interpreted by considering that when the system size decreases, the amplitude of the overdensity δ = δ(R_{0}) increases and so the value of the characteristic time scale τ_{c} = τ_{c}(R_{0}) (see Table 2). Thus, to compare the behaviours for different R_{0} 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_{0}, of the corresponding value of τ_{c} = τ_{c}(R_{0}). The residual difference in the observed behaviours is due to the effect of the different ratio R_{0}/L: when R_{0} > 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_{0} (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_{0} and in a coherent way with the behaviour of A(a).
Fig. 18 Density profile in physical coordinates at a = 2.5 for different values of R_{0} (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 R_{0} converges precisely to the same behaviour characterised by a x^{3} decay. As discussed in Sect. 5, at fixed timestepping 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_{0}, for long enough times, all curves collapse on each other.
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 R_{0} (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 b_{0} → − 1 the halo is already almost virialized while when b_{0} → 0 the halo is far from such a configuration. However, for b_{0} ∈ [ − 1,0] the structure undergoes to a collapse, which can be less or more violent, that drives the system, in about the same timescale τ_{c} (see Eq. (2)), into a virialized QSS that can be well described as a collisionless equilibrium configuration (Sylos Labini 2012).
Figure 20 shows the behaviour of A(a) for simulations with different initial virial ratio b_{0} ∈ [ − 1,0], i.e. with different initial velocity dispersion. We find that a_{LI} = a_{LI}(b_{0}): in particular, the smaller is b_{0} 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_{0} at the higher is the density of the core reached after the collapse.
Fig. 20 Behaviour of A(a) (Eq. (26)) for for simulations with different initial size b_{0}. 

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 b_{0} > − 1/2. In this situation the system evolution is very similar to the b_{0} = 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} = 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_{0} < − 1/2 the virialized structure after the collapse does not present a density profile with a 1/r^{4} tail. However, for long enough times the density profile approaches the 1/r^{3} behaviour for any value of b_{0}. This shows that the 1/r^{3} 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_{0} and the slower is the evolution. In particular, the smaller b_{0} the longer it takes to form the 1/r^{3} tail, a behaviour that is approached at long enough times for any value of b_{0}.
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 a_{LI} = a_{LI}(R_{0},N,b_{0}), 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^{4}, 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^{4},10^{5}], 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 R_{0} we find that a_{LI} decreases as well: this is coherent with the fact that the smaller is R_{0} 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_{0}) 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_{0} 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} → 0 the collapse is much more violent than for b_{0} → − 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 overdensity in a cosmological Nbody 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 overdensity 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 Nbody simulation, to an isolated structure in the stable clustering regime.
Joyce & Sylos Labini (2013) considered the evolution of an isolated overdensity 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 overdensity with initial virial ratio in the range b_{0} − ∈ [ − 1,0]. When b_{0} → − 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} → 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} = 0, and spherical overdensity 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 powerlaw density profiles.
We found that, for b_{0} → 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 τ_{2} 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} < τ_{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 timestepping η. 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 interparticle 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 timestepping 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): twobody 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 timestepping 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 Nbody 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 nonlinear 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_{0} and the initial velocity dispersion (and thus the initial virial ratio b_{0}). As result we found that When decreasing R_{0} we find that a_{LI} decreases: the initial interparticle distance gets smaller, as ℓ ∝ R_{0}, and thus hard scatterings are more probable. The dependence of a_{LI} on the particle number clearly shows the effect of twobody collisions. For small particle number, i.e. N < 10^{4}, 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^{4},10^{5}] 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 τ_{2} in Eq. (1) becomes almost independent on N.
Finally by increasing the initial velocity dispersion, i.e. decreasing b_{0}, 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_{0} = −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} = 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 twobody scatterings. However, we have found that for b_{0} ∈ [ − 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^{3} behaviour for any value of b_{0} 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 twobody 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 overdensity) and numerical (softening length, timestepping) 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^{3},10^{5}], 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 selfgravitating 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 quasistationary configuration, given the same numerical parameter of the numerical integration (softening length, timestepping, 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 Nbody 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
 Aarseth, S., Lin, D., & Papaloizou, J. 1988, ApJ, 324, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Aarseth, S. J. 2003, Gravitational NBody Simulations: Tools and Algorithms (Cambridge University Press) [Google Scholar]
 Bertin, G. 2000, Dynamics of Galaxies (Cambridge University Press) [Google Scholar]
 Boily, C., & Athanassoula, E. 2006, MNRAS, 369, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1 [NASA ADS] [CrossRef] [Google Scholar]
 David, M., & Theuns, T. 1989, MNRAS, 240, 957 [NASA ADS] [Google Scholar]
 Diemand, J., Moore, B., Stadel, J., & Kazantzidis, S. 2004, MNRAS, 977 [Google Scholar]
 Gabrielli, A., Baertschiger, T., Joyce, M., et al. 2006, Phys. Rev., 74, 021110 [NASA ADS] [Google Scholar]
 Gabrielli, A., Joyce, M., & Marcos, B. 2010, Phys. Rev. Lett., 105, 210602 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Heggie, D., & Hut, P. 2003, The Gravitational MillionBody Problem: A Multidisciplinary Approach to Star Cluster Dynamics, Cambridge [Google Scholar]
 Hénon, M. 1964, Ann. Astrophys., 27, 1 [Google Scholar]
 Heitmann, K., Ricker, P. M., Warren, M. S., & Habib, S. 2005, ApJS, 160, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Knebe, A., Kravtsov, A. V., Gottlober, S., & Klypin, A. A. 2000, MNRAS, 317, 630 [NASA ADS] [CrossRef] [Google Scholar]
 Knebe, A., Green, A., & Binney, J. 2001, MNRAS, 325, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Iannuzzi, F., & Dolag, K. 2011, MNRAS, 417, 2846 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., & Sylos Labini, F. 2013, MNRAS, 429, 1088 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., Marcos, B., & Baertschiger, T. 2009a, MNRAS, 394, 751 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., Marcos, B., & Sylos Labini, F. 2009b, MNRAS, 397, 775 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., Marcos, B., & Sylos Labini F. 2009c, J. Stat. Mechanics, P04019 [Google Scholar]
 Levin, Y., Pakter, R., & Rizzato, F. 2008, Phys. Rev., E78, 021130 [NASA ADS] [CrossRef] [Google Scholar]
 Miniati, F., & Colella, Ph., 2007, J. Comput. Phys., 227, 400 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P. J. E. 1980, The LargeScale Structure of the Universe (Princeton University Press) [Google Scholar]
 Power, C., Navarro, J. F., Jenkins, A. R., et al. 2003, MNRAS, 338, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Romeo, A. B., Agertz, O., Moore, B., & Stadel, J. 2008, ApJ, 686, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Splinter, R. J., Melott, A. L., Shandarin, S. F., & Suto, Y. 1998, ApJ, 497, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., Yoshida, N., & White, S. D. M. 2001, New Astron., 6, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., White, S. D. M., & Jenkins, A., et al. 2005, Nature, 435, 629 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Stiavelli, M., & Bertin, G. 1987, MNRAS, 229, 61 [NASA ADS] [Google Scholar]
 Sylos Labini, F. 2012, MNRAS 423, 1610 [Google Scholar]
 Sylos Labini, F. 2013, MNRAS, 429, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2001, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Theis, C., & Spurzem, R. 1999, A&A, 341, 361 [NASA ADS] [Google Scholar]
 Theuns, T., & David, M. 1990, Astrophys. Sp. Sci., 170, 276 [NASA ADS] [CrossRef] [Google Scholar]
 van Albada, T. 1982, MNRAS, 201, 939 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Upper panel: energy in physical coordinates (K_{p} Eq. (15) and W_{p} 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 
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 
Fig. 3 Upper panel: absolute value of the potential energy in physical coordinates W_{p}(a) (Eq. (16)) as a function of the scale factor. The approximate value of a_{LI}, separating the regime where W_{p} ≈ const. from the one where  W_{p}  ~ 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 
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 
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 
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 
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 
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 
Fig. 9 Behaviour of A(a) (Eq. (26)) for simulations with different ε (see Table 1 for details). 

Open with DEXTER  
In the text 
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 
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 
Fig. 12 Behaviour of A(a) (Eq. (26)) for simulations with different η and ε (see labels). 

Open with DEXTER  
In the text 
Fig. 13 Behaviour of a_{LI} as a function of η for different values of the softening length ε. 

Open with DEXTER  
In the text 
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 
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 
Fig. 16 Density profile in physical coordinates for simulations with 10^{3} and 10^{6} 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 = 10^{3} and N = 10^{5} are reported as reference. 

Open with DEXTER  
In the text 
Fig. 17 Behaviour of A(a) (Eq. (26)) for simulations with different initial size R_{0} as a function of time (normalised to τ_{c}(R_{0})). 

Open with DEXTER  
In the text 
Fig. 18 Density profile in physical coordinates at a = 2.5 for different values of R_{0} (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 
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 R_{0} (see labels). A reference line with slope − 4 is plotted as reference. 

Open with DEXTER  
In the text 
Fig. 20 Behaviour of A(a) (Eq. (26)) for for simulations with different initial size b_{0}. 

Open with DEXTER  
In the text 
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 