Issue 
A&A
Volume 621, January 2019



Article Number  A8  
Number of page(s)  14  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201833460  
Published online  19 December 2018 
Phasespace structure analysis of selfgravitating collisionless spherical systems
^{1}
UPMCCNRS, UMR7095, Institut d’Astrophysique de Paris, 98 bis boulevard Arago, 75014 Paris, France
email: anaelle.halle@obspm.fr; colombi@iap.fr
^{2}
Max Planck Institut für Astrophysik, KarlSchwarzschildStrasse 1, 85741 Garching bei München, Germany
^{3}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Bd de l’Observatoire, CS 34229, 06304 Nice Cedex 4, France
Received:
21
May
2018
Accepted:
3
October
2018
In the mean field limit, isolated gravitational systems often evolve towards a steady state through a violent relaxation phase. One question is to understand the nature of this relaxation phase, in particular the role of radial instabilities in the establishment/destruction of the steady profile. Here, through a detailed phasespace analysis based both on a spherical Vlasov solver, a shell code, and a Nbody code, we revisit the evolution of collisionless selfgravitating spherical systems with initial powerlaw density profiles ρ(r) ∝ r^{n}, 0 ≤ n ≤ −1.5, and Gaussian velocity dispersion. Two subclasses of models are considered, with initial virial ratios η = 0.5 (“warm”) and η = 0.1 (“cool”). Thanks to the numerical techniques used and the high resolution of the simulations, our numerical analyses are able, for the first time, to show the clear separation between two or three wellknown dynamical phases: (i) the establishment of a spherical quasisteady state through a violent relaxation phase during which the phasespace density displays a smooth spiral structure presenting a morphology consistent with predictions from selfsimilar dynamics, (ii) a quasisteadystate phase during which radial instabilities can take place at small scales and destroy the spiral structure but do not change quantitatively the properties of the phasespace distribution at the coarse grained level, and (iii) relaxation to a nonspherical state due to radial orbit instabilities for n ≤ −1 in the cool case.
Key words: gravitation / Galaxy: kinematics and dynamics / dark matter
© ESO 2018
1. Introduction
Dark matter in the universe and stars in galaxies behave like a selfgravitating collisionless fluid of which the dynamics can be described by the VlasovPoisson system:
where f(x, v, t) is the phasespace density of the fluid at position x, velocity v, and time t, ρ is the mass density and ϕ is the gravitational potential.
A major issue when considering the dynamics of gravitational systems such as darkmatter halos, elliptical galaxies, or star clusters in the noncollisional regime is to understand the main processes underlying the creation of the quasistationary states that build up after a number of dynamical times, for instance the universal profiles of darkmatter halos (Navarro et al. 1996, 1997).
One way to relate initial to quasiequilibrium state is to assume that the system reaches some maximum entropy state after a violent relaxation phase with strong mixing (LyndenBell 1967). However, the maximum entropy approach is at best partly successful (see, e.g. Yamashiro et al. 1992; Arad & Johansson 2005; Arad & LyndenBell 2005; Yamaguchi 2008; Joyce & Worrakitpoonpon 2011, and references therein) and the only way to improve the results is to introduce additional constraints and ad hoc ingredients (see, e.g. Hjorth & Williams 2010; Pontzen & Governato 2013; Carron & Szapudi 2013). Indeed, relaxation might be incomplete, or a maximum entropy state might not even exist, although one can generalise the concept of entropy by introducing Hfunctions (Tremaine et al. 1986), for which there exist stationary points related to actual stable equilibria even if an entropy maximum does not exist.
Another popular alternative method for understanding the establishment of quasistationary profiles consists in investigating the subspace of selfsimilar solutions (see, e.g. Fillmore & Goldreich 1984; Bertschinger 1985; Henriksen & Widrow 1995; Sikivie et al. 1997; Mohayaee & Shandarin 2006; Alard 2013, but this list is far from exhaustive). While it is difficult to actually demonstrate the onset of selfsimilarity, it seems to be a natural outcome of gravitational dynamics, at least in the absence of any characteristic scale. Although they exist as well in the warm case (see, e.g. Henriksen & Widrow 1995), selfsimilar solutions have been mostly studied in the cold case, for which the phasespace distribution function is of zero initial velocity dispersion. In this configuration, a Ddimensional phasespace sheet evolves in 2D phasespace and builds up a spiral pattern.
The above approaches, along with perturbation theory in cosmological systems (see, e.g. Bernardeau et al. 2002), provide some partial analytical framework to study the VlasovPoisson system. However, in general, these equations usually require a numerical approach, which consists in decomposing the phasespace distribution function on an ensemble of macroparticles interacting with one another with a softened gravitational force (see, e.g. Hockney & Eastwood 1988; Bertschinger 1998; Colombi 2001; Dolag et al. 2008; Dehnen & Read 2011, for reviews on the subject). An alternative way, easily tractable in a small number of dimensions or for systems with a high level of symmetry, consists in using direct Vlasov solvers where the phasespace distribution function is generally sampled on a Eulerian mesh. Most direct Vlasov solvers have been developed in plasma physics and are of semiLagrangian nature. They directly exploit Liouville theorem, namely that the phasespace density is constant along characteristics. In the standard semiLagrangian scheme, a test particle is associated to each grid site where f has to be calculated. This particle is followed back in time during a timestep and the value of f is given by the interpolation of the phasespace density at the previous timestep at the root of the trajectory. In the seminal implementation by Cheng et al. (1976), this is performed in a split fashion between velocities and positions. Many improvements and modifications have been added over time to the splitting algorithm of Cheng & Knorr (see, e.g. the extensive review in the introduction of Sousbie & Colombi 2016), mainly by plasma physicists. The splitting scheme was first applied to astrophysical systems by Fujiwara (1981), Nishida et al. (1981), and Watanabe et al. (1981). In the classical implementation that we use below, reinterpolation of the phasespace distribution function is performed at each timestep using a thirdorder spline.
In all the cases, validating the results obtained from numerical resolution of VlasovPoisson equations remains difficult, particularly if one aims to remain in the mean field limit. In particular, Nbody results are often debated. For instance, the close Nbody encounters and collective effects due to particle shot noise can have some dramatic, possibly cumulative effects (see, e.g. Aarseth et al. 1988; Kandrup & Smith 1991; Boily et al. 2002; Binney 2004; Joyce et al. 2009; Colombi et al. 2015; Beraldo e Silva et al. 2017; Romero & Ascasibar 2018), particularly when initial conditions are cold or close to cold (see, e.g. Melott et al. 1997; Melott 2007). While Vlasov codes do not use particles, they are still subject to nontrivial numerical effects, because a phasespace grid still remains a discrete representation of the system (see, e.g. Colombi et al. 2015; hereafter C15). But since the numerical implementation is still different from Nbody codes, a comparison between Vlasov and Nbody codes seems appropriate and timely, especially when trying to analyse in detail the quasistationary state reached in the fluid limit by gravitational systems.
One question indeed remains open. What are the main processes involved in the violent relaxation phase leading to a quasistationary state? We propose here to approach this question by studying the evolution of a number of initially spherical systems with various initial density profiles and velocity dispersions, focusing on the phasespace structure. Spherical symmetry allows us to compare highresolution Vlasov simulations to Nbody simulations. Our analyses focus on the detailed structure of the phasespace distribution function and comparisons with predictions from selfsimilarity.
The advantage of systems with initial spherical symmetry is that they have been studied in great detail in the literature both from the theoretical and the numerical points of view. One major question for instance concerns the role of radial instabilities and radial orbit instability in the establishment of the quasisteady state observed after violent relaxation. Some spherical equilibria or quasisteady states are known to be unstable to radial perturbations (see, e.g. Hénon 1973; Henriksen et al. 1997) as well as angular perturbations that translate into radial orbit instability as well studied in the literature (see, e.g. Polyachenko & Shukhman 1981; Merritt & Aguilar 1985; Barnes et al. 1986; Cannizzo & Hollister 1992; Barnes et al. 2009; Maréchal & Perez 2011; Vogelsberger et al. 2011; Polyachenko & Shukhman 2015). These perturbations are usually induced by shot noise due to the discrete nature of the distribution of stars or particles in the system, which directly relates to the discussion above about the validity of numerical simulations. Here, it is interesting to see what happens in the mean field limit, or at least in a regime that tries to approach it by using a direct Vlasov code and Nbody simulations with a very large number of particles.
More specifically, assuming G = 1 and following the footsteps of Burkert (1990), Hozumi et al. (1996), and Hozumi et al. (2000), we perform a number of controlled numerical experiments of unity total mass systems, initially spherical with a powerlaw density profile and a Gaussian isotropic velocity dispersion:
with R_{0} = 2 being the initial radius of the sphere and the initial slope spanning the range n = 0 to n = −1.5. We consider “warm” and “cool” cases defined by their respective values of the virial ratio, η = 0.5 and 0.1, with
where T is the total kinetic energy and W the total potential energy of the system.
In these simulations, we aim to study in detail the evolution of the phasespace distribution function, the onset of instabilities, and the consequence of these at the coarsegrained level. We also relate our measurements of fine details of the spiral pattern of the phasespace distribution to expectations from selfsimilar dynamics (see, e.g. Alard 2013; hereafter A13). To perform our simulations, we use three kinds of codes: a spherical semiLagrangian Vlasov solver, VlaSolve, presented in C15, the Nbody public treecode Gadget2 (Springel 2005), and a standard spherical shell Nbody code (see, e.g. Hénon 1964). Importantly, while the systems are forced to remain spherical in VlaSolve and in the shell code, it is not the case for Gadget2, which allows for the development of angular anisotropies. The variety of the codes employed in this work helps us to understand in detail the nature of different sources of instabilities, whether physical or numerical. In particular, we study the influence of finite spatial resolution in Vlasolve and of finite number of particles in the shell code and in Gadget2.
This paper is organised as follows. In Sect. 2, we present the numerical codes used to perform the simulations and provide details on the various runs we performed. In Sect. 3, we perform a detailed visual inspection of the phasespace distribution function and discuss the various dynamical phases at play. Subsequently, Sect. 4 deals with selfsimilarity: we show how the calculations of A13 can be extended to spherical systems in a very simple way, and compare theoretical predictions on the shape of the phasespace distribution function to our numerical experiments. Finally, Sect. 5 summarises and discusses the results.
2. Simulations
In spherical symmetry, Vlasov equation can be written as
with r being the spherical radius, υ_{r} the radial velocity, j the conserved angular momentum and M(< r) the mass contained in a sphere of radius r. For a given value j, the evolution is therefore driven by the interplay between the gravitational force, dominating at large radii, and the centrifugal force j^{2}/r^{3}, dominating at small radii.
To solve Eq. (6), we have resorted to two numerical methods.
Firstly and mainly, we employ the spherical Vlasov solver VlaSolve presented in C15. This semiLagrangian code is similar to that of Fujiwara (1983) and thus uses the splitting algorithm of Cheng et al. (1976) to compute the evolution of the phasespace density on a mesh. The phasespace is divided into threedimensional (3D) cells along radius, radial velocity, and angular momentum. A logarithmic scaling is used for radius to properly resolve the dynamics of the collapse at low radii, making the use of a minimum radius R_{min} necessary. To compute accurately the dynamics at low radii, the exact time spent by matter elements inside the sphere of radius R_{min} is computed assuming that gravitational force is negligible, which is an improvement over previous implementations which used the reflecting sphere method (see, e.g. Gott & Thuan 1976; Fujiwara 1983). More algorithmic details and tests of the code can be found in C15. Unless specified below, the grid for all the simulations is such that (N_{r}, N_{v}, N_{j}) = (2048, 2048, 128), where N_{r} is the number of vertices of radius in log scale, N_{v} is the number of vertices of radial velocity in linear scale, and N_{j} is the number of slices of angular momentum such that the kth slice contains fluid of angular momentum j_{max}((k−1/2)/N_{j})^{2}, corresponding formally to the interval [j_{max}((k−1)/N_{j})^{2}, j_{max}(k/N_{j})^{2}]. The computation domain is log_{10} R_{min} ≡ −2 < log_{10}(r)< 1.4, −υ_{l} < υ_{r} < υ_{l} and 0 < j < 1.6, where υ_{l} = 2 for the runs with a virial ratio η = 0.5 and υ_{l} = 3 for the runs with a virial ratio η = 0.1. The limits are chosen such that almost all the mass of the system is contained in the computing domain during the simulated time, except for matter elements passing inside the sphere of radius R_{min} or those escaping from the system at large radius. For all the simulations, the timestep was chosen to be constant, equal to Δt = 0.005, a value larger than in C15 to avoid excessive diffusion due to overfrequent resamplings of the phasespace distribution function, but we checked it is still on the safe side.
To avoid excessive aliasing effects in the VlaSolve runs, we apodize the initial profile given by Eq. (4) as follows,
with Δ = 1/2, exactly as in C15. This apodization slightly changes the value of the virial ratio, of the order of 10% at most.
The second code we use relies on the standard spherical shells approach as in for example Hénon (1964) and Gott & Thuan (1976). We note that we employ it only for the most critical cases, namely (η, n) = (0.1, −1) and (η, n) = (0.1, −1.5), when the results from Gadget2 differ too much from VlaSolve. In this Nbody code, each particle represents a shell in configuration space and interacts with the other particles through gravitational force, −GM(< r)/r^{2}, which can be obtained very easily with a sorting procedure. The resolution of the Lagrangian equations of motion of the particles is performed simply with a leapfrog integrator with a constant timestep Δt = 0.001, five times smaller than the one chosen for the Vlasov code. Similarly as in Vlasolve, the Leapfrog algorithm is implemented using a decomposition of the Hamiltonian of the motion into a fully analytical drift part including centrifugal force and a kick part including solely gravitational force (see, e.g. Colombi et al. 2008). Initial shells distribution simply consists in taking the initial conditions of the Gadget2 simulations described below, with the radius of each shell being equal to the magnitude of the position of each particle and their respective radial velocities and angular momenta directly derived from the three coordinates of the particle velocities.
Finally, we perform simulations using the public 3D Nbody treecode Gadget2 (Springel 2005) in its noncosmological setup and with the treecode part only. The positions and velocities of particles are generated in a random way using a standard rejection method.
In the Gadget2 simulations, spherical symmetry is no longer imposed, which leaves room for the development of angular anisotropies due to Poisson fluctuations in the initial particle distribution, even though the profile is initially spherical in the statistical sense. In particular, we see that the (η, n) = (0.1, −1) and (η, n) = (0.1, −1.5) simulations are subject to radial orbit instability. All the simulations in this work involved 10 million particles, except for one additional one performed with 100 million particles to examine the case (η, n) = (0.1, −1.5) more closely. The parameters for the Gadget2 runs, in terms of softening length, force accuracy, and timestep control are otherwise the same as in C15.
Table 1 summarises the parameters used to perform the simulations of this work. In particular, the last two columns give the value of 25 t_{dyn} where
corresponds to the duration of a full radial orbit in a harmonic potential corresponding to a fixed density ρ. We estimate ρ from the simulations themselves, once the system has reached a quasistationary state, either directly at the centre of the system or as the average of the density in a sphere containing 10% of the total mass. Our estimates are rather crude but justify the choice of final time equal to 80 and 40 for the “warm” and “cool” cases, respectively. However, we must remain aware of the fact that increasing the magnitude of the slope n decreases the value of t_{dyn}.
Parameters used for the simulations performed in this article.
3. Visual inspection: phasespace structure and density profiles
Figures 1 and 2 display, for a typical slice of fixed angular momentum, the phasespace distribution function of the “warm” and “cool” VlaSolve simulations, respectively. Additionally, Fig. 3 examines in more detail the case (η, n) = (0.1, −1.5), which is subject to radial orbit instability (ROI), while Fig. 4 provides projected density profiles. To supplement our discussion about ROI, we study in Fig. 5 the velocity anisotropy parameter for all the simulations as well as deviation from sphericity for the runs which experience ROI.
Fig. 1. Snapshots of the phasespace density for the simulations with “warm” initial conditions, η ≃ 0.5. A typical slice of f(r, υ_{r}, j) with j = 0.16 is shown in (r, υ_{r}) space for the VlaSolve simulations at t = 0, at an early time t = 15, at an intermediate time used to perform tests of selfsimilarity of the phasespace spiral and at the final time. 
Fig. 2. As in Fig. 1 but for the simulations with “cool” initial conditions, η ≃ 0.1 and for a slice with j = 0.06. 
Fig. 3. Closer examination of the onset of instabilities in phase space for the (η, n) = (0.1, −1.5) simulations: effects of spatial and mass resolutions. For the same angular momentum slice, j = 0.06, as in lower panels of Fig. 2, the phase space density is represented in (r, υ_{r}) space. The two first lines of panels correspond to two VlaSolve simulations with respective resolutions (N_{r}, N_{v}, N_{j}) = (2048, 2048, 128) and (1024, 512, 512). The next two lines of panels correspond to two Gadget2 simulations with respective numbers of particles N = 10^{7} and N = 10^{8} and the last line of panels gives, for N = 10^{7}, the result obtained for the shell code. For the Nbody simulations, the phasespace density is sampled on grids with resolution (N_{r}, N_{v}, N_{j}) = (1024, 1024, 32). One can notice that the phasespace sheet is less well defined or “fuzzier” in the Nbody simulations than in the Vlasov code at low radius; this is because what is actually plotted is the distribution of particles (or shells) in a relatively large interval of angular momentum j ∈ [0.056, 0.077] to have sufficient number of particles to trace the phasespace distribution function, while for the Vlasov simulation, we selected only the slice corresponding to the value of j of interest. This figure illustrates the effect of radial instabilities and their dependence on spatial resolution (for the Vlasov code) or mass resolution (for the Nbody code). We note the close agreement between the shell simulation and the Gadget2 runs with 10 million particles in the middle panels of the third and fifth lines, while they diverge in the right panels when radial orbit instability effects become prominent in Gadget2. 
Fig. 4. Radial density profile measured at various times in the VlaSolve simulations, namely initial conditions (blue dashes), intermediate time used to perform tests of selfsimilarity of the phasespace spiral (green), and final time (red dashes). For the final time, the results are also compared to a Gadget2 run (thick grey), as well as the output of the shell code (thick grey dashes) for (η, n) = (0.1, −1) and (0.1, −1.5). In addition, the logarithmic slopes −4 and −2.1 (as measured in Hozumi et al. 2000) are shown respectively as a thin solid and a thin dashed line. 
Fig. 5. Velocity anisotropy and deviation from sphericity. Left panel: Velocity anisotropy parameter as a function of time for the VlaSolve (solid lines) and Gadget2 simulations we performed (dashed curves). Right panel: Evolution of the departure from spherical symmetry for the two kind of initial conditions experiencing radial orbit instability in Gadget2, namely (η, n) = (0.1, −1.0) (blue curves) and (η, n) = (0.1, −1.5) (red curves, solid and dashed for the 10 and 100 million particles simulations, respectively). The quantity indicates the deviation from unity of r_{axis} (in log scale and with a negative sign for r_{axis} > 1), where r_{axis} = b/c (upper curves) or b/a (lower curves), and a ≤ b ≤ c are the principal axis lengths of the Gadget2 particle distribution derived from the inertia tensor. 
Thanks to the high resolution of our simulations, when examining these figures, one can clearly separate, for the first time, two or three wellknown dynamical phases, depending on initial conditions: (i) a violent relaxation phase during which the system converges to a quasisteady state by building a very regular spiral structure in phasespace, (ii) a quiescent phase during which the quasisteady state is preserved against smallscale radial instabilities which can destroy the spiral and (iii) relaxation to a nonspherical state through ROI when the system is prone to develop it. The novelty in our measurements is obviously not the discovery of the various phases of the dynamics, which are heavily discussed in the literature, but instead the clear articulation between them for the systems we study. We now discuss these three phases in detail.
3.1. Violent relaxation
In a first phase, the system undergoes violent relaxation that leads quickly to the establishment of a quasisteady state. During this phase, spherical symmetry is preserved and the phasespace distribution function presents a very regular spiral structure in all cases, even in the Nbody runs, thanks to the large number of particles we used to perform them. A movie of the evolution of the system shows that it is also subject to a global pulsation that introduces at some point irregular features in the phasespace distribution function in the outer parts of the spiral and at large radii, for example the dark region in upperright panel of Fig. 1. During this violent relaxation phase, Gadget2 agrees very well with VlaSolve, even in regions where these irregular features develop, as already noticed by C15 for the (η, n) = (0.5, 0) case, which shows that these features are intrinsic to the physical system and are not related to some additional instability due to some numerical noise.
3.2. Quasisteady regime with smallscale radial instabilities
In a second phase, the system stays in quasiequilibrium and preserves its spherical symmetry. However, some radial instabilities perturb it at small scales, irrespective of the numerical technique used, and can destroy the spiral structure. The time of the appearance of these instabilities is related to spatial resolution in the VlaSolve simulation and to the number of particles in the Nbody simulations. This is well illustrated by Fig. 3 for (η, n) = (0.1, −1.5): for instance, increasing the number of angular momentum slices in the VlaSolve simulation reduces the magnitude of the perturbations of the spiral (compare middle insert of first and second line of panels), and similarly when increasing the number N of particles in the Gadget2 runs (compare middle insert of third and fourth line of panels).
In the Nbody case, these collective instabilities are induced by small, random but correlated errors on the gravitational force due to Poisson fluctuations in the particle distribution. In the VlaSolve code, they are related to coherent errors on the force due to the representation of the phasespace density on a grid, but the effect is analogous to the Nbody case. These instabilities become naturally more significant when the initial velocity dispersion is reduced, as explained in C15. We also notice here that they take place earlier for larger n, in agreement with our calculation of dynamical times in the two right columns of Table 1. As a result, during the interval of time covered by our simulations, they can be seen in the cool runs and in the warm case for (η, n) = (0.5, −1.5), but they are not present or are negligible in other cases.
An important observation is that these instabilities intervene only at the fine level: they do not change the structure of the system at the coarse level, even quantitatively. To be more specific, if the phasespace density was smoothed at scales larger than the interfilament separation – by filament, we mean for example some fold of the spiral structure – and larger than the size of the fluctuations introduced by radial instabilities, there would be no significant difference between late times, where these instabilities can destroy the spiral structure, and earlier times, when the spiral structure is still well defined. A good way to illustrate this consists in examining the projected density profile measured in the VlaSolve simulations, as displayed in Fig. 4, and to compare red curves to the green ones, which correspond to these aforementioned late and earlier times, respectively. The calculation of the projected density, by integrating the phasespace density over velocities, indeed corresponds to some coarsegraining procedure, although such anisotropic smoothing does not erase the quasicaustic structures seen on the green curves of Fig. 4. These bumps correspond to projection of parts of the spiral (or any filament) that are orthogonal to configuration space. However, with proper (adaptive) smoothing at scales larger than the space between successive spiral folds, one can be convinced that agreement between the green curves and the red curves, already very good in most cases, should improve further.
The quasiequilibrium built cinematically by the spiral is therefore stable against radial perturbations, but not necessarily the spiral structure.
3.3. Deviation from sphericity: radial orbit instability
In a third phase, small angular anisotropies induced by numerical noise become amplified through ROI for n ≤ −1 in the cool cases and the system deviates from spherical symmetry by acquiring a prolate shape (right panel of Fig. 5). A consequence of ROI is the reduction of the magnitude of the spherically averaged density profile ρ(r) at small radii, as can be seen in the two bottomright panels of Fig. 4.
Radial orbit instability signature is best seen in the velocity anisotropy parameter
(see, e.g. Hozumi et al. 1996), where υ_{r} and υ_{⊥} are respectively the radial and transverse velocities, as plotted in left panel of Fig. 5. Due to the dominant nature of radial infall during the very first phase of violent relaxation, cool initial conditions induce a strong velocity anisotropy, which is known, when exceeding some (still not fully known) threshold, to trigger ROIs in the presence of small perturbations to spherical symmetry (see, e.g. Polyachenko & Shukhman 1981; Merritt & Aguilar 1985; Barnes et al. 1986, 2009; Maréchal & Perez 2011; Polyachenko & Shukhman 2015, and references therein). In this case, the onset of ROI significantly reduces the value of α, as illustrated by the left panel of Fig. 5 for (η, n) = (0.1, −1) and (η, n) = (0.1, −1.5).
In our Gadget2 simulations, small perturbations from spherical symmetry are related to shot noise, meaning that the onset of ROI is particle number dependent (see, e.g. Benhaiem et al. 2018), as illustrated in Fig. 5 by our two runs with 10 and 100 million particles in the (η, n) = (0.1, −1.5) case. On this figure, one also notices that ROI takes place later for n = −1 than for n = −1.5, but this is roughly consistent with the dynamical times given in the last two columns of Table 1.
The conditions of establishment of radial orbit instability are however not yet fully understood: some theoretical calculations and numerical experiments show that it should take place when α > α_{critical} with α_{critical} ranging between 1 and 2.9 (see, e.g. Polyachenko & Shukhman 2015): this condition is clearly satisfied for (η, n) = (0.1, −1.0) and (0.1, −1.5) when examining the left panel of Fig. 5. Strictly speaking, given the limited amount of time we run the simulations, the other cases remain undecided even though we do not detect any ROI. The results obtained elsewhere in the literature, in particular by Merritt & Aguilar (1985), and Barnes et al. (2009), suggest that our “warm” systems are probably not prone to ROI, while, for (η, n) = (0.1, −0.5) and (η, n) = (0.1, 0), there is still a chance that ROI develops after some time. Clearly, our simulations are not run long enough to have all the details of the history of the system, which might evolve further to another interesting state.
Finally, we note that whether or not pure radial instability takes place before ROI is difficult to decipher from our simulations. Using linear analysis during collapse phase, Aarseth et al. (1988) argue that in the cold case, angular anisotropies introduced by Poisson noise are subdominant compared to radial ones, which suggests that a radial instability phase could take place before ROI. This argument is partly supported by Fig. 3, where excellent agreement is found between Gadget2 and the spherical shell code at t = 4 (middle inserts of third and fifth lines of panels), when radial perturbations already significantly disrupt the phasespace spiral while ROI has not yet developed.
4. Selfsimilarity in phase space
In practice, seeking selfsimilar solutions to VlasovPoisson equations consists in finding solutions invariant with respect to some homothetic transforms, for example,
which requires in this case f to be of the form of
Solving the Vlasov equation provides the solution for the function F. While a rigorous framework can be set to define selfsimilarity through Lie derivatives (Carter & Henriksen 1991), there is no unique way to express it. For instance, as studied in for example Henriksen & Widrow (1995), instead of Eq. (10), one can, in the spherically symmetric case, introduce anisotropy in the selfsimilar solution by separating the radial velocity from the angular momentum, or, in an extreme but wellknown case, one can simply assume, as in the cold case, pure radial motions with zero angular momentum.
Strictly speaking, selfsimilarity implies a pure powerlaw behaviour for the projected density (see, e.g. Henriksen & Widrow 1995), which is obviously not the case for the systems we study here when examining Fig. 4, except to some extent at sufficiently large radii. However, selfsimilarity usually takes place only in a limited domain of phasespace. This can for instance simply be due to the finite extension of the system, which may induce deviations of the density profile from a pure power law, even though phasespace density is perfectly self similar. As a clear illustration of this statement, Alard (2013) argues that despite the cutoffs due to the finite nature of the system, local selfsimilarity in phasespace implies a powerlaw behaviour for the quantity Q(r) = ρ(r)/σ(r)^{3} where ρ(r) and σ(r) are respectively the projected density and the velocity dispersion, even if these latter are not found to be exact power laws of radius due to the cutoffs. Such a powerlaw property for Q(r) is verified to a great accuracy by darkmatter halos (see, e.g. Taylor & Navarro 2001), of which the density profiles are known to deviate from pure power laws (Navarro et al. 1996, 1997). As shown by Dehnen & McLaughlin (2005), solving the spherical Jeans equation assuming that Q(r) is a power law can indeed lead to nonpure powerlaw density profiles.
Here, we are clearly not in the case of a single selfsimilar regime (with possible cutoff effects). Indeed, for each value of the angular momentum j, our systems can trivially be separated into two distinct regions of phase space. In the first region, corresponding to a smallenough radius, the centrifugal acceleration j^{2}/r^{3} dominates. Its powerlaw nature is expected to induce selfsimilarity in some domain of the considered phasespace slice. In the second region, corresponding to a largeenough radius, the gravitational acceleration −GM(< r)/r^{2} dominates. Provided that it is also a powerlaw of radius, one expects another selfsimilar behavior. Given the discussion above about the possible effects due to the finite extent of the system, selfsimilar properties may be found even if the force is not exactly a power law. The transition between these two regions is sharp, as illustrated by Fig. 6, which is important to make our approach meaningful.
Fig. 6. Small and large radii regimes in the (η, n) = (0.5, −1) VlaSolve run at t = 40. Left panel: plots separately the magnitude of centrifugal force j^{2}/r^{3} and of the gravitational force as functions of radius as well as the magnitude of the sum of both forces. Right panel: isocontours of the specific energy (white curves) superposed on the phasespace distribution function for j = 0.16. In addition, the white dashes and dotdashes assume that only the gravitational or the centrifugal force contributes, respectively. 
Selfsimilarity also predicts the formation of a spiral in phasespace, of which the structure is defined by the selfsimilar parameters (see, e.g. Fillmore & Goldreich 1984, A13). We note, as already mentioned, that the onset of selfsimilarity does not require the assumption of cold initial conditions as is often supposed. For instance, in the calculations of A13, no such hypothesis is made, and the existence of a spiral structure in phase space is clearly seen when assuming selfsimilarity. Here, we cannot rigorously demonstrate the existence of such a spiral structure but can postulate it, predict its local properties in each of the supposed selfsimilar regimes mentioned above, and compare the predictions to our simulation measurements, which we do below, closely following Alard (2013; and priv. comm.).
When examining a slice of fixed angular momentum, we notice that the Vlasov equation for a spherical system is exactly analogous to the 1D case:
except that the force is derived from the following scalar field
where ϕ(r) is the gravitational potential. A13 derived detailed selfsimilar solutions in the 1D case that we extend below in the regimes where the centrifugal force dominates ψ(r)≃j^{2}/(2r^{2}) and in the regime where gravitational potential dominates and is a power law, ψ(r)≃ϕ(r)∝r^{β + 2}.
Assuming that the conserved angular momentum j is a dummy variable, the selfsimilar solution for the phasespace distribution function can be expressed as follows:
Setting
we can express both the gravitational and the centrifugal force as functions of these new variables.
Starting from the gravitational force, we define a function U such that
and we assume that
with
Then,
and we obtain, from the Poisson equation,
Similarly, the centrifugal force can be written
which implies
and, if we assume again that , we obtain
By injecting these various expressions in the Vlasov equation in a regime where either the gravitational or centrifugal force dominates, one obtains
Eliminating time dependence in this equation imposes
Enforcing stationarity of the force t^{α}U(r/t^{α1}) with the powerlaw (19) implies
Hence the only viable solution is
which is of course consistent with Eq. (23) and leaves α_{0} as a free parameter if the centrifugal force dominates; while if the gravitational force dominates, it fixes α_{0} = −(β + 2)/β. We note that, in general, total mass is not conserved. Indeed, enforcing total mass conservation (or mass conservation per angular momentum slice, as well) imposes α_{0} + α_{1} + α_{2} = 0, a condition which is fulfilled only for β = −3/2 and in this case α_{0} = 1/3. This is not a problem because selfsimilarity is expected to take place only in a finite dynamical range.
To follow the notations of A13 as closely as possible, we now make the following change of variables,
Equation (26) becomes
which is exactly the same equation as Eq. (11) of A13, except that the first term α_{2} + 2 is replaced here with −α_{0}. Hence, the solution of this equation is very similar to the expressions given in A13. The main difference here is that the values of β we consider are outside the domain of validity of the calculations of A13, which implies that the isocontours of the solutions are closer to hyperbolic curves than to a spiral. Here, however, we have to take into account the fact that we have two distinct supposed selfsimilar regimes, one dominated by centrifugal force, say for r ≲ r_{crit}(j), and the other dominated by gravitational force, say for r ≳ r_{crit}(j). Hence, the actual solution is the connection between too partial solutions following selfsimilar properties.
Rescaling variable so that k = 1/2 in Eq. (18), and introducing, exactly as in A13, the new variables
we obtain nearly exactly Eq. (12) of A13, but the parameters of this equation change according to whether the value of R cos Ψ is above or below a threshold fixed by r_{crit}. With we write, following exactly the footsteps of A13, the general solution for H when the powerlaw force is stationary,
with Ψ ∈ ] − π/2, π/2[ and where Q is some function. At this point, introducing the same concept of spiral as in A13 is not simple, because u > 0 does not allow Ψ to make a full excursion on the circle. Furthermore, the values of the logarithmic density profile slope β that we have to consider range in the interval −4 ≲ β < −2, which implies, from Eq. (30), −1/2 ≲ α_{2} < 0, hence some divergence of the integral
when Ψ approaches π/2. This is however not a real problem, because the objective is to connect two selfsimilar solutions. Here, we are unable to demonstrate the existence of the spiral structure in phasespace; we must postulate it. We therefore define a new angular variable θ and
where g(θ) is a function with period 2π verifying
and the symbols − and + correspond to the regimes dominated by the centrifugal and the gravitational force, respectively. Function g(θ) makes a smooth transition between g_{−} and g_{+}. The only, trivial but important fact we have to know is that I(θ) defined this way is roughly proportional to θ which allows us now to explicitly define the concept of a spiral across both selfsimilar domains. Interestingly, the subsequent calculations of A13 are not changed at all when taking this new definition of I, and Eq. (19) of A13 still stands in each selfsimilar domain, with now an unknown constant instead of a welldefined integral as in A13.
Therefore, in the situation where there are many folds, we have the following expected relationship for the interfold distance in each selfsimilar region:
In particular, coming back to standard variables (r, υ_{r}), the interfold distances along the axis υ_{r} = 0 read
To test this property directly, we analyse, at a time where the spiral structure is still well defined, the function f(r, υ_{r} = 0, j) for a fixed value of angular momentum, as illustrated by Fig. 7. We determine the positions of the folds using local parabolic fits. For each fold i, we determine two semiheight radial positions log(r_{l, i}) and log(r_{r, i}) (the computational grid being logarithmic in radius) on the left and right of the peak, respectively, and define the error on the position of the peak as δlog(r_{i}) = log(r_{r, i})−log(r_{l, i}).
Fig. 7. Illustration of the method used to determine the positions of the folds and corresponding interfold distance law at small radius in a phasespace slice. Top panel: zoom performed around the axis υ_{r} = 0 in the region dominated by centrifugal force for the (η, n) = (0.5, −1.0) simulation at t = 40. The corresponding phasespace distribution function f(r, υ_{r} = 0, j = 0.16) is plotted in the lower panel. The black dots give the positions of local maxima estimated with our local quadratic fit, while the blue and red dots provide upper and lower bounds to compute the (very conservative) error bars shown on Fig. 8. 
Fig. 8. Interfold distance at null radial velocity vs. selfsimilar predictions for a fixed value of angular momentum. The distance dr between local maxima of the function f(r, υ_{r} = 0, j) is plotted as a function of r for the VlaSolve runs with various initial conditions. The time considered corresponds to the middle column of panels in Figs. 1 and 2. The two left columns of panels correspond to the “warm” case, η = 0.5 with j = 0.16, and the two right ones to the “cool” case, η = 0.1 with j = 0.06. Then, odd column numbers (1 and 3) and even column numbers (2 and 4) correspond to the regime where centrifugal/gravitational force dominates, respectively. On each panel a red line indicates the logarithmic slope predicted by selfsimilarity. When the gravitational force dominates, two additional dashed curves provide a bracket of the estimated slope taking into account deviations from selfsimilarity, i.e. variations of the effective logarithmic slope of the gravitational force. 
Figure 8 summarises the results of our measurements in the VlaSolve simulations. At a “small” radius, the system is dominated by the centrifugal force, β = −4, hence dr ∝ r^{3}. This prediction is compared to measurements in the simulations in the first and third columns of Fig. 8, which correspond to the simulations with “warm” and “cool” initial conditions, respectively. At a large enough radius, where the system is dominated by gravity, the force is only approximately a power law but an average slope can nevertheless be inferred in some interval of scales E_{grav}, corresponding to the regime where the gravitational force remains at least ten times larger than the centrifugal force and for r smaller than the turnaround radius. In the second and fourth columns of panels of Fig. 8, the slope of the red line is given by the corresponding value of 1 − β/2. There are also two dashed cyan and green lines corresponding to the minimum and maximum value of β found in E_{grav}, which gives an idea of deviation from a pure power law. Globally, the simulations agree rather well with selfsimilar predictions, except maybe at very small r in the first and third columns of panels of Fig. 8 and in the top panel of the second column in the same figure. We note however that measurements of the interfold distance for small values of r might be partly spurious, because we are in a regime where the phasespace distribution function is small and can be affected by aliasing. Also, we note that the spiral structure survives only a short time for (η, n) = (0.1, −1.5) which leaves only a small number of folds to deal with. Nevertheless, the agreement with selfsimilarity remains good when taking into account the limitations found in all the cases at very small r, already after only a few dynamical times.
Another interesting property than can be derived directly from Eq. (39) is the local shape of the spiral near the axis υ_{r} = 0, hence Ψ ≃ 0, and for small r, hence small R. Following the unnumbered equation after Eq. (16) of A13 and taking into account the fact that α_{2} < 0, we expect an isocontour of the function H(R, Ψ) to have the following shape in the regime R ≪ 1, Ψ ≃ 0,
We note therefore that because of the form of the interfold law (44), the spiral actually locally coincides with a curve defined by R ∝ θ^{−α2} with θ playing the same role as Ψ, but no longer restricted to ] − π/2, π/2[, belonging instead to this interval and its multiples ] − π/2 + 2kπ, π/2 + 2kπ[. The left panel of Fig. 9 nicely illustrates how this prediction closely matches the local spiral shape of the simulated phasespace distribution function, including its local curvature, in the regime dominated by the centrifugal force (β = −4) for (η, n) = (0.5, −1.0).
Fig. 9. Spiral shape vs. selfsimilar predictions for the VlaSolve run with (η, n) = (0.5, −1.0). Left panel: comparison, at t = 40, of the local shape of the spiral structure in the region dominated by centrifugal force to the curve given by Eq. (46) (red curve). Middle panel: test for selfsimilarity in time. The green spiral shape obtained at t = 30 is rescaled according to Eqs. (47) and (48) to be compared to the blue one, in the regime dominated by centrifugal force. Right panel: using only the determination of the position of the folds in the region dominated by centrifugal force, it is possible to draw the full shape of the spiral if the gravitational potential is known by interpolating the specific energy in angle coordinate defined by Eqs. (49)–(51). 
In middle panel of this figure, we also check, in the regime dominated by centrifugal force, for selfsimilarity in time of the spiral shape, namely that if one considers the system at two different times, t_{1} and t_{2}, the state at t = t_{2} should superpose to the state at t = t_{1} rescaled as follows
Of course, since we have two distinct selfsimilar regimes, this property works well only in the neighbourhood of υ_{r} ≃ 0 and for values of r where the gravitational force is subdominant compared to the centrifugal force.
Finally, the right panel of Fig. 9 shows that if gravitational potential is known, the spiral shape of the phasespace distribution function can be fully reconstructed accurately simply by knowing its intersection with the υ_{r} = 0 axis in the regime dominated by angular momentum (or reversely, in the regime dominated by gravitational force) by simple linear interpolation of the specific energy E along the spiral during an orbit in the following angle coordinate 𝒜,
where s is a curvilinear coordinate. This technique was actually used to draw the spiral of the middle panel. Of course, this result is kind of trivial from a dynamical point of view. However, it suggests that passing to actionangle space or energyangle as performed here may represent the best way to smoothly connect both selfsimilar regimes and therefore to have a full analytic description of the fine grained structure of the phasespace distribution function. To do this, one needs to relate locally the angular variable Ψ intervening in the selfsimilar solutions to the angle given by Eq. (49), but this is left for future work.
5. Conclusion
In this article, we carried out a detailed analysis of the phasespace structure of various systems with spherical initial conditions, consisting in a powerlaw density profile with a Gaussian velocity dispersion. Two cases were considered: the “warm” setup with virial ratio η ≃ 0.5 and the “cool” one with η ≃ 0.1. The choice of such initial conditions is not really new but the numerical setup is different from what can be found in the literature. Firstly, we compared three kind of codes, a Vlasov code, a treecode, and a shell code. Secondly, we performed this comparison with unprecedented numerical resolution, namely (N_{r}, N_{v}, N_{j}) = (2048, 2048, 128) for the Vlasov code and 10 million particles for the Nbody simulations, up to 100 million for one Gadget2 run.
The high resolution of our simulations allowed us to study all the fine details of the phasespace distribution function and to really distinguish, for the first time, three wellknown dynamical phases of the evolution of these systems, namely, (a) a violent relaxation phase to a quasisteady state where the phasespace density can be mainly described by a smooth spiral structure winding with time, (b) a steadystate phase during which radial instabilities can destroy the spiral structure but do not affect the macroscopic properties of the system, and (c) relaxation to a nonspherical state due to radial orbit instability in the cool cases with n ≤ −1. Obviously, we did not push the simulations long enough to approach the socalled gravothermal catastrophe regime, where a corehalo structure can appear due to collisional relaxation (see, e.g. Antonov 1962; LyndenBell & Wood 1968; LyndenBell & Eggleton 1980; Makino 1996; Baumgardt et al. 2003).
While the concept of violent relaxation phase to a quasisteady state is a wellknown process studied heavily in the literature, the fact that it is expressed as a welldefined smooth spiral structure in phase space is not trivial. Subsequent radial instabilities that can appear indeed do not introduce sufficient disorder to significantly disturb the steady state initially built from the kinematic evolution of the spiral structure. Only radial orbit instability changes the properties of the system at the coarse level. But even in this case, this happens only at small radii, the outer part of the system still being given by the quasisteady state solution obtained previously.
These results seem to diverge from what can be obtained in the pure cold case. For instance, a similar analysis was done by Henriksen et al. (1997), but for spherical initially cold systems with powerlaw density profiles using a shell code: in this case, Henriksen et al. (1997) found that radial instabilities are sufficiently strong to destroy the quasisteady selfsimilar state obtained during the violent relaxation phase and produce a density profile close to ρ(r)∝r^{−2}, hence, changing the properties of the system at the macroscopic level. However, the number of shells employed by these authors was rather small and it is possible, despite the convergence tests they did, that they missed an intermediary phase where radial instability would be sufficient to destroy the spiral while preserving, as in our case, the coarse grained properties of the distribution function. We note that this question remains rather academic as radial orbit instability is expected to be prominent for such systems when allowed to deviate from spherical symmetry, although they still present some selfsimilar properties (Vogelsberger et al. 2011).
In a second part of our analysis, in order to at least partly understand the dynamical processes taking place during the violent relaxation phase, we examined the properties of the spiral structure of our systems in the framework of selfsimilar solutions. Obviously, our systems are not fully selfsimilar, but we showed that they follow selfsimilar properties in welldefined domains of phase space. Indeed, each slice of phase space of given angular momentum can be trivially decomposed into two regions, one where centrifugal force dominates, and the other one where gravitational force dominates. While the centrifugal force j^{2}/r^{3} is a pure power law, this is only approximately the case for the gravitational force. Nevertheless, this approach allowed us to partly predict the properties of the spiral structure, for instance the interfold distance at zero radial velocity. Even if that is not enough by itself to be able to fully predict the steady state abinitio, selfsimilarity in phase space remains an interesting path of investigation.
Acknowledgments
We thank C. Alard for providing us the main insights about the analyses performed in Sect. 4, T. Sousbie and S. Hozumi for stimulating discussions. This work was supported in part by ANR grant ANR13MONU0003. We also acknowledge the support of YITP in organising the workshop “VlasovPoisson: towards numerical methods without particles” in Kyoto, funded by grant YITPT1502, ANR grant ANR13MONU0003 and by Institut Lagrange de Paris (ANR10LABX63 and ANR11IDEX000402).
References
 Aarseth, S. J., Lin, D. N. C., & Papaloizou, J. C. B. 1988, ApJ, 324, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Alard, C. 2013, MNRAS, 428, 340 [NASA ADS] [CrossRef] [Google Scholar]
 Antonov, V. A. 1962, Vest. Leningrad Univ., 7, 135 [NASA ADS] [Google Scholar]
 Arad, I., & Johansson, P. H. 2005, MNRAS, 362, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Arad, I., & LyndenBell, D. 2005, MNRAS, 361, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, J., Hut, P., & Goodman, J. 1986, ApJ, 300, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, E. I., Lanzel, P. A., & Williams, L. L. R. 2009, ApJ, 704, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H., Heggie, D. C., Hut, P., & Makino, J. 2003, MNRAS, 341, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Benhaiem, D., Joyce, M., SylosLabini, F., & Worrakitpoonpon, T. 2018, MNRAS, 473, 2348 [NASA ADS] [CrossRef] [Google Scholar]
 Beraldo e Silva, L., de Siqueira Pedra, W., Sodré, L., Perico, E. L. D., & Lima, M. 2017, ApJ, 846, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rev., 367, 1 [Google Scholar]
 Bertschinger, E. 1985, ApJS, 58, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Bertschinger, E. 1998, ARA&A, 36, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J. 2004, MNRAS, 350, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Boily, C. M., Athanassoula, E., & Kroupa, P. 2002, MNRAS, 332, 971 [NASA ADS] [CrossRef] [Google Scholar]
 Burkert, A. 1990, MNRAS, 247, 152 [NASA ADS] [Google Scholar]
 Cannizzo, J. K., & Hollister, T. C. 1992, ApJ, 400, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Carron, J., & Szapudi, I. 2013, MNRAS, 432, 3161 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, B., & Henriksen, R. N. 1991, J. Math. Phys., 32, 2580 [NASA ADS] [CrossRef] [Google Scholar]
 Cheng, C. Z., & Knorr, G. 1976, J. Comput. Phys., 22, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Colombi, S. 2001, New Astron., 45, 373 [CrossRef] [Google Scholar]
 Colombi, S., & Touma, J. 2008, Commun. Nonlinear Sci. Numer. Simul., 13, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Colombi, S., Sousbie, T., Peirani, S., Plum, G., & Suto, Y. 2015, MNRAS, 450, 3724 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & McLaughlin, D. E. 2005, MNRAS, 363, 1057 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Dehnen, W., & Read, J. I. 2011, Eur. Phys. J. Plus, 126, 55 [CrossRef] [Google Scholar]
 Dolag, K., Borgani, S., Schindler, S., Diaferio, A., & Bykov, A. M. 2008, Space Sci. Rev., 134, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Fillmore, J. A., & Goldreich, P. 1984, ApJ, 281, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Fujiwara, T. 1981, PASJ, 33, 531 [NASA ADS] [Google Scholar]
 Fujiwara, T. 1983, PASJ, 35, 547 [NASA ADS] [Google Scholar]
 Gott, III, J. R., & Thuan, T. X. 1976, ApJ, 204, 649 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. 1964, Ann. d’Astrophys., 27, 83 [NASA ADS] [Google Scholar]
 Hénon, M. 1973, A&A, 24, 229 [Google Scholar]
 Henriksen, R. N., & Widrow, L. M. 1995, MNRAS, 276, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N., & Widrow, L. M. 1997, Phys. Rev. Lett., 78, 3426 [NASA ADS] [CrossRef] [Google Scholar]
 Hjorth, J., & Williams, L. L. R. 2010, ApJ, 722, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles (Bristol: Hilger) [CrossRef] [Google Scholar]
 Hozumi, S., Fujiwara, T., & KanYa, Y. 1996, PASJ, 48, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Hozumi, S., Burkert, A., & Fujiwara, T. 2000, MNRAS, 311, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., & Worrakitpoonpon, T. 2011, Phys. Rev. E, 84, 011139 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, M., Marcos, B., & Sylos Labini, F. 2009, MNRAS, 397, 775 [NASA ADS] [CrossRef] [Google Scholar]
 Kandrup, H. E., & Smith, Jr., H. 1991, ApJ, 374, 255 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 1967, MNRAS, 136, 101 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Eggleton, P. P. 1980, MNRAS, 191, 483 [NASA ADS] [Google Scholar]
 LyndenBell, D., & Wood, R. 1968, MNRAS, 138, 495 [NASA ADS] [Google Scholar]
 Makino, J. 1996, ApJ, 471, 796 [NASA ADS] [CrossRef] [Google Scholar]
 Maréchal, L., & Perez, J. 2011, Transp. Theory Stat. Phys., 40, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Melott, A. L. 2007, ArXiv eprints [arXiv:0709.0745] [Google Scholar]
 Melott, A. L., Shandarin, S. F., Splinter, R. J., & Suto, Y. 1997, ApJ, 479, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D., & Aguilar, L. A. 1985, MNRAS, 217, 787 [NASA ADS] [CrossRef] [Google Scholar]
 Mohayaee, R., & Shandarin, S. F. 2006, MNRAS, 366, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Nishida, M. T., Yoshizawa, M., Watanabe, Y., Inagaki, S., & Kato, S. 1981, PASJ, 33, 567 [NASA ADS] [Google Scholar]
 Polyachenko, V. L., & Shukhman, I. G. 1981, Sov. Astron., 25, 533 [NASA ADS] [Google Scholar]
 Polyachenko, E. V., & Shukhman, I. G. 2015, MNRAS, 451, 601 [NASA ADS] [CrossRef] [Google Scholar]
 Pontzen, A., & Governato, F. 2013, MNRAS, 430, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Romero, M., & Ascasibar, Y. 2018, MNRAS, 479, 4225 [NASA ADS] [CrossRef] [Google Scholar]
 Sikivie, P., Tkachev, I. I., & Wang, Y. 1997, Phys. Rev. D, 56, 1863 [NASA ADS] [CrossRef] [Google Scholar]
 Sousbie, T., & Colombi, S. 2016, J. Comput. Phys., 321, 644 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, J. E., & Navarro, J. F. 2001, ApJ, 563, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S., Hénon, M., & LyndenBell, D. 1986, MNRAS, 219, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Mohayaee, R., & White, S. D. M. 2001, MNRAS, 414, 3044 [NASA ADS] [CrossRef] [Google Scholar]
 Watanabe, Y., Inagaki, S., Nishida, M. T., Tanaka, Y. D., & Kato, S. 1981, PASJ, 33, 541 [NASA ADS] [Google Scholar]
 Yamaguchi, Y. Y. 2008, Phys. Rev. E, 78, 041114 [NASA ADS] [CrossRef] [Google Scholar]
 Yamashiro, T., Gouda, N., & Sakagami, M. 1992, Prog. Theor. Phys., 88, 269 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1. Snapshots of the phasespace density for the simulations with “warm” initial conditions, η ≃ 0.5. A typical slice of f(r, υ_{r}, j) with j = 0.16 is shown in (r, υ_{r}) space for the VlaSolve simulations at t = 0, at an early time t = 15, at an intermediate time used to perform tests of selfsimilarity of the phasespace spiral and at the final time. 

In the text 
Fig. 2. As in Fig. 1 but for the simulations with “cool” initial conditions, η ≃ 0.1 and for a slice with j = 0.06. 

In the text 
Fig. 3. Closer examination of the onset of instabilities in phase space for the (η, n) = (0.1, −1.5) simulations: effects of spatial and mass resolutions. For the same angular momentum slice, j = 0.06, as in lower panels of Fig. 2, the phase space density is represented in (r, υ_{r}) space. The two first lines of panels correspond to two VlaSolve simulations with respective resolutions (N_{r}, N_{v}, N_{j}) = (2048, 2048, 128) and (1024, 512, 512). The next two lines of panels correspond to two Gadget2 simulations with respective numbers of particles N = 10^{7} and N = 10^{8} and the last line of panels gives, for N = 10^{7}, the result obtained for the shell code. For the Nbody simulations, the phasespace density is sampled on grids with resolution (N_{r}, N_{v}, N_{j}) = (1024, 1024, 32). One can notice that the phasespace sheet is less well defined or “fuzzier” in the Nbody simulations than in the Vlasov code at low radius; this is because what is actually plotted is the distribution of particles (or shells) in a relatively large interval of angular momentum j ∈ [0.056, 0.077] to have sufficient number of particles to trace the phasespace distribution function, while for the Vlasov simulation, we selected only the slice corresponding to the value of j of interest. This figure illustrates the effect of radial instabilities and their dependence on spatial resolution (for the Vlasov code) or mass resolution (for the Nbody code). We note the close agreement between the shell simulation and the Gadget2 runs with 10 million particles in the middle panels of the third and fifth lines, while they diverge in the right panels when radial orbit instability effects become prominent in Gadget2. 

In the text 
Fig. 4. Radial density profile measured at various times in the VlaSolve simulations, namely initial conditions (blue dashes), intermediate time used to perform tests of selfsimilarity of the phasespace spiral (green), and final time (red dashes). For the final time, the results are also compared to a Gadget2 run (thick grey), as well as the output of the shell code (thick grey dashes) for (η, n) = (0.1, −1) and (0.1, −1.5). In addition, the logarithmic slopes −4 and −2.1 (as measured in Hozumi et al. 2000) are shown respectively as a thin solid and a thin dashed line. 

In the text 
Fig. 5. Velocity anisotropy and deviation from sphericity. Left panel: Velocity anisotropy parameter as a function of time for the VlaSolve (solid lines) and Gadget2 simulations we performed (dashed curves). Right panel: Evolution of the departure from spherical symmetry for the two kind of initial conditions experiencing radial orbit instability in Gadget2, namely (η, n) = (0.1, −1.0) (blue curves) and (η, n) = (0.1, −1.5) (red curves, solid and dashed for the 10 and 100 million particles simulations, respectively). The quantity indicates the deviation from unity of r_{axis} (in log scale and with a negative sign for r_{axis} > 1), where r_{axis} = b/c (upper curves) or b/a (lower curves), and a ≤ b ≤ c are the principal axis lengths of the Gadget2 particle distribution derived from the inertia tensor. 

In the text 
Fig. 6. Small and large radii regimes in the (η, n) = (0.5, −1) VlaSolve run at t = 40. Left panel: plots separately the magnitude of centrifugal force j^{2}/r^{3} and of the gravitational force as functions of radius as well as the magnitude of the sum of both forces. Right panel: isocontours of the specific energy (white curves) superposed on the phasespace distribution function for j = 0.16. In addition, the white dashes and dotdashes assume that only the gravitational or the centrifugal force contributes, respectively. 

In the text 
Fig. 7. Illustration of the method used to determine the positions of the folds and corresponding interfold distance law at small radius in a phasespace slice. Top panel: zoom performed around the axis υ_{r} = 0 in the region dominated by centrifugal force for the (η, n) = (0.5, −1.0) simulation at t = 40. The corresponding phasespace distribution function f(r, υ_{r} = 0, j = 0.16) is plotted in the lower panel. The black dots give the positions of local maxima estimated with our local quadratic fit, while the blue and red dots provide upper and lower bounds to compute the (very conservative) error bars shown on Fig. 8. 

In the text 
Fig. 8. Interfold distance at null radial velocity vs. selfsimilar predictions for a fixed value of angular momentum. The distance dr between local maxima of the function f(r, υ_{r} = 0, j) is plotted as a function of r for the VlaSolve runs with various initial conditions. The time considered corresponds to the middle column of panels in Figs. 1 and 2. The two left columns of panels correspond to the “warm” case, η = 0.5 with j = 0.16, and the two right ones to the “cool” case, η = 0.1 with j = 0.06. Then, odd column numbers (1 and 3) and even column numbers (2 and 4) correspond to the regime where centrifugal/gravitational force dominates, respectively. On each panel a red line indicates the logarithmic slope predicted by selfsimilarity. When the gravitational force dominates, two additional dashed curves provide a bracket of the estimated slope taking into account deviations from selfsimilarity, i.e. variations of the effective logarithmic slope of the gravitational force. 

In the text 
Fig. 9. Spiral shape vs. selfsimilar predictions for the VlaSolve run with (η, n) = (0.5, −1.0). Left panel: comparison, at t = 40, of the local shape of the spiral structure in the region dominated by centrifugal force to the curve given by Eq. (46) (red curve). Middle panel: test for selfsimilarity in time. The green spiral shape obtained at t = 30 is rescaled according to Eqs. (47) and (48) to be compared to the blue one, in the regime dominated by centrifugal force. Right panel: using only the determination of the position of the folds in the region dominated by centrifugal force, it is possible to draw the full shape of the spiral if the gravitational potential is known by interpolating the specific energy in angle coordinate defined by Eqs. (49)–(51). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.