Cusps in the center of galaxies: a real conflict with observations or a numerical artefact of cosmological simulations?
A. N. Baushev, L. del Valle, L. E. Campusano, A. Escala, R. R. Muñoz, G. A. Palma
aa r X i v : . [ a s t r o - ph . GA ] M a y Cusps in the center of galaxies: a real conflict with observations or a numericalartefact of cosmological simulations?
A. N. Baushev
Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 141980 Dubna, Moscow region, Russia andDepartamento de Astronom´ıa, Universidad de Chile, Casilla 36-D, Correo Central, Santiago, Chile
L. del Valle, L.E. Campusano, A. Escala, and R.R. Mu˜noz
Departamento de Astronom´ıa, Universidad de Chile, Casilla 36-D, Correo Central, Santiago, Chile
G.A. Palma
Departamento de F´ısica, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile (Dated: May 12, 2017)Galaxy observations and N-body cosmological simulations produce conflicting dark matter halodensity profiles for galaxy central regions. While simulations suggest a cuspy and universal densityprofile (UDP) of this region, the majority of observations favor variable profiles with a core in thecenter. In this paper, we investigate the convergency of standard N-body simulations, especiallyin the cusp region, following the approach proposed by [1]. We simulate the well known Hernquistmodel using the SPH code Gadget-3 and consider the full array of dynamical parameters of theparticles. We find that, although the cuspy profile is stable, all integrals of motion characterizingindividual particles suffer strong unphysical variations along the whole halo, revealing an effectiveinteraction between the test bodies. This result casts doubts on the reliability of the velocitydistribution function obtained in the simulations. Moreover, we find unphysical Fokker-Planckstreams of particles in the cusp region. The same streams should appear in cosmological N-bodysimulations, being strong enough to change the shape of the cusp or even to create it. Our analysis,based on the Hernquist model and the standard SPH code, strongly suggests that the UDPs generallyfound by the cosmological N-body simulations may be a consequence of numerical effects. A muchbetter understanding of the N-body simulation convergency is necessary before a ’core-cusp problem’can properly be used to question the validity of the CDM model.
I. INTRODUCTION
Results of N-body simulations come into increasingconflict with observations of the dark matter (DM) dis-tribution in the central regions of dwarf galaxies. Astro-nomical observations favor relatively soft cored densityprofiles [2–8]. On the contrary, N-body simulations ofcold dark matter tell us that dark matter halos have auniversal shape, independent of the halo mass and initialdensity fluctuation spectrum, and that the central uni-versal density profile (hereafter UDP) is cuspy. The firstworks on the subject proposed the Navarro-Frenk-Whiteprofile (hereafter NFW) that behaves as ρ ∝ r − at thecenter. Later simulations [9, 10] favor the Einasto pro-file with a finite central density. However, the obtainedEinasto index is so high (typically n ∼ −
6) that theprofile is still cuspy and close to the NFW one.For a time, there was a hope that the ’core-cusp prob-lem’ would disappear once the baryon contribution istaken into account. However, recent simulations includ-ing baryon matter have rather amplified the problem [11]:apart from the profile disagreement, a more fundamentaldifficulty was found. Of course, the presence of baryonsin simulations changes the density profile, but it remainsalmost universal for all the halos, while the profiles ofreal galaxies are extremely varied. The conflict betweensimulations and observations might suggest that the coldDM paradigm is wrong. However, before reaching this conclusion, the accuracy and convergence of the simula-tions should be scrutinized. For instance, the overesti-mation of the energy exchange between the test bodiesthat may occur in the N-body simulations leads to thecusp formation [12]. If the energy evolution during thehalo formation is limited, then the density profile of theformed halo resembles more closely the observed one [13].As an example, the overestimation of the particle en-ergy exchange may be due to the unphysical pair colli-sions of the test bodies. Its importance may be charac-terized by the relaxation time [14, eqn. 1.32] τ r = N ( r )8 ln Λ · τ d (1)where N ( r ) is the number of test bodies inside a sphereof radius r , ln Λ is the Coulomb logarithm, τ d =(6 πG ¯ ρ ( r )) − / is the characteristic dynamical time of thesystem at radius r , ¯ ρ ( r ) is the average density inside r .Equation 1 has two important consequences. First, τ r de-pends on the smoothing radius of the N-body simulationsonly through Λ, i.e., only logarithmically [1]. Therefore,the influence of the unphysical collisional relaxation can-not be decreased much by the smoothing of the test bodypotentials. Second, since the number of dark matter par-ticles is huge ( ∼ , if dark matter consists of elemen-tary particles), the collisional relaxation plays no role innature, being a purely numerical effect.The algorithm stability is the critical point of N-bodysimulations: the Miller’s instability makes the Liapunovtime comparable with the dynamical time of the sys-tem [15]. Even if we take into account the specificityof N-body algorithms, like the potential smoothing, theinstability development time is much shorter than τ r and remains comparable with the dynamical time at thegiven radius τ d ( r ) [16, 17]. However, different N-bodycodes, with various versions of the Poisson solvers, inte-gration algorithms etc., lead to final halos with the above-mentioned UDP, which is almost the same and close toNFW. Therefore, it is widely believed that the universalprofile is physically meaningful and that it describes realhalos, even though the orbits of individual test bodieshave no physical significance [14, section 4.7.1(b)]. Theaim of this paper is to question this opinion.Indeed, the convergency criteria of N-body simulationsused at present are exclusively based on the density pro-file stability. [18] found that the cusp of the UDP re-mains stable at least until t = 1 . τ r and then a coreforms. On this basis [18] supposed that the core forma-tion is the first sign of the collision influence and offeredthe most extensively used criterion for simulation conver-gency t < . τ r . The acceptance that the collisions haveno effect even if the simulation time exceeds τ r seemssurprising. However, later convergency tests (also basedonly on the stability of the density profile) suggested evensofter criteria [19, 20]. In this paper we perform a moresophisticated convergency test, going beyond the densityprofile analysis and considering the full array of the dy-namical parameters of the particles. II. CALCULATIONSA. The main idea
In order to test the N-body convergency, we follow themethod offered in [1]. We simulate the well-known Hern-quist model with the density profile ρ ( r ) = M a/ [2 πr ( r + a ) ] (where a is the scale radius and M is the total halomass), and with the isotropic velocity distribution at eachpoint [21]. The model is spherically symmetric and fullystable, i.e., the density and velocity profiles should notchange with time. We chose the Hernquist model becauseit is close to the NFW and behaves exactly as the NFW( ρ ∝ r − ) in the central region, but it has a known an-alytical solution for the stationary velocity distribution,contrary to the NFW one. The region of the cusp ( r < a )is of main concern to us.Since the gravitational potential φ ( r ) is constant, thespecific energy ǫ = φ ( r ) + v / ~K of each particle should be conserved. In-stead of ǫ , it will be more convenient to use the apocenterdistance of the particle r (i.e. the maximum distance onwhich the particle can move off the center, which can befound from the implicit equation ǫ = φ ( r ) + K / r ).Being an implicit function of the integrals of motion ǫ and K , r is an integral of motion as well. Thus, any timevariation of ǫ , ~K , or r is necessarily a numerical effect, % / r r/a FIG. 1: The ratio between the time t when the mass insidesome radius r drops on more than 20%, to the relaxation time τ r ( r ). and we may judge the simulation convergency followingthe behavior of these quantities.We need to clarify two important points of our work.Some properties of the perfectly symmetrical model weconsider (like the exact conservation of the angular mo-mentum for every particle) are unstable and not realisticfor real astrophysical DM halos that are always triaxialas a result of tidal perturbations etc. The application ofperfectly spherical models to real systems may give riseto false conclusions [22]. However, the use of the spher-ical model for our purposes is well founded. We are notconsidering the task of comparison of simulation resultswith observations. Our aim is just to check if the ’N-body matter’ behaves as a collisionless matter, which isthe principle question of the dark matter modelling.Second, there is a frequent belief that it is much eas-ier to converge on the spherically averaged density dis-tribution than on the full properties of the phase spacedistribution function. Indeed, we need not correctly re-produce each individual particle trajectory. Moreover, itis not even necessarily desirable since real dark matterhalos are not spherically symmetric and therefore hostchaotic orbits. However, it would be completely wrongto disregard the phase evolution of the system or considerits evolution as a ’second order effect’ with respect to thedensity profile shape. As we will show in the Results:the simulation convergency section, correct simulationsof the energy and angular momentum of each particle(contrary to individual particle trajectories) are of criti-cal importance for correct simulations of the density pro-file.
B. The simulations
We simulate a single separate Hernquist halo. Theaim of this work is to perform a sophisticated test of thestandard convergency criteria, therefore we do not try tomodel any real astronomical object. Since the standardN-body units [23] are used, the results are independenton the choice of a and halo mass. However, we choosesome values of the parameters, for illustrative purposes.Let us set a = 100 pc, which roughly corresponds tothe well-known dwarf spheroidal satellite of the MilkyWay, Segue 1. This is one of the most popular objectsfor the indirect dark matter search, since it is close tothe Solar System; its present-day mass can be estimatedas 3 · M ⊙ [24]. Segue 1 experienced strong tidal dis-ruption, and we do not know its initial mass. We con-sider two limiting cases. In the body of the paper weaccept the halo mass M = 10 M ⊙ , which is compara-ble with the present-day mass of a larger dwarf satellite,Fornax [25] and almost certainly exceeds the initial massof Segue 1. Thus we consider the case of a compact andvery dense dwarf spheroidal galaxy. However, since allthe simulations are performed in the dimensionless N-body units, a reader may easily extend the results forany value of M . If a is fixed, the only value that is sen-sitive to the choice of M is time: all the time intervalsscale as ∆ t ∝ M − . (while the ratios of time intervalsremain the same). As an illustration, we also consid-ered the case of M = 10 M ⊙ , which is certainly lower,then the present-day mass of Segue 1. The only differ-ence is that all the time intervals get ten times larger,and we everywhere specify the values corresponding to M = 10 M ⊙ in the footnotes. Anticipating events, wesay that the results shown in all the figures in this paperare not sensitive at all to the choice of M .We use N = 10 test bodies [35]. They are placedrandomly, in accordance with the analytically obtainedspace and velocity distributions [21]. The relaxationtime at r = a is τ r ( a ) ≃ . · s ≃ . · years.Therefore, we make 200 snapshots with the time interval∆ t = 10 s ≃
30 mln. years, covering the time from 0 to t max = 2 · s ≃ . · years (for the case of the halomass M = 10 M ⊙ , τ r ( a ) ≃ . · years, ∆ t = 10 s ≃
300 mln. years, t max = 2 · s ≃ . · years). Werecord the positions and the velocities of each particle oneach snapshot.We evolve the system using one of the most exten-sively employed in cosmological simulations SPH codes,Gadget-3, an update version of Gadget-2 [26, 27]. Thegravitational interactions in Gadget-3 are computed us-ing a hierarchical tree [28, 29]. In this algorithm the spaceis divided in different cells and the gravitational force act-ing on a particle is computed using a direct summationfor particles that are in the same cell and by means ofmultipole (up to the quadrupole) expansion for the par-ticles that are in a different cell. The minimum distancebetween particles to be part of a different cell is con-trolled by a tree opening criterion. Gadget-3 uses theBarnes-Hut tree opening criterion for the first force com-putation. This criterion is controlled by an opening angle µ , which determines the maximum ratio between the dis-tance to the center of mass of the cell ( d ) and the size r /a FIG. 2: The averaged relative variations h b ∆ K/K circ i (squares) and h b ∆ r /r i (crosses) of the integrals of motionin a single time step ∆ t . See the definition of the averaging b ∆ in the section Results: the integrals of motion . of the cell ( l ). If the cell is too close to the particle, d/l will be greater than µ , and new cells have to be openedto maintain the accuracy on the force computation. Inthe further evolution of the system a dynamical updat-ing criterion (controlled by the fractional error f acc ) isused. We use the standard set of parameters µ = 0 . f acc = 0 . ∼ . . a = 2 pc, in accordance with [19, 32]. III. RESULTS: THE INTEGRALS OF MOTION
Initially we convert each of the 201 snapshots into thecenter-of-mass frame of references at the moment whena snapshot is made.First of all, we try to reproduce the results of [18]. Thedensity profile indeed remains quite stable, and then acore in the center appears. Exactly following [18], weconsider the moment t when the mass inside someradius drops on 20% comparing to the initial value asthe moment of the core formation. The ratio of t tothe relaxation time τ r at the same radius r is representedin Fig. 1. We see that our data by and large confirm theresults of [18], the core really appears at t ≃ τ r .Before proceeding any further, two important com-ments relating to all the subsequent text should be made.First, our convergency tests are mainly oriented on theradius interval [0 . a ; 1 . a ] where they are the most pre-cise. This choice of the working interval might appearstrange at first sight: typically the convergency prob-lems occur much closer to the halo center. However, ifwe had chosen a realistic area ( r ≤ . a ), then it wouldhave contained only ∼
100 test particles, and the statisticwould have been poor. On the other hand, the density r /a FIG. 3: The ratios K circ τ r D b ∆ K ∆ t E − (squares) and τ r D b ∆ r r ∆ t E − (crosses) of the time, in which the particles completely ’forget’the initial values of the integrals of motion, to the relaxationtime τ r ( r ). profile between 0 . a and 0 . a remains much the sameas in the center, since a power-law profile ρ ∝ r − is self-similar. The lower border of the region under consider-ation r = 0 . a is defined by our choice of the timestep∆ t = 10 s (for the case of the halo mass M = 10 M ⊙ ,∆ t = 10 s ≃
300 mln. years). At r = 0 . a , τ r ≃ ∆ t , andthe Hernquist profile is certainly corrupted by the colli-sions even on the first timestep. However, as we will seefrom the discussion of fig. 4, the core formation becomesvisible in phase portrait at much larger distances thanin the density profile itself. Therefore, only the resultsrelated to r ≥ . a can be totally trusted.Second, we want to consider variations of the integralsof motion as a function of radius. However, each particlecontributes to the density profile on an interval betweenits pericenter radius r min and apocenter radius r . Here-after we will consider r as the characteristic radius cor-responding to the particle. Indeed, if the particle orbitis elongated, the particle spends almost all the time nearthe apocenter, in accordance with the Kepler’s secondlaw. On the contrary, if the orbit is circular, the parti-cle moves along almost uniformly, but its radius alwaysremains close to r .In order to study the behavior of the integrals of mo-tion (theoretically they should conserve), we order all10 particles according to their r in the initial snap-shot, and then divide the particles into 200 groups of5000 particles each. All the particles in the same grouphave similar r , and the group may be characterizedby the average initial r of its members. We calcu-late ∆ r /r = ( r ( i + 1) − r ( i )) /r and ∆ K/K circ =( K ( i +1) − K ( i )) /K circ for each particle on each timestep.Here i is the number of the snapshot, K circ is the angu-lar momentum corresponding to the circular orbit at r ;apparently, this is the maximum value of K any particlewith the apocenter distance r may possess. Then we r/a N / N (r) FIG. 4: The upward ∆ N + ( r ) / ∆ t (squares) and downward∆ N − ( r ) / ∆ t (crosses) Fokker-Planck streams of particles di-vided by the total number of the particles N ( r ) inside r . find the root-mean-squares of ∆ r /r and ∆ K/K circ av-eraged over each group and for each snapshot. Our analy-sis shows that the root-mean-squares do not significantlydepend on time until the moment when the core formsat the radius corresponding to r of the group. There-fore, we then average the root-mean-squares of ∆ r /r and ∆ K/K circ over all the timesteps where the core hadnot formed yet. We denote the values averaged in sucha complex manner by h b ∆ r /r i and h b ∆ K/K circ i .The dependance of h b ∆ K/K circ i (squares) and h b ∆ r /r i (crosses) from the dimensionless radius r /a isrepresented in Fig. 2. We see that even in a single timestep ∆ t = 10 s ≃
30 mln. years (for the case of thehalo mass M = 10 M ⊙ , ∆ t = 10 s ≃
300 mln. years)the integrals (that should be constant) vary significantly.Fig. 3 represents the values K circ τ r D b ∆ K ∆ t E − (squares) and τ r D b ∆ r r ∆ t E − (crosses) that are the ratios of the time in-tervals in which an average particle totally ’forgets’ itsinitial values of K and r to the relaxation time τ r ( r ).Everywhere in the region of reliability ( r ≥ . a ) theratios are much less than 1.It means that the particles totally ’forget’ their inte-grals of motion in a time much shorter than τ ( r ). Ingeneral one could not expect a reliable simulation of thevelocity distribution at t ∼ τ ( r ) under such conditions. IV. RESULTS: THE SIMULATIONCONVERGENCY
Thus, the integrals of motion of the particles are notconserved at all, while the density profile remains sta-tionary in quite good agreement with the theory in oursimulations (we should mention, however, that the ab-sence of the collision influence up to almost 2 relaxation r/a
FIG. 5: The fractions of N ( r ) particles inside radius r thatis carried away and in by the upward and downward Fokker-Planck streams in the time 1 . τ r offered by [18] (1 . τ r ∆ N + ( r ) N ( r )∆ t (squares) and 1 . τ r ∆ N − ( r ) N ( r )∆ t (crosses), respectively). times looks surprising [1]. The same stability and repro-ducibility of the cusps in cosmological modelling leads tothe wide acceptance of the idea that, though no signif-icance can be attached to the trajectories of individualparticles in the N-body simulations, the cuspy densityprofile is meaningful and should correctly describe theprofiles of real halos. Let us use our results to illustratethe vulnerability of the profile stability as the only con-vergency criterion of the N-body simulations.Indeed, if a Hernquist halo consists of real DM (wesuppose that it is cold and noninteracting), the values of ǫ , ~K , and r of each particle must conserve, the particledistribution function f should depend only on ǫ and K [14] and obey the collisionless kinetic equation df /dt = 0.It means that there are no particle fluxes in the phasespace ( ǫ, K ).However, figures 2 and 3 doubtlessly reveal an inten-sive energy and angular momentum exchange betweenthe particles, i.e., the test bodies interact. Then the sys-tem may be described by the Fokker-Planck (hereafterFP) equation [33] dfdt = ∂∂q α (cid:26) ˜ A α f + ∂∂q β [ B αβ f ] (cid:27) (2)where q α is an arbitrary set of generalized coordinates,˜ A α = δq α δt B αβ = δq α δq β δt (3)We may choose q = ǫ , q = K , and figure 2 shows thatat least coefficients B and B in the equation (2) differessentially from zero. Thus we model real DM halos thatare believed to be collisionless, by a system of test bodiesgoverned by the kinetic equation with a significant colli-sional term, i.e., by an essentially collisional equation. An important point must be underscored: the densityprofile in our simulation indeed holds its shape (closeto the NFW one in the center), in gratifying agreementwith the theoretical predictions. The variations of the in-tegrals of motion, that we found, mainly touches on thevelocity distribution of the particles. Together with theUDP stability in cosmological simulations, it can producea dangerous illusion that N-body simulations might ade-quately model the density profiles of dark matter struc-tures, despite the fact that the velocity distribution wasdistorted. We should emphasize that it cannot be true.Indeed, let us consider a stationary spherically symmetrichalo for the sake of simplicity. The particle distributionin the phase space f d xd v is a function of only the parti-cle energy ǫ and three components if its angular momen-tum ~K [14]. If the velocity distribution of the particlesis anisotropic in each point (which is the case under con-sideration in this paper), f depends only on the particleenergy f ( ep ) = f ( φ ( r ) + v /
2) = f ( φ ( r ) + K / (2 r )).The particle speed distribution at some radius r is there-with equal to 4 πv f ( φ ( r ) + v / dv , and the densityis R πv f ( φ ( r ) + v / dv . These relationships clearlydemonstrate the impossibility of a reliable determinationof the density profile without a reliable determination ofthe velocity distribution. The distributions over ~v and r are not just bound, there are actually a sort of projectionsof the same distribution f on the velocity or space coor-dinates. Apparently, this conclusion is very general anddoes not depend on the assumption about the sphericalsymmetry that we made.We are able to compare the results with the theoreticalprediction and check their agreement in the model casethat we consider. However, it is impossible in the caseof real cosmological simulations. Therefore, any numer-ical effect influencing on the velocity distribution f ( ~v )or on the integrals of motion of the test particles putsthe density profiles obtained in the simulations in doubt.Moreover, the fact that the cuspy profile ρ ∝ r − turnsout to be very stable in our simulation, despite of r and K variations, is probably not a coincidence, but a directresult of the numerical effects we discuss.Indeed, the fact that we model collisionless systemswith the test bodies governed by an essentially collisionalequation is surprising per se , but the main consequenceis that the profile stability does not guarantee the sim-ulation correctness. The FP streams in the phase spacecreated by the particle interaction may form stable den-sity profiles (corresponding to the stationary solutions ofthe Fokker-Planck equation), but these profiles and theirpersistence are at odds with the behavior of real colli-sionless systems. As the first and crude illustration, thecollisions lead to the contraction of the central region ofany realistic profile and finally to the core collapse. Thedensity profile outside the core approaches a power law ρ ∝ r − . and then remains quite stable for a long time[14]. Of course, this distribution is already formed bythe unphysical test body collisions, and the immutabil-ity of the ρ ∝ r − . profile says nothing either aboutthe simulation convergency or about the behavior of realcollisionless systems.The core collapse appears at t ≫ τ r and has nothingto do with the Hernquist or the UDP profiles. However,the Fokker-Planck equation has an another stationarysolution close to the NFW one [1, 34].A question appears: if we obtain a stable cuspy den-sity profile, how can we differentiate cusps correlatingwith the properties of real collisionless systems from thesolutions created by the numerical effects? In a collision-less system, the values of r of the particles in the cuspshould remain constant. If collisions are significant, thevalues of r should experience a random walking, and theparticles move up and down in the cusp forming a down-ward stream (of the particles with decreasing r ) and anupward stream (of the particles with increasing r ). Forthe cusp to be stable, the streams should compensateeach other, which corresponds to a stationary solutionof the Fokker-Planck equation. Thus, if the cusp is cre-ated by the FP diffusion, we should see two significantstreams of particles with decreasing and increasing r ,and the streams should compensate each other in orderto provide the cusp stability.We chose two adjacent (i.e., divided by a single ∆ t )snapshots at the beginning of the simulations, in order tominimize the core formation effects. For an array of radii r , we calculated the number ∆ N + ( r ) of particles thathad r < r at the first snapshot and r > r at the secondone, and the number ∆ N − ( r ) of particles that had r > r at the first snapshot and r < r at the second one. Ofcourse, ∆ N + ( r ) = ∆ N − ( r ) = 0 in the collisionless case,since r is an integral of motion.Fig. 4 represents ∆ N + ( r ) (squares) and ∆ N − ( r )(crosses) divided by the total number of the particles N ( r ) inside r . As we can see, the FP streams ex-ist, though they compensate well each other outside of r = 0 . a . As we approach the center, the upward streambecomes increasingly stronger than the downward one.This is the first sign of the core formation, that is stillinvisible in the density profile at this radius, being al-ready quite clear at the phase picture of the system.A question appears: are the discovered fluxes ∆ N + ( r )and ∆ N − ( r ) real and important? May they be just asmall noise, produced by particles near the boundary,crossing and recrossing it and thus giving the impressionof flows that do not exist? One can readily see thatthis is not the case. First of all, ∆ N + ( r ) and ∆ N − ( r )apparently give only the lower bounds on the upwardand downward FP streams: the value of r of a particlecould have crossed r an odd number of times (and thenit is counted only once) or an even number of times (andthen it is not counted at all). Since we count each particleno more than once on a timestep, we totaly avoid therecrossing effect.Second, the flows are just too strong to be just anoise. For instance, Fig. 4 shows that, though ∆ N + ( r )and ∆ N − ( r ) are just the lower bounds on the streams,∆ N + ( r ) ≃ ∆ N − ( r ) ≃ · at r = a , i.e. ∼
2% of the total halo mass crosses this radius because of this unphys-ical effect on each timestep. This is approximately thetotal number of particles in the layer of thickness ∼ a/ r = a . The value a/
12 by far exceedsthe smoothing radius or any reasonable numerical noisethat may occur in the computing scheme.The surprisingly high intensity of the Fokker-Planckdiffusion is the main result of this work. Approximately8% of particles are renewed even inside r = a . It meansthat in only 10∆ t ≃
300 mln. years (for the case of thehalo mass M = 10 M ⊙ , 10∆ t ≃ · years), i.e., in 5%of the simulation time, all the particles inside the sphere r = a (which contains a quarter of the total mass ofthe system) can be substituted by a purely numerical ef-fect. The fractions of particles inside radius r that can becarried away or in by the upward and downward Fokker-Planck streams in the Power’s time 1 . τ r are 1 . τ r ∆ N + ( r ) N ( r )∆ t and 1 . τ r ∆ N − ( r ) N ( r )∆ t . Fig. 5 shows that they always signif-icantly exceed 1. The cusp in our simulations was cre-ated in the initial conditions, but its shape is similar tothe UDP. We can see that the unphysical FP streams arestrong enough to arbitrarily change the shape of the cusp(and therefore the shape is defined by the FP diffusionrather than by the properties of the collisionless system)and even to create it.Another argument in support of the numerical natureof the cusps in cosmological simulations is the profile uni-versality. The similarity contradicts the observational re-sults [11], but is quite natural if the cusps are formed bythe FP diffusion. An UDP-like stationary solution is in-nate for the FP equation: the suppositions in [34] and[1] are quite different, but the results are similar. Theproperties of the solution of the FP equation are totallydefined by only a few coefficients ˜ A α and B αβ that canbe similar for different N-body codes using similar al-gorithms, and are almost certainly the same within thesame simulation. As a result, the resulting halos are alsoself-similar, while nature is much more variable.The second important conclusion of this section is thatthe profile stability cannot be used as the simulation con-vergency criterion: the first unquestionable signs of theinfluence of the test particle interaction appear in thephase portrait much earlier than the density profile evo-lution and the beginning of the core formation.The third conclusion is that, since the variations of K and r are very significant in Fig. 2 even at r > a , wherethe role of collisions or potential softening is minor, theintegral variations there are most likely due to the po-tential calculating algorithm. But whatever the reasonof the variations may be, the ill effect on the simulationsis the same from the point of view of the kinetic equa-tion: variations if the integrals of motion reveal the colli-sional influence and suggest that the system behavior isno longer described by the correct collisionless equation.A convergence study with varying the opening angle µ ofthe Barnes-Hut tree opening criterion, softening scale, aswell as other parameters of the gravitational force com-putation, is essential to understanding the origin of thenon-conservation of integrals of motion and find the opti-mal parameter set to decrease these undesirable numeri-cal effects.A much better understanding of the N-body simula-tion convergency is necessary to cast doubts on the CDMmodel on the basis of ’cusp vs. core’ contradiction. V. ACKNOWLEDGEMENTS
The work is supported by the CONICYT Anilloproject ACT-1122 and the Center of Excellence in As- trophysics and Associated Technologies CATA (PFB06).We used the HPC clusters Docorozco (FONDE-CYT1130458) and Geryon(2) (PFB06, QUIMAL130008and Fondequip AIC-57). [1] A. N. Baushev, Astroparticle Physics , 47 (2015),1312.0314.[2] L. Chemin, W. J. G. de Blok, and G. A. Mamon, AJ , 109 (2011), 1109.4247.[3] W. J. G. de Blok, S. S. McGaugh, and V. C. Rubin, AJ , 2396 (2001).[4] W. J. G. de Blok and A. Bosma, A&A , 816 (2002),arXiv:astro-ph/0201276.[5] S.-H. Oh, W. J. G. de Blok, E. Brinks, F. Walter, andR. C. Kennicutt, Jr., AJ , 193 (2011), 1011.0899.[6] F. Governato, A. Zolotov, A. Pontzen, C. Christensen,S. H. Oh, A. M. Brooks, T. Quinn, S. Shen, and J. Wad-sley, MNRAS , 1231 (2012), 1202.0554.[7] E. J. Tollerud, R. L. Beaton, M. C. Geha, J. S. Bullock,P. Guhathakurta, J. S. Kalirai, S. R. Majewski, E. N.Kirby, K. M. Gilbert, B. Yniguez, et al., Astrophys. J. , 45 (2012), 1112.1067.[8] A. Del Popolo and F. Pace, Ap&SS , 162 (2016),1502.01947.[9] J. Stadel, D. Potter, B. Moore, J. Diemand, P. Madau,M. Zemp, M. Kuhlen, and V. Quilis, MNRAS , L21(2009), 0808.2981.[10] J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vo-gelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, andA. Helmi, MNRAS , 21 (2010), 0810.1522.[11] K. A. Oman, J. F. Navarro, A. Fattahi, C. S. Frenk,T. Sawala, S. D. M. White, R. Bower, R. A. Crain,M. Furlong, M. Schaller, et al., MNRAS , 3650(2015), 1504.01437.[12] A. N. Baushev, Astrophys. J. , 65 (2014), 1205.4302.[13] A. N. Baushev, A&A , A114 (2014), 1309.5162.[14] J. Binney and S. Tremaine, Galactic Dynamics: SecondEdition (Princeton University Press, 2008).[15] R. H. Miller, Astrophys. J. , 250 (1964).[16] M. Valluri and D. Merritt, Advanced Series in As-trophysics and Cosmology , 229 (2000), astro-ph/9909403.[17] P. Hut and D. C. Heggie, Journal of Statistical Physics , 1017 (2002), astro-ph/0111015.[18] C. Power, J. F. Navarro, A. Jenkins, C. S. Frenk, S. D. M. White, V. Springel, J. Stadel, and T. Quinn, MNRAS , 14 (2003), astro-ph/0201544.[19] E. Hayashi, J. F. Navarro, J. E. Taylor, J. Stadel,and T. Quinn, Astrophys. J. , 541 (2003), astro-ph/0203004.[20] A. Klypin, F. Prada, G. Yepes, S. Hess, and S. Gottlober,ArXiv e-prints (2013), 1310.3740.[21] L. Hernquist, Astrophys. J. , 359 (1990).[22] A. Pontzen, J. I. Read, R. Teyssier, F. Governato,A. Gualandris, N. Roth, and J. Devriendt, MNRAS ,1366 (2015), 1502.07356.[23] M. H. H´enon, Ap&SS , 151 (1971).[24] A. N. Baushev, S. Federici, and M. Pohl, Phys. Rev. D , 063521 (2012), 1205.3620.[25] M. G. Walker and J. Pe˜narrubia, Astrophys. J. , 20(2011), 1108.2404.[26] V. Springel, N. Yoshida, and S. D. M. White, New As-tronomy , 79 (2001), astro-ph/0003162.[27] V. Springel, MNRAS , 1105 (2005), astro-ph/0505010.[28] J. Barnes and P. Hut, Nature (London) , 446 (1986).[29] L. Hernquist and N. Katz, ApJS , 419 (1989).[30] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow,A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, andS. D. M. White, MNRAS , 1685 (2008), 0809.0898.[31] M. Boylan-Kolchin, V. Springel, S. D. M. White,A. Jenkins, and G. Lemson, MNRAS , 1150 (2009),0903.3041.[32] E. van Kampen, ArXiv Astrophysics e-prints (2000),astro-ph/0002027.[33] L. D. Landau and E. M. Lifshitz, Statistical physics. Pt.1,Pt.2 (1980).[34] N. W. Evans and J. L. Collett, ApJL , L103 (1997),astro-ph/9702085.[35] All the data, as well as results of simulations ofa Plummer sphere of mass 10 M ⊙ we used asan auxiliary test model, are publicly available atwe used asan auxiliary test model, are publicly available at