Long-lived transient structure in collisionless self-gravitating systems

TL;DR

This paper demonstrates through simulations that purely self-gravitating systems can form long-lived, complex, non-stationary structures like spirals and rings, arising from simple initial conditions and out-of-equilibrium dynamics.

Contribution

It reveals that collisionless gravitational relaxation can produce persistent, intricate galactic-like structures without dissipative processes, challenging traditional views.

Findings

Structures can last longer than dynamical times but shorter than relaxation times.

Simulations show formation of spiral arms, bars, and rings in non-spherical, rotating systems.

Outer regions exhibit radial motions while cores are virialized and rotating coherently.

Abstract

The evolution of self-gravitating systems, and long-range interacting systems more generally, from initial configurations far from dynamical equilibrium is often described as a simple two phase process: a first phase of violent relaxation bringing it to a quasi-stationary state in a few dynamical times, followed by a slow adiabatic evolution driven by collisional processes. In this context the complex spatial structure evident, e.g., in spiral galaxies is understood either in terms of instabilities of quasi-stationary states, or a result of dissipative non-gravitational interactions. We illustrate here, using numerical simulations, that purely self-gravitating systems evolving from quite simple initial configurations can in fact give rise easily to structures of this kind of which the lifetime can be large compared to the dynamical characteristic time, but short compared to the collisional relaxation time scale. More specifically, for a broad range of non-spherical and non-uniform rotating initial conditions, gravitational relaxation gives rise quite generically to long-lived non-stationary structures of a rich variety, characterized by spiral-like arms, bars and even ring-like structures in special cases. These structures are a feature of the intrinsically out-of-equilibrium nature of the system's collapse, associated with a part of the system's mass while the bulk is well virialized. They are characterized by predominantly radial motions in their outermost parts, but also incorporate an extended flattened region which rotates coherently about a well virialized core of triaxial shape with an approximately isotropic velocity dispersion. We discuss the possible relevance of these simple toy models to the observed structure of real galaxies emphasizing the difference between dissipative and dissipationless disc formation.

Tables

Table 1. Table 1: Values of the relevant parameters (see text) for the specific simulations we report here.

Name	$\frac{a_1}{a_3}$	$\frac{a_2}{a_3}$	$\frac{R}{a_3}$	$\left|b_{rot}\right|$	$\frac{M_e}{M_s}$	$N_c$	$\frac{{\ell }_c}{{\Lambda }_c}$	$\lambda$
A1	2	1	–	0.2	–	–	–	0.19
A2	2	1.25	–	1.0	–	–	–	0.30
B1	1.5	1	0.1	0.25	1	–		0.12
B1a	1.5	1	0.01	0.25	1	–		0.29
B2	1.5	1.5	0.1	0.50	3	–		0.29
C1	–	–	–	1	–	5	0.5	0.28
C2	–	–	–	0.8	–	10	0.5	0.17

Equations

$$b_{rot}=\frac{2K_{rot}}{W_0},$$

$$\lambda =\frac{J|E{|}^{1/2}}{G{M}^{5/2}},$$

$$R^g(t)=\frac{G{M}_b(t{)}^2}{|W_b(t)|},$$

$${\lambda }_i=\frac{1}{5}{M}_b(a_j^2+a_k^2)$$

$$\iota =\frac{{\lambda }_3}{{\lambda }_1}-1,$$

$$\tau =\frac{{\lambda }_3-{\lambda }_2}{{\lambda }_3-{\lambda }_1}$$

$$\varphi =\frac{{\lambda }_3-{\lambda }_1}{{\lambda }_2+{\lambda }_1}.$$

$$v_r=\frac{\vec{v}\cdot \vec{r}}{|\vec{r}|}$$

$${\vec{v}}_t(r)=\frac{\vec{r}\times \vec{v}(r)}{|\vec{r}|},$$

$$\beta (r)=1-\frac{⟨v_t^2⟩}{2⟨v_r^2⟩},$$

$$a_r={\dot{v}}_r-\frac{v_t^2}{r}={\dot{v}}_r-a_c,$$

$$\zeta =\frac{⟨{\dot{v}}_r⟩}{⟨a_r⟩}.$$

$${\tau }_d=\sqrt{\frac{{\pi }^2}{8GM}},$$

$$n_{rev}=\frac{T}{2\pi R/v}\approx 30\frac{T_{10}{v}_{200}}{R_{10}}$$

Full text

I Introduction

The dynamical evolution of many particles solely interacting by Newtonian gravity is a fundamental paradigmatic problem in physics, which is essential for the modeling and interpretation of astrophysical structures. The study of this problem can also be placed in the broader framework of long-range interactions, which, from the point of view of statistical mechanics, share essential features (for a review, see, e.g., Campa et al. (2014); Chavanis (2013); Levin et al. (2014)). An approach based on equilibrium statistical mechanics (for the case of gravity, see, e.g., Padmanabhan (1990); Chavanis et al. (2002)) is physically relevant to such systems only on time scales that are very long compared to the those characteristic of the mean-field dynamics driven by the collective force fields sourced by many particles, and described by Vlasov equations. This dynamics leads typically to dynamical equilibria referred to variably as “virial equilibria”, or “collisionless equilibria” or “quasi stationary states” (QSS) — we will adopt the latter term here. These states are interpreted as stationary states of the appropriate Vlasov equation. Such states can also manifest instabilities leading to further evolution, giving rise, for example, to changes in symmetry (through “radial orbit instability”) and/or to the development of complex spatial structure (through, e.g., spiral wave instabilities). On longer time scales, diverging in particle number when expressed in terms of the “dynamical time” characteristic of the mean-field time, and described theoretically by a broader framework than the Vlasov equations incorporating “collisional” terms, such systems then typically evolve adiabatically through QSS, finally evolving to thermal equilibrium if such a state is well defined (see, e.g., Dauxois et al. (2002); Yamaguchi et al. (2004); Campa et al. (2009); Joyce and Worrakitpoonpon (2010); Marcos (2013); Benetti et al. (2014); Campa et al. (2014). )

For the case of gravity in three dimensions, numerical study indicate (see, e.g., Marcos et al. (2017) and references therein) that collisional evolution is driven primarily by two-body collisions on a time scale of the order of $\tau_{coll}\sim(N/\ln(N))\tau_{dyn}$ where $\tau_{dyn}\sim 1/\sqrt{G\rho}$ (where $\rho$ is the system’s average density) as originally argued by Chandrasekhar (1943). In most astrophysical systems the timescale for such relaxation is this much longer than the Hubble time and their dynamics, on relevant timescales, is thus expected to be accurately described by the collisionless (mean-field) dynamics. The framework for the study of stellar and galactic dynamics (see, e.g., Binney and Tremaine (2008)) is thus that of the collisionless Boltzmann equation (i.e,. the Vlasov equation coupled to the Poisson equation). Because of the extremely long time scale of collisional relaxation, the assumption of stationarity is often used as a working hypothesis for studying the structure of galaxies, with possible secular evolution through collisionless dynamics. In this context the striking spatial structuration of galaxies — most evidently spiral structure — has been described analytically in terms of instabilities of such QSS (see, e.g., Dobbs and Baba (2014) and references therein). With modern numerical studies of galaxy formation (see, e.g., Gott (1977); Lake and Carlberg (1988); Katz (1991); Katz and Gunn (1991); Katz (1992); Steinmetz and Mueller (1994); Steinmetz and Muller (1995); Marinacci et al. (2014); Naab and Ostriker (2017) and references therein) there has been extensive study of the origin of such structure, but an essential role in its generation is then played by dissipative non-gravitational processes. In this paper we report results showing that structures of this kind can arise purely within the restricted framework of self-gravitating systems. Further these structures, despite the fact that they are generated by a collisionless dynamics, are the result of a far from equilibrium state which persists for times which are very long compared to the dynamical time on which the initial approximate virialization of the system occurs.

To study these transient configurations we evolve numerically, using gravitational N-body simulations, a set of toy models; this allows us to understand the basic physical processes which give rise to such out-of equilibrium structures. The relaxation of isolated self-gravitating systems from simple initial conditions (IC), of the kind we consider here, has been extensively studied in numerical simulations over several decades (see, e.g., Henon (1973); van Albada (1982); Aarseth et al. (1988); David and Theuns (1989); Aguilar and Merritt (1990); Theuns and David (1990); Boily et al. (2002); Barnes et al. (2009)). The focus of such studies has, almost exclusively, been on the formation and properties of the virialized structures which are very efficiently produced by the collapse’s dynamics. In particular, in the context of the theory of galaxy formation, there was much interest in the capacity of such IC to produce structures resembling elliptical galaxies. The focus of our study here is, instead, on an aspect of these systems which has been overlooked: the production of rich transient structures from the small, but non-negligible, fraction of loosely bound and ejected mass which is very generically also produced by the relaxation process (Joyce et al., 2009; Sylos Labini, 2012, 2013). This phenomenon is, we believe, of basic theoretical interest and potentially of considerable relevance to understanding astrophysical structures. In a recent article (Benhaiem et al., 2017), we have shown that, starting from a specific class of simple idealized IC — uniform rotating ellipsoidal configurations — the relaxation of the system under its self-gravity typically leads to extended transient structures resembling qualitatively that of the outer parts of spiral galaxies. In the study reported in this paper, we investigate a much broader range of IC, and also with a greater range of particle number, whether these phenomena occur more generically. Our principle finding is that the generation of such structures — which, although transient in nature, may be very long-lived (in units of the system’s characteristic dynamical time) — is indeed a quite robust and generic feature of violently relaxing systems. Further the morphologies of the structures produced are even more rich and diverse than we had anticipated.

The spiral-like structures generated far out-of-equilibrium in the systems we study arise in a manner very different to that usually envisaged in the modeling of such structures in the astrophysical literature. The mechanisms proposed are perturbative in nature, envisaging the spiral-like structure as the result of instability of an equilibrium disc configuration (see, e.g., Dobbs and Baba (2014); Binney and Tremaine (2008)). In addition, it is well known that the formation of flat disc-like configurations can be obtained with simulations of collapsing spheroids where, in addition to gravity, the dissipative effects of several astrophysical processes are taken into account (see, e.g., Gott (1977); Lake and Carlberg (1988); Katz (1991); Katz and Gunn (1991); Katz (1992); Steinmetz and Mueller (1994); Steinmetz and Muller (1995) and references therein). Indeed, in this context, an indispensable element for the appearance of such structures is the inclusion of gas with a dissipative dynamics (e.g., cooling). The principle motivation for including this dynamics was initially to allow such structures to be produced. The central finding of our study, in contrast, is that disc-like configurations with transient spiral arms and with bars and/or rings may be produced by a purely ( i.e., non dissipative) gravitational dynamics. We will discuss below in further detail this essential difference with respect to previous works in the literature. We also stress the peculiar features of the complex velocity fields generated by the gravitational dynamics we have investigated, which are different from the structures of this kind generated by dissipative dynamics. Indeed, while in the latter case the velocity fields are essentially dominated by rotational motions, in the former case these show different regimes and are essentially radial in the systems’ outermost regions.

The article is organized as follows. In Sect.II we describe the implementation of the numerical simulations, discussing our choice of numerical parameters, of IC, and the specific quantities we have measured. In Sect.III we present our results, focusing in particular on the spiral-like, bars and rings transients produced. We describe in considerable detail the mechanism producing them, which varies in detail from one kind of IC to the other. We turn in Sect.IV to a brief discussion of the possible relevance to real galactic structures of our simple models. The differences between dissipative and non-dissipative dynamics in the process of the formation of a disc-like structure is discussed in Sect.V. Finally in Sect.VI we summarize our findings and outline some of the many further questions raised by them.

II Numerical Simulations

Our numerical experiments consist of molecular dynamics simulations of an Hamiltonian system of $N$ point particles in an infinite (three dimensional) domain, interacting by a pair potential $\phi$ which is that of Newtonian gravity with a short distance regularization. They have been performed using the publicly available (and widely used) code Gadget-2 (Springel, 2005), in the appropriate version (non-expanding universe, open boundary conditions). The regularization of the potential — the so-called “gravitational softening” — in this code is implemented by solving the Poisson equation for spherical continuous clouds, with compact support of characteristic radius $\varepsilon$ (and total mass equal to that of the particle). The force thus takes exactly its Newtonian value at separations greater than $\varepsilon$. The functional form of the regularized potential for $r<\varepsilon$, of which the exact expression can be found in (Springel, 2005), is a cubic spline interpolating between the exact Newtonian potential at $r=\varepsilon$ and a constant value at $r=0$, with vanishing gravitational force at this point.

II.1 Precision and softening

For simulations of the dynamics from the IC we consider, which give rise typically to very significant contraction of the system — leading to very high densities and short characteristic time-scales — the choice of the force softening length and numerical parameters controlling the accuracy of the time-stepping and force calculation requires particular care.

We have performed simulations, which we call “low precision” (LP), with the numerical parameters of the code set at the values suggested by the Gadget-2 user guide 111See https://wwwmpa.mpa-garching.mpg.de/gadget/users-guide.pdf.. We monitor the total energy and find that these runs typically conserve it to within about $0.5\%$. We have then also run “high precision” (HP) simulations using values of the relevant parameters which lead to energy conservation to within one part in $10^{5}$ 222Specifically the parameter ErrTolTheta is $10^{-10}$ for HP compared to $0.7$ for LP, while ErrTolForceAcc is $10^{-10}$ for HP and $5\times 10^{-4}$ for LP.. We also monitor conservation of total linear and angular momentum conservation and observe similar results.

We have run both LP (with $N$ in the range $[10^{4},5\times 10^{5}]$) and HP (with $N$ in the range $10^{6}$) simulations of realizations of our IC (described in detail below), and found no apparent significant differences for the macroscopic quantities which we study. In what follows we will report only results for LP simulations, which likewise give results completely consistent with the lower $N$ simulations at HP and LP. Thus it appears that our results are $N$ independent — and represent those of an $N\rightarrow\infty$ limit. However, as we discuss further below, such a conclusion should be treated with caution.

Our default value of the force softening length (the code’s parameter $\varepsilon$), and that used in the specific simulations reported below, is $10^{-3}$ in the units of length we define further below. In practice this means it is $10^{-3}$ of the smallest characteristic length in the IC (e.g., shortest semi-principal axis in the case of an initial ellipsoid), and is such that it is always considerably smaller than the typical size of the system when it is most contracted. We have also performed extensive tests to control the effect of varying the force-softening parameter, and found that we indeed obtain very stable results, provided this condition is respected on the comparison of $\varepsilon$ and the minimal size reached by the whole structure during the collapse.

We have also run some test simulations in which the particles have two different masses. Specifically in a reference simulation in which all the particles have mass $m$, we have re-sampled randomly with particles of mass $\alpha\times m$ and $m/\alpha$, with $\alpha=2,5,10$, determining their number so that the total mass is fixed. We have then studied the spatial and velocity properties of the two species separately, and found them to be in good agreement both with one another and with the properties found in the original single mass simulation of the IC. This provides strong evidence that the dynamics producing the distributions we have described are indeed representative of the mean-field (or collisionless) dynamics of the continuum IC we have sampled. Indeed this is consistent with what one would expect as, even for the longest times we simulate (at very most $200$ dynamical times), collisional two body effects would be expected to be negligible for the particle numbers we consider even in the virialized core of the system.

II.2 Initial conditions

We wish to study evolution of self-gravitating $N$ body systems starting from IC which are sampled from spatial mass distributions which break rotational symmetry, and also mass distributions which are non-uniform (i.e. for which the mass density is not constant in the compact region where it is non-zero). Clearly the space of parameters that describe a generic IC of this type is infinite. We have chosen three very simple few parameter families of such IC. In order of increasing complexity, we consider mass distributions which are (i) uniform ellipsoids, (ii) uniform ellipsoids with a central spherical region of higher density, and, (iii) a collection of uniform spheres of varying radii with centers randomly placed in a spherical region. For the initial velocities we consider only two very simple cases: coherent rigid body type rotation of the whole structure, and in a few cases only, some additional random uncorrelated motion. The precise details of our choices for how the different IC are then parametrized, are given below.

With these choices of IC we aim to single out the effect of the initial shape of an isolated mass distribution on its collapse and subsequent evolution. This aspect appears to have been largely overlooked in the literature, which has focused instead mostly on the effect of internal correlations between density fluctuations, for the case of a spherical cloud (Gott, 1977; Lake and Carlberg, 1988; Katz, 1991; Katz and Gunn, 1991; Katz, 1992; Steinmetz and Mueller, 1994; Steinmetz and Muller, 1995). As we discuss below, one of our main results is that the spherical IC represent a very peculiar and pathological case that give rise a to dynamics that suppresses a common characteristic feature of the gravitational collapse of isolated over-densities: the formation of a disc-like structure with spiral type arms that occurs when the IC significantly break spherical symmetry.

II.2.1 Uniform ellipsoidal clouds

For this class of systems the particles are

• •

randomly distributed, without correlation, with uniform probability inside an ellipsoid, and

• •

given a coherent rigid body rotation about the shortest semi-principal axis.

The three parameters which we choose to characterize them are:

• •

The ratios, $a_{1}/a_{3}$ and $a_{2}/a_{3}$, of the two longer semi-principal axes ($a_{1}\geq a_{2}$) to the shortest one $a_{3}$.

• •

The virial ratio associated with the rotation

[TABLE]

where $K_{rot}$ is the kinetic energy of the rotational motion, and $W_{0}$ is the initial gravitational potential energy. As mentioned above, we have considered some numerical experiments in which we add a random motion in addition to the coherent rotational one. The amplitude of this motion has been taken to be a fraction of the order of $10-20\%$ of the rotational one. We have not found any significant difference on a macroscopic scale, but the transient features like the spiral arms are indeed more diffuse when random motions have the highest amplitude we have tested.

We have performed simulations for a large variety of prolate (i.e., $a_{1}>a_{2}=a_{3}$), oblate (i.e., $a_{1}=a_{2}>a_{3}$) and triaxial (i.e., $a_{1}>a_{2}>a_{3}$) shapes, as well as the special case when the cloud is spherical (and thus rotational symmetry is broken only by the finite $N$ Poissonian fluctuations). We have explored the range of values defined by the ratio $a_{1}/a_{3}$, $a_{2}/a_{3}$ $\in[0.05,2]$. For the velocities we have explored the range $b_{rot}\in[-1,-0.05]$. We report in detail results for just two representative cases, A1 and A2, for which the parameter values are given in the Table 1. We report also the value of the so-called “spin parameter” as defined in Peebles (1969, 1971):

[TABLE]

where $J$ is the total angular momentum of the system, $E$ the binding energy, $M$ the total mass and $G$ Newton’s constant. This parameter is widely used in the astrophysical context when characterizing the angular momentum of astrophysical systems (see also, e.g., Barnes and Efstathiou (1987); Bullock et al. (2001)). It simply provides a dimensionless measure of the angular momentum in the natural units of a self-gravitating system (given by dimensional analysis by the combination of $GM^{5/2}/E^{1/2}$). It thus provides an indication of the degree of rotational support of a self gravitating system provided by its angular momentum.

II.2.2 Non-uniform ellipsoidal clouds

In this class of systems we consider a coherently rotating ellipsoidal configuration exactly as described above, but then modify it in a spherical region around its center, ascribing

• •

a different number (and thus) mass density of the uniformly distributed particles, and

• •

only random uncorrelated velocities, sampled uniformly in velocity space up to a maximal magnitude.

We consider only the case that the kinetic energy $K_{s}$ of the particles in the spherical region is such that $2K_{s}=-W_{s}$ where $W_{s}$ is their potential energy, i.e., this central region is initially approximately virialized. The initial conditions are thus chosen to probe the dynamics of the collapse of a rotating ellipsoidal cloud which already has a smaller virialized structure inside it.

There are then five parameters characterizing this family of IC, which we choose to be: $a_{1}/a_{3}$, $a_{2}/a_{3}$, $R/a_{3}$ (where $a_{i}$ are again the lengths of the semi-principal axes, and $R$ the radius of the sphere), the ratio $M_{e}/M_{s}$ of mass (i.e. particle densities) in the ellipsoidal region to that in the sphere, and the ratio $b_{rot}$ as given in (1), with $W_{0}$ the total potential energy of the configuration of the ellipsoid alone.

We have considered simulations in which $b_{rot}\in[-1,-0.1]$ which means that in almost all cases all particles are initially bound. For the other parameters we have explored $a_{1}/a_{3}$, $a_{2}/a_{3}$ $\in[0.05,2]$, $R/a_{3}\in[0.01,0.1]$, and $M_{e}/M_{s}\in[1/3,3]$.

We note that the characteristic time for collapse of the rotating cloud, in units of the dynamical time of the sphere, is $\sim(a_{3}/R)^{3/2}\dot{(}M_{s}/M_{e})^{1/2}$. Thus we explore the range in which this ratio is between about ten to thirty. We report in detail results for just three representative cases, B1, B1a and B2, for which the parameter values are given in Table 1.

II.2.3 Non-spherical and non-uniform clouds

The third class of models mass distributions which still rotate coherently, but in which rotational symmetry is broken in a more random and less idealized manner.

Specifically:

• •

We choose $N_{c}$ points randomly with uniform probability in a sphere of radius $R_{0}$.

• •

Taking each of these $N_{c}$ points as centers, we distribute randomly $N_{p}$ points in spheres of radius $\ell_{c}$ centered on them.

• •

We calculate the moment of inertia tensor and use it to determine the direction of the principal axis with the largest eigenvalue. We give a coherent rigid body rotation to the whole cloud about this axis.

For a given total particle number, this is thus a three parameter family of configuration. We take these parameters to be $N_{c}$, $\ell_{c}//\Lambda_{c}$ and $b_{rot}$, where the latter is given again by equation (1) with $W_{0}$ the total initial potential energy. $\Lambda_{c}$ here is the mean distance between neighboring clouds, $\Lambda_{c}\approx 0.55(4\pi R_{0}^{3}/3N_{c})^{1/3}\;$, and the parameter $\ell_{c}//\Lambda_{c}$ thus characterizes the degree of overlap of the individual clouds. We are in practice interested in values not so different from unity so that the initial density fluctuations are not so large, and there is a global collapse of the whole structure (rather than separate collapses for the sub-structures and the whole structure). We have studied the range of parameters $N_{c}\in[3,50]$, $\eta=l_{c}/\Lambda_{c}\in[0.1,2]$and $b_{rot}\in[-1,-0.1]$.

We again report results for just two cases, C1 and C2, specified in the Table 1, whose features are representative of this class of models.

It is important to note that our random sampling with a finite number $N$ of particles introduces mass density fluctuations, which are additional to those intrinsic to our continuum IC. Such density fluctuations can of course play a role in the dynamics, and as their amplitude is $N$-dependent (with $\frac{\delta\rho}{\rho}\sim 1/\sqrt{N}$), one might expect this to induce potentially an $N$-dependence in our results even for macroscopic results. However, as detailed above, the parameters of our IC are in fact chosen so that $\frac{\delta\rho}{\rho}\sim 1$, so one might expect the finite $N$ fluctuations to be negligible. We discuss in more detail finite $N$ effects and the problem of taking the continuum limit in Sect.III.11.

As noted above this is indeed consistent with our numerical findings over a range of $N$ between $10^{4}$ and $10^{6}$. For this reason we present in the paper only results for the case $N=10^{6}$ and do not discuss any further the role of $N$ in our results. Nevertheless we caution that it is quite possible that the physical processes we simulate might have subtle dependencies on $N$ which do not show up clearly in the range we have simulated. Notably, as illustrated by the study in (Benhaiem et al., 2018) of the special case of spherical IC, such dependencies may arise because the finite $N$ fluctuations break symmetries of our continuum IC.

II.3 Physical quantities measured

A useful basic quantity to monitor the global evolution of the system is the “gravitational radius” defined by

[TABLE]

where $W_{b}(t)$ is the gravitational potential energy of the bound particles and $M_{b}(t)\leq mN$ is their mass.

To characterize the system’s shape we compute the three eigenvalues $\lambda_{i}$ (with $i=1,2,3$) of its inertia tensor, which, in the case of an ellipsoid, are related to the lengths $a_{i}$ of the three semi-principal axes by

[TABLE]

where $i\neq j\neq k$ and $i,j,k,=1,2,3$. It follows that $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$. It is standard then to introduce three different combinations of the $\lambda_{i}$: the flatness parameter,

[TABLE]

the triaxiality index,

[TABLE]

and the disk parameter,

[TABLE]

These parameters allow one to distinguish not only between different type of ellipsoids (e.g., prolate, oblate and triaxial) but also between other shapes ( i.e., bars vs. disks). For instance a sphere has (0, –, 0) a disk (1,0,0.5) and a narrow cylinder (i.e., a bar) ($\gg 1$,1, $\approx 1$).

In addition to the radial component of a particle’s velocity $\vec{v}$,

[TABLE]

we define the vectorial “transverse velocity” as

[TABLE]

i.e., the vector of which the magnitude is that of the non-radial component of the velocity, but oriented parallel to the particle’s angular momentum relative to the origin.

We will denote the average of a quantity in a spherical shell about radius $r$ by $\langle\cdots\rangle$. Coherent rotation of all the particles in the shell about the same axis then corresponds to $\langle|\vec{v}_{t}|\rangle$= $|\langle\vec{v}_{t}\rangle|$. Furthermore we consider the anisotropy parameter

[TABLE]

where $\langle v_{t}^{2}\rangle$ and $\langle v_{r}^{2}\rangle$ are respectively the average square value of the transversal and radial velocity. The anisotropy parameter has the following limiting behaviors: $\beta\rightarrow 0$ for an isotropic velocity distribution, and $\beta\rightarrow 1$ when ${\langle v_{t}^{2}\rangle}\ll{2\langle v_{r}^{2}\rangle}$, i.e. when the motion is predominately radial.

To characterize the kinematics further, we also consider the different components of the radial acceleration, which can be decomposed as

[TABLE]

where $a_{c}$ is the magnitude of the centripetal acceleration associated with the transverse component of the velocity. To quantify in a simple manner the degree of circular vs. radial motion, we will consider the ratio

[TABLE]

When particles’ motion is purely circular we thus have $\zeta=0$, while if it is purely radial we have $\zeta=1$: thus $\zeta$ captures different properties of the velocity field than $\beta$, although they both tend to unity when the motion is purely radial.

III Results

The phenomenology of the gravitational collapse of the IC we study is, in many respects, very similar to that of non-rotating isolated clouds discussed at length in previous works. We first summarize these behaviors, and in particular recall the mechanism by which particles may gain enough energy even to be ejected from the system. We then subsequently focus on the features of the evolving system which are specifically associated with the initial rotational motion, and in particular the emergence of long-lived spiral-like structure, as well as transient bars and/or rings, in the spatial configuration. We recall that all the results presented explicitly in the article are for simulations with $N=10^{6}$ particles, but that all the quantities considered have shown no apparent $N$ dependence for simulations of the same IC with $N$ ranging from $10^{4}$ to $5\times 10^{5}$. As noted above, this corresponds to there being no apparent dependence on the fluctuations associated with the particle sampling of the continuous mass distributions characterizing the initial conditions.

III.1 Units

As unit of length we take $a_{3}=1$ for our first two families of IC, and $R_{0}=1$ for the third one. As unit of time we then take

[TABLE]

i.e., the characteristic time for the collapse of a sphere of radius unity (where $M$ is the total mass of the system). Finally particle energies will be given in units of ${GMm}$ where $m$ is the particle mass (and $M=Nm$).

III.2 Collapse and re-expansion

Fig. 1 (upper panel) shows the evolution of the gravitational radius (Eq.3), for three different initial conditions. These are those which show the largest and the smallest variation among the ones we have selected. In the case of A1 and C1 the system, which is initially far from equilibrium, contracts globally reaching a minimum on a time scale of order $\tau_{d}$ ($\approx 10\tau_{d}$ for the case of B1 because of the lower density of the external ellipsoid), then it re-expands and, after a number of damped oscillations which varies, rapidly settles down to a fairly stable value. A similar behavior is manifested by the virial ratio (lower panel of Fig. 1) and thus the stabilization of $R_{g}$ reflects the relaxation of the system to a state close to virial equilibrium. As we will see below, this inference is only approximately true, as a small fraction of the mass remains in a time-dependent configuration on much longer time scales: indeed, it is this fraction of the mass distribution which we will discuss at length below.

For the cases of A1 and C1 a fraction of the mass is ejected after the collapse, and as a result the global virial ratio stabilizes around a value smaller than -1; in the case of B1 the collapse is less violent and there is no mass ejected.

We note that for B1 and C1, as all other simulations of these classes, the gravitational radius is reduced by a smaller factor than for the case of homogeneous ellipsoids. Compared to this case, the fluctuations of the gravitational field generated are thus much weaker.

III.3 Particle energies distributions: before and after collapse

Fig. 2 shows the distribution $P(E)$ of particle energies $E$ at two different times in the same three simulations as in the previous figure. Plotting these distributions at longer times than the last one shown we find no noticeable evolution, i.e., these distributions represent well in principle a final stationary distribution. We see that their qualitative behavior divides them clearly as in the previous figure: in the cold simulation, A1, the change in the energy distribution brought by the dynamics is much more marked, with (i) a much more widely spread energy distribution compared to that in B1 and C1, and (ii) a significant fraction of the mass with positive energy, while there is a much smaller fraction of (or even no) such particles in the other cases. In the case of B1 there is a small but non-negligible evolution of the particle energy distribution while C1 represents an intermediate case with respect to A1 and B1.

The correlation between the behavior in Fig. 1 and Fig. 2 is simple to understand: changes in particle energy are driven by the variation of the gravitational field, and this is much more violent in the cold cases. In the warmer case the variation of the field is relatively gentle, and mixing in phase space has time to play a greater role in the relaxation process. Nevertheless, as we will see, the generation of even a small number of particles with positive energy or indeed a significant fraction of bound mass with energy close to zero, is sufficient to produce a considerable and non-trivial evolution in configuration space.

Let us recall another relevant feature of these energy changes, which we have shown in previous study of IC without rotation (Joyce et al., 2009; Sylos Labini, 2012): the particles which have a large energy gain, and which thus constitute predominantly both the unbound and loosely bound mass, are those which lie initially in the outer part of the structure. In the spherical case, these are the particles initially in the outer shells, and in the ellipsoidal case, particles which are at large radii and close to the longest semi-principal axis.

The reason for this correlation between the energy gain/loss and particle initial position, is related to the times at which particles first pass through the center of the structure: particles which pass through the center after the bulk of the particles, experience the intense gravitational field which changes their energy at a time when it is weakening, simply because the bulk of mass generating it is already re-expanding. Such particles thus fall into a potential which is deeper than the one they subsequently climb out of, and they thus have a greater net boost of their energy. In the ellipsoidal case it is quite evident why the initially outermost particles arrive late: in the evolution from such an initial condition, collapse occurs first along the shortest axis and last along the longest axis (Lin et al., 1965). In the quasi-spherical case the reason for the average late arrival of particles in the outermost shells is more subtle as it is due to a boundary effect: particles near the boundary experience a lower average density and thus have a longer fall time (see Joyce et al. (2009) for a detailed discussion).

III.4 Mean density profiles

Fig. 3 shows the mean density averaged in radial shells of equal logarithmic thickness $\Delta\left(\ln(r)\right)$ as a function of radius, for the same simulations as in the previous figures 333We take as center the particle which has the lowest gravitational potential. Note that we will use the same radial binning everywhere below when we compute averages.. We again obtain results very similar to that for non-rotating IC, with cold IC producing (i) a more compact core than the warmer IC, and (ii) a characteristic $1/r^{4}$ decay of the density at large radii. As discussed for example in Sylos Labini (2012), this latter behavior is associated with the very loosely bound particles on highly radial orbits, and can be explained in a simple manner by considering that the outermost particles move approximately in a central and stationary potential. In these plots the system appears to settle down to a stationary form on the same times scales as inferred from the plots in the previous figures (and which are slightly longer for the warm cases). We can just make out the signature of some continuing evolution at the largest radii. It is this which we now focus on.

III.5 Non-stationary features at longer times

We consider here the non-stationary features of the mass distribution at longer times, of which the very non-trivial distinctive space and velocity structure is a result of the initial rotation of the IC. In space this can be seen very evidently in the (linear scale) snapshots of the evolving spatial configurations projected on the plane orthogonal to the initial axis of rotation of the IC, shown in Fig.4. We focus first on how to quantify this non-stationary evolution which shows up the features common to all the different IC, and then discuss subsequently the genesis of the large variety of different forms which are evident.

Fig. 5 shows the particle energies averaged in shells as a function of radius, for the same three simulations as in Figs. 1-3, and for two different times. In all cases, the non-stationarity is now clearly visible in the propagation to larger radii of the outermost mass, both loosely bound and unbound. This non-stationarity is essentially a consequence of the fact that the mass in this energy range has a characteristic time for virialization with the rest of the mass which diverges as $E\rightarrow 0$ from below: to a first approximation the potential they move is central and stationary, with an associated Keplerian period $\tau\sim 1/(-E)^{3/2}$.

Fig. 5 shows clearly that the structure can be divided in an inner stationary part and an outer non-stationary part. Examination of the properties in velocity space, illustrated in Figs.6 -8, show that a further refinement into three distinct regions is warranted:

• •

For $r<R_{1}$, $\beta(r)\approx 0$, corresponding to an isotropic velocity dispersion, there is neither net radial flow nor net rotation (i.e. both $\langle v_{r}\rangle\approx 0$ and $|\langle\vec{v}_{t}\rangle|\approx 0$). This part of the distribution is the virialized core showing the approximately flat density profile discussed above.

• •

For $R_{1}\leq r\leq R_{2}$, there is no net radial flow (i.e., $\langle v_{r}\rangle\approx 0$) but there is a significant net rotation. Indeed, $|\langle\vec{v}_{t}\rangle|$ grows monotonically as a function of radius until it reaches a value where it is comparable to $\langle|\vec{v}_{t}|\rangle$, and this remains so as both quantities slowly decline over the radii up to $R_{2}$. This latter scale is defined such that $\langle v_{r}\rangle>|\langle\vec{v}_{t}\rangle|$ for $r>R_{2}$. Thus the rotational motion grows until it is close to a completely coherent one around a single axis. Correspondingly $\zeta$ is much smaller than unity as $\dot{v}_{r}$ is small compared to $a_{c}$: this region corresponds to the flattened part of the distribution where rotational motions dominate.

• •

For $r>R_{2}$, net outward radial motion dominates (i.e. $\langle v_{r}\rangle>\langle|\vec{v}_{t}|\rangle\approx 0$, $\beta\rightarrow 1$, $\zeta\rightarrow 1$ (as $a_{r}\approx\dot{v}_{r}$) and the sub-dominant transversal component of the velocity decays monotonically towards zero. We note that in B1 this region is negligible as there is almost no ejected mass.

The scale $R_{1}$ does not evolve significantly with time, corresponding to the stationarity of the region inside it. The scale $R_{2}$, which corresponds approximately to the transition from bound to unbound mass, increases monotonically with time. Indeed, the mass distribution in both the outer part of the central region, and in the entire outer region, is manifestly non-stationary on these long time scales, and remains so for arbitrarily long times.

Comparison of Fig.4 with Figs.6-8 shows a rather interesting property of this class of models: the more the shape of the structure formed after the collapse deviates from axisymmetry, the larger is the radial velocity at large distance from its center. In addition we stress that Figs.6-8 show the average components of the velocity in spherical shells, while the actual velocity fields are characterized by large anisotropies: most notably, the amplitude of the radial velocity is correlated with the semi-major axis. These features, which will be studied in greater detail in a forthcoming work, must be considered when comparing results of this class of simulations with observations (see Sect.VI).

III.6 Emergence and evolution of spiral structure

Detailed study of the temporal evolution of the configurations confirms that the mechanism for generation of the spiral-like structure evident in Fig.4 is indeed the strong injection of energy given to particles which pass through the center of the system just after the time of the maximal compression. This gives these particles a significant radial component of their motion added to the initial rotational motion. The radial distance these particles subsequently travel, once they are outside the core, because of approximate conservation of angular momentum, is correlated with the angle they move through: particles which have larger radial velocities initially thus “trail” behind particles with smaller radial velocities. For this process to generate a spiral-like structure, the only necessary additional ingredient is that the distribution of the directions of motion of these particles is anisotropic. This is indeed the case starting from these IC. As we have discussed above, as for cold non-rotating IC, the particles which gain energy are those which lie furthest from the origin initially. In the case of an ellipsoid, the radial motion arising from the energy injection is thus preferentially correlated with the longest semi-principal axis. Asymptotically, the motion of the particles with positive energy, which are furthest out, becomes purely ballistic and radial. As a result the spiral-like structure is “frozen” and stretches more and more. Nevertheless this intrinsic non-axisymmetry of the disc is typically not yet so marked even at the quite long times we show in our plots.

We consider now in more detail the different forms arising from the range of IC that we have selected. We recall that the first two simulations, A1 and A2, in Fig. 4 show the evolution of two prolate ellipsoids (see Table 1). In both cases the collapse is quite violent, with the system undergoing a very strong contraction. This leads to the production of a significant fraction (about 10 %) of positive energy particles. The collapse of most of the mass occurs along the direction of the initially shortest semi-principal axis, while the subsequent marked expansion of the system is along the direction of the longest semi-principal axis.

The simulations B1 shows a very similar morphology, but the spiral-like arms are markedly more “wound up” than in A1 and A2. Nevertheless the basic process leading to the formation of this structure is essentially the same, with an anisotropic energy injection correlated with the axis along which particles arrive latest. This difference in the winding in particular is a direct result of a much less violent collapse, in which only a relatively small fraction of the total mass participates. As a result the radial velocities injected into particles along the longest semi-principal axis are smaller, and the particles thus travel a smaller radial distance as they rotate with their initial rotational motion. The outer shells initially collapse toward the center in a contraction which is fastest along the shortest semi-axis. This contraction enhances the initial anisotropy of the system which also leads to a transient bar structure which we discuss in greater detail below.

In B1a, in which the core is a factor ten less extended than in B1, the evolution is very similar, except for the morphological details of the bar and spiral arms. This test shows that the precise morphology of the system formed after the collapse depends sensitively on the features of the IC.

The simulation B2 likewise is characterized by a very anisotropic collapse, but as it is an oblate ellipsoid, the contraction occurs along the direction of the shortest semi-principal axis, and approximate symmetry about this axis is maintained. The contraction does not change the particle energy distribution very significantly, but is enough to give to some particles a radial velocity component directed outwards. In this case the particles which “arrive late” during the collapse are those initially close to the outer shells of this plane of symmetry, and as they re-expand outward they give rise to a quite different spiral-like structure with “flocculent” multiple arms. These appear to be seeded by the growth of the density fluctuations in the system during its collapse phase: for this reason we expect that these features depend on the amplitude of the initial density fluctuations (see discussion in Sect.III.11). In addition, a transient ring structure emerges, which we discuss further below.

In B1, B1a and B2 there is, because of the much more gentle collapse, very little or no ejected mass (for B1, see Fig. 7) and the motions are only predominantly radial at longer times, and only at the very largest distances. Likewise the region where the coherent rotational motion is dominant with respect to the radial motion is much more extended than for the case of, e.g., the simulation A1. Differently to the other cases, we observe in B2 that there appear to be two distinct phases in which we see quite different spiral arms emerge.

The simulations C1 and C2 show the typical behavior we observe in this class of more inhomogeneous and anisotropic IC. Transient structures similar to those in the other cases again emerge, and the same basic physical mechanism is at play. The spiral arm structure which forms after the collapse is clearly less axisymmetric, reflecting the lower symmetry of the IC. Both the visual appearance and the structure in velocity space (see Fig. 8) show that the resultant structures are more similar to A1/A2 than B1/B2. This reflects the fact that the collapse is indeed stronger in these cases.

III.7 Formation of bar structures

We have noted that the structures we observe are generically non-axisymmetric (except in specific cases like B2 where the axisymmetry of the IC survives the collapse phase better). Thus the configuration which results is not just very flattened along the axis parallel to the initially shortest semi-principal axis, but a coherent anisotropic structure resembling a bar emerges in the plane defined by the two longer semi-principal axes in the IC. This anisotropy is, as one would anticipate, present not only in configuration space but also in the velocity field.

Fig.9, showing several snapshots of B1 projected in the $x-y$ plane, illustrates how these transient features form and grow. The particle distribution has been coarse-grained onto a $32^{3}$ grid, and the average velocity determined in each cell. We see that, up $t\approx 15$ when the global collapse of the system occurs, the velocities are essentially rotational, but progressively develop an inward radial component, simply because of the collapse dynamics. During the phase of maximal contraction, the particles originally furthest out, i.e. along the initial semi-principal axis, gain a radial velocity component directed outwards, but their total velocity remains predominantly rotational. The particles closer in initially, on the other hand, have lower energy and remain more bound around the central structure. Thus the two arms which emerge clearly are formed from two groups of particles which initially were located at opposite ends of the longest principal semi-axis. For longer time scales, i.e. $t>50$, the spiral arms and the bar structure start to be washed out by the radial component of the motion.

III.8 Formation of ring structures

We have noted the formation of a time dependent ring-like structure in simulation B2 at $t\approx 10$, in the plane corresponding to the two largest initial principal semi-axes. As can be seen in Fig.10, this is indeed a local density enhancement which expands outwards in time. Investigation confirms that it is generated by a fraction of particles moving outward at higher than average radial velocity. These are particles which were initially in or near the outermost radial shells in the plane of the oblate IC, and which received a strong energy injection from the time dependent potential generated by the collapse along the shortest axis. As these particle carry also the initial velocities of the rotation about this axis, the ring also rotates coherently. It persists up to the end of our simulation at $t=50$. As noted, we have varied the ratios $a_{1}/a_{3}$ and $a_{2}/a_{3}$, for the more general case of a triaxial ellipsoid, in a relatively wide range. We have found that, whenever we are close to an oblate ellipsoid, such a ring-like structure is formed.

III.9 Shape parameters

As we have discussed the particles which gain most energy are those which are initially furthest away from the center. For IC like the ones we consider here, this leads to a very anisotropic distribution for the loosely bound and ejected mass. The details of this anisotropic distribution depend on the IC. We recall that for cold prolate ellipsoidal IC without rotation these particles are focused into two broad “jets” in opposite directions around the initially longest semi-principal axis (Benhaiem and Sylos Labini, 2015). This is the case simply because the particles farthest from the center are indeed close to this axis initially. With additional coherent rotation, as considered here, this jet-like structure, is, as we have seen, transformed into a spiral-like structure as the particles propagate outward, and the axis defined by the farthest out particles thus remains closely correlated with the initially longest semi-principal axis.

In all our simulations the outer part of the structure is, correspondingly, very flattened in the plane orthogonal to the initially shortest semi-principal axis. Fig. 11 show the evolution of the three shape parameters $\iota,\tau,\phi$ in simulation A1, separately for the ’internal’ particles constituting the core of the structure (upper panels) and for the ’external’ particles constituting the spiral-like arms (bottom panels), where this division is defined by the radius at which the measured averaged radial density $n(r)$ discussed above reaches half its value in the core. In this case the core is a triaxial ellipsoid with, respectively, $\iota\approx 0.3\;,\phi\approx 0.1\;,\tau\approx 0.4$ and $\iota\approx 0.7\;,\phi\approx 0.2\;,\tau\approx 0.5$. On the other hand, the external particles are much more flattened with $\iota\approx 3\;,\phi\approx 0.5\;,\tau\approx 0.5$ in both cases.

A similar behavior is shown by the evolution of the flatness parameter for the simulations B1 and B2: in these cases $\iota$ is computed by considering only the external, low energy, bound particles. In simulation B1, after the contraction phase during which $\iota$ is large and fluctuating, it becomes of order one for $t>20$. In simulation B2, on the other hand, the flatness parameter remains of order one at all times without any significant variation. Simulations C1 and of C2 are in this respect similar to A1/A2.

III.10 Sub-structures

When the initial ellipsoidal deformation is large enough (e.g., $a_{1}/a_{3}\approx 2$) the collapsing system may fragment into two or more clumps which eventually merge long after the collapse : this is the case of A2 (see Fig.4). Clearly when the IC is made up of clumps with small amplitude fluctuations, as in C1 and C2, this effect is enhanced. On the other hand, during the collapse, initial density fluctuations may evolve due to gravitational instability forming small aggregates which, after the collapse, may correspond to sub-structures. These sub-structures, which are typically not virialized as they are subject to strong tidal fields, lie in the same plane defined by the jets. Eventually they will fall into the largest virialized object and be destroyed by the interaction with the core. The formation of sub-structures, anisotropically distributed on a planar configuration, around the main virialized object appears to be a generic result of the evolution from cold IC of this type (Benhaiem and Sylos Labini, 2015, 2017).

III.11 Role of density fluctuations and the continuum limit

Let us discuss further the role of density fluctuations in the collapsing cloud in an appropriate continuum limit defined by taking $N\rightarrow\infty$. How this limit is taken must be specified, as it is not unique. Indeed here there are (at least) two evident ways of taking such a limit, and the role of density fluctuations is different in each case.

First if we consider finite $N$ configurations as Poisson samplings of continuum configurations with fixed mass density $\rho_{0}$ i.e., without any intrinsic density fluctuations , we can take the limit $N\rightarrow\infty$ together with the particle mass $m\rightarrow 0$, so that $\rho_{0}=N\times m=const.$,. In this case the density fluctuations, that are proportional in amplitude to $\delta\rho/\rho_{0}\propto N^{-1/2}$, also vanish. In this case all substructures generated by the growth of density fluctuations must disappear in the continuum limit. They can in this sense then be interpreted as finite $N$ effects.

On the other hand, we can also take the continuum limit in a different way, by taking $N\rightarrow\infty$ and $m\rightarrow 0$ with $\rho_{0}=N\times m=const.$ and with $\delta\rho/\rho_{0}=const.$ (or, more precisely, keeping the statistical properties of the latter fixed): in this case the internal fluctuations of the cloud grow in the same way independently on $N$, and thus they give rise to the same substructures. Indeed, as it was shown in detail by Joyce et al. (2009) for the case of the spherical collapse model, while the whole system collapses small density fluctuations inside the cloud grow and form substructures of increasing size with time. The collapse is halted when the size of non linear structures formed inside the collapsing cloud becomes of the same order of magnitude of the cloud itself (that meanwhile is collapsing). Thus in this process there are two competing effects: the global monolithic collapse, which is a top down process, and the bottom-up mechanism of structure formation. This latter mechanism is regulated by the properties of density fluctuations (in our case a simple Poisson distribution).

As final remark we note that the question of how many particles are in practice needed to simulate accurately the collisionless limit (up to some specified time) can only be answered in the context of a given problem and it is in general a very difficult task to be sorted out. For instance, it was recently shown through the study of the evolution a simple class of initial conditions (initially spherical density profiles), that there is a distinct $N$ dependence associated with the presence of instabilities in the collisionless dynamics that arises because the initial seeds for the instability themselves depend on $N$ Benhaiem et al. (2018): this dependence on $N$ is very difficult to find, as it manifests itself only in a very weak dependence of the time of triggering of the instability, and not, at sufficiently large $N$, in the properties of the state to which the instability drives the system.

The continuum limit for the dynamics of our finite $N$ systems is given by the Vlasov-Poisson system and in principle could be simulated numerically directly. While much progress has been made on the numerical solution of these equation (see e.g. Sousbie and Colombi (2016)), such an approach is, for the present, feasible numerically only for systems with high degrees of symmetry, and not for those here in which the breaking notably of rotational symmetry plays a crucial role.

IV Models and observed structures

Both our initial conditions and the dynamics of the systems we are studying are highly idealized. In particular we expect that in the formation of most astrophysical structures that non-gravitational processes will play a crucial role (e.g., gas dynamics, star-formation and feedback from it in galaxy formation). Thus our models are intrinsically not suitable to provide a detailed quantitative model for the formation of astrophysical objects. On the other hand, focusing on the specific case of galaxy formation, we note that there is in fact little direct constraint on initial conditions for it, as the fluctuations measured in the cosmic microwave background constrain strongly only significantly larger scales. Further the relative importance of gravity and other forces in shaping galaxies is very uncertain. As the structures we have seen in our simulations bear a striking resemblance to spiral galaxies, we believe it is worth looking more carefully at whether the qualitative features of these structures are compatible or incompatible with the observed qualitative properties of galaxies.

IV.1 Morphological features

Firstly we note that there are several very common and non-trivial features of spiral galaxies, which are accounted for in this kind of model while they are problematic in the usual theoretical approaches, in which spiral structures emerge through instability of an equilibrium disc (see, e.g., Dobbs and Baba (2014); Binney and Tremaine (2008)).

Our models lead, as we have seen, very easily to two armed spiral structure, which is observationally the most common kind. Its predominance has been considered puzzling and is difficult to account for in the usual theoretical approaches. Further the “pitch angle” $\alpha$, defined as the angle between the tangent to the arm and the tangent to a circle at the same angle, gives values in the range $10^{\circ}-40^{\circ}$ in our simulations (at $t=50$) except for the very cold cases like A1 in which the particles are ejected with very high radial velocities (leading to an $\alpha$ approaching $90^{\circ}$). Pitch angles of this order are typical observationally, while theoretical models predict much smaller angles and have great difficulty accounting for those observed Dobbs and Baba (2014). Further, as is almost invariably the case observationally, the spiral arms formed in our models are trailing, i.e. the outer tip points in the direction opposite to rotation (Binney and Tremaine, 2008). As we have described above, the spiral-like structure arises precisely because the particles which are furthest out have, by angular momentum conservation, smaller transverse velocities and thus lag behind in the angle they rotate through up to a given time. This is again an observational fact which does not have an apparent explanation in the usual theoretical approaches. Finally our mechanism produces, by construction, structures which are non-axisymmetric and often the central core is bar-like. Further, the spiral arms start at the end of the bar. These properties likewise appear not only to be compatible with observations (see, e.g., Jog and Combes (2009); Dobbs and Baba (2014)), but to potentially explain them in a very natural manner which eludes the usual theoretical approaches.

It is interesting to note that the mass of stars and gas in the spiral arms and in the outer part of the disc in general do not exceed the $\sim 20\%$ of the luminous mass of the galaxy (see, e.g., Sofue (2017)). In our simulations we also find that most of the mass is concentrated in the central bulge and that the fraction of particles with energy close to, or larger than, zero represent a small fraction, typically of the order of $10\div 20\%$, of the total mass.

IV.2 Time-scales

We have discussed in Benhaiem et al. (2017) some simple considerations about the compatibility of time and length scales with real spiral galaxies, for the most idealized case of a single ellipsoidal cloud. Making simple assumptions linking the final size and velocity scale to those of real galaxies, the collapse process which generated the structure must be assumed to occur on a timescale of the order of $\sim$ 1 Gyr, that is the characteristic time scale of all out of equilibrium transient structures that are formed in our simulations, i.e. spiral arms, bars and ring-like structures. This is much shorter than the age of the oldest stars in these galaxies ($\sim$ 10 Gyr) which is usually assumed to correspond also to the age of such structures. From the observational point of view, however, there is no definitive evidence establishing the age of spiral arms; rather some observations have suggested that spiral arms are not long-lived (see e.g. Dobbs and Baba (2014); Binney and Tremaine (2008) and references therein). Indeed the oldest stars and the galaxy are formed by very different dynamical processes occurring on very different length scales (the size of the clouds where star formation occurs is of the order of $10^{-2}-10^{-1}$ kpc, while the size of a typical galaxy is $10-10^{2}$ kpc) and thus it is not at all evident that these two time scales must be of the same order of magnitude.

The second family of IC we have studied here illustrate clearly that this time-scale inferred from the space and velocity scales represents only that of the violent collapse leading to this outer structure, which could quite possibly be dissociated from that of the formation of the central part of the galaxy, which could occur much before a secondary collapse of surrounding matter giving rise to the disc and extended “halo” structure we have described.

In this respect it is perhaps useful to recall that the usual assumption in modelling galaxies as dynamical equilibria (i.e. as QSS) is intimately linked to these considerations of time and length scales. Indeed if we suppose a star orbits the galactic centers at a distance $R$, the number of revolutions it has made in a time $T$ is

[TABLE]

where $T_{10}$ is the time-scale in units of $10$ Gyr, $v_{200}$ is the velocity in units of $v=200$ km/sec, and $R_{10}$ the radius in units of $10$ kpc.

For stars (or other emitters) moving on closed Keplerian orbits, with a circular velocity $v\sim\sqrt{GM/r}$), this assumption (of “stationarity”) appears to be reasonable only if $n_{rev}\gg 1$, i.e., if these bodies have characteristic crossing times in the system considerably longer than its estimated lifetime. For smaller values they cannot have had the time to attain orbits in which there is an equilibrium between centrifugal and centripetal acceleration. The precise value of the number of revolutions needed to reach a relaxed configuration cannot be constrained in a simple way theoretically because, as we have discussed above, it depends on the time in which the relaxation from an out of equilibrium configuration to a QSS takes place. However, from a qualitative point of view, a reasonable requirement is that $n_{rev}\gg 1$: only if this condition is satisfied can the assumption of stationarity possibly be justified.

For a typical disc galaxy of the size of the Milky Way, with a characteristic velocity $v_{200}\sim 1$ the number of revolutions is $n_{rev}\geq 10$ for $R\leq 30$ kpc only if $T_{10}\sim 1$, i.e. if the age of the galaxy structure is of order of the oldest stars. If $T_{10}\sim 1$, i.e. if the age of the galaxy structure is of order of the oldest stars then objects in the inner part of the galaxy (i.e., $R<10$ kpc) had enough time to make $n_{rev}\geq 10^{2}$ while at larger distances $n_{rev}\leq 10$ and thus the assumption that emitters move on closed Keplerian orbits appears very difficult to justify for the outermost regions of the galaxy. On the other hand, if the age of the galaxy structure is $T_{10}\sim 0.1$, the assumption of stationarity clearly cannot be valid at larger radii because $n_{rev}\ll 10$.

The out of equilibrium scenario of our models for the outer parts of a galaxy of the size of the Milky Way is then, however, coherent with the observed velocity and distance scales.

IV.3 Velocities

The compatibility of the velocity space structure of our transient structures with observed properties of galaxies is much less evident. Indeed a generic feature of the structures in our simulations, arising from the nature of the mechanism, is that velocities becomes predominantly radial at large radii. Extensive observational study over decades, using different tracers, has placed much constraint on such motions (Sofue, 2017), and indicates that motion in the outer parts of such galaxies is in fact very predominantly rotational (Kalberla and Dedes, 2008; Sofue, 2017), although a significant radial motions have been detected in many objects (López-Corredoira and González-Fernández, 2016). As discussed in Benhaiem et al. (2017) for the results based on simple ellipsoidal IC, it turns out that the naive expectation that such motions are excluded by observations is not confirmed. The reason is that our velocity fields have a very particular spatial (anisotropic) spatial structure which makes it difficult to distinguish them in projection from rotating disc models. For the broader class of IC we have explored here, the same considerations are valid, as there is a similar kind of correlation between the velocities and the spatial configuration. Further, we have seen that different IC can produce a less violent evolution than in the pure ellipsoidal model, leading to less radial motion and a much more extended region in which there is predominantly rotational motion. Thus the compatibility with observations of galaxy kinematics depends also on the identification of the time and length scales of the models with those of real galaxies. In our conclusion below we will comment about the radial motions observed in our own Galactic disc, that are relevant to understand the possible non stationary nature of the outer parts of the disc and of the spiral arms.

V Dissipationless and dissipative disc formation

Let us now briefly discuss the difference between the models that we have presented, where the formation of a disc like flat structure is originated solely by a gravitational, and thus dissipationless, dynamics and models in which instead the formation of a disc like structure is driven by dissipative effects. In standard models of galaxy formation the key element in the formation of a disc galaxy, is the dissipation associated with non-gravitational processes — gas dynamics, cooling, star formation, etc. The models used when gas dynamics is added to gravitational physics consider the collapse of isolated and rotating clouds, like the one we studied here, but solely with a spherical initial configuration. The central finding of our study is that disc-like configurations with transient spiral arms and with bars and/or rings in our simulations are formed by a purely dissipationless gravitational dynamics if the initial conditions break spherical symmetry. There have been attempts Katz (1991) to study the formation of quasi equilibrium configurations through a purely gravitational and dissipationless collapse dynamics in which the initial conditions are represented by isolated, spherically symmetric top hats in solid body rotation and in Hubble flow. These initial conditions differ from the ones we considered in this work by (i) the initial spherical shape and (ii) the small-scale fluctuations, intended to model the fluctuations in standard “cold dark matter cosmologies”. Starting from these initial conditions the QSS formed are slowly rotating and are supported by an anisotropic velocity dispersion and closely resemble elliptical galaxies, and do not resemble at all spiral galaxies.

The seminal work by Gott (1977) described a scenario — then developed in many other subsequent works — in which elliptical galaxies are products of a purely gravitational dissipationless collapse at high redshift, while spirals formed later with considerable dissipation. To simulate such a scenario dissipative gas dynamics was included in numerical simulations Katz and Gunn (1991) in a system with the same class of initial conditions as in Katz (1991). It is precisely the dynamics of the gas in this two-component system, which leads then to structures resembling spiral galaxies: a thin disc made of gas and surrounded by purely gravitational matter. Indeed since the gas can shock and dissipate energy it can develop a much flatter distribution than the dissipationless matter. By adding to the same initial conditions of Katz (1991), other non-gravitational effects, as star formation, supernova feedback Katz (1992) and metal enrichment due to supernovae Steinmetz and Mueller (1994); Steinmetz and Muller (1995), the crucial mechanism for the formation of the disc and of spiral arms is again played dissipative non-gravitational processes. These scenarios are thus completely different to how this structure emerges in our simulations, which are pure gravitational and dissipationless.

VI Discussion and conclusions

We have described the results of numerical experiments exploring the evolution under their self-gravity of non-spherical uniform and non-uniform clouds with a coherent rigid-body rotation about their shortest semi-principal axis. We have focused in particular on the very rich spatial and velocity space organization of the outer parts of these structures, which is a result of the combination of the violent relaxation, which leads to a high energy tail in the energy distribution, and the coherent rotation. Under very general conditions spiral-like structure arises, while more or less evident bars and/or rings appear depending on the properties of the initial conditions. These outer parts of the structure are intrinsically non-stationary and continue indefinitely to evolve in time. Although it will disappear asymptotically, such structure is very long lived and, in most cases, is still clearly defined at the longest times we simulate to, of order $100$ to $200$ times the dynamical time (characterizing the time for the formation, and characteristic time, of the virialized core). The particles forming the spiral arms will escape from the system, if they have positive energy, or will form an extremely dilute and anisotropic halo; eventually some of them, those having negative energy, will return back toward the core. On the other hand the largest fraction of the mass is bound in a triaxial system.

It is perhaps relevant to remark why this simple path to producing such structure in a self-gravitating system has, apparently, been overlooked in the literature. Indeed the dynamics of self-gravitating clouds of various forms, and with a wide variety of initial velocity distributions, with and without initial angular momentum, has been studied in depth in the literature, and it may seem surprising that the phenomena we have discussed have not been noticed. We believe that the explanation is probably linked to, on the one hand, the small fraction of mass involved, and, on the other hand, the relatively long times scales on which the system must be monitored. Indeed most studies of this kind of system focus on the relatively short times on which the system appears to virialize as indicated by global parameters. Further most studies of this kind date back two to three decades, and simulations in which the number of particles was typically of $\sim 10^{4}$. In this case, the high energy particles which are typically of order $5-10\%$ are too few to resolve the structures we have studied.

The spatial distributions of these transient structures are all “spiral-like” , but even within the very circumscribed and idealized set of IC we have considered, they show a wide variety of forms, from ones qualitatively resembling grand design spiral galaxies, to multi-armed and flocculent spiral galaxies and barred spiral galaxies. The mechanism for producing these structures is completely different in its physical principle to the mechanisms widely considered as potentially explaining observed spiral structure. Indeed while such mechanisms treat the spiral structure as a perturbative phenomenon — produced by the perturbation of an equilibrium rotating disc (Hohl, 1971; Zang and Hohl, 1978; Sellwood, 1985; Sellwood and Moore, 1999; Binney and Tremaine, 2008; Fujii et al., 2011; Dobbs and Baba, 2014). — the mechanism at play in our simulations is intrinsically far out-of-equilibrium. While our model is too simple and idealized to provide a quantitative model for real spiral galaxies, we have noted that, in many respects, it apparently reproduces very naturally many of their noted qualitative features.

For what concerns the problem of cosmological galaxy formation, we note that the monolithic collapse discussed here is compatible with a top-down structure formation of the kind that occurs, e.g., in the so-called warm dark matter models. On the other hand, in models where dark matter is cold, structure formation proceeds in a hierarchical bottom-up manner, so that galaxies are formed through an aggregation (i.e., merging) of smaller sub-structures. In this respect we note that when the initial conditions break spherical symmetry and are inhomogeneous (as our models C1 and C2), before the complete monolithic collapse of the whole cloud there are substructures forming and merging: in this scenario, because all substructures take part to the whole system collective dynamics, they can finally form coherent structures, like bars, rings and spiral arms, which are as large as the system itself. Thus, the scenario we have discussed is not in contradiction with the various observational evidences that merging was efficient in the early universe, but clearly a comprehensive theory of cosmological structure formation must be specified by the whole power-spectrum of density fluctuations: this does beyond the scope of the present work but will be addressed in forthcoming papers.

We will explore in future work some of the questions opened up by our results. One such question is of course whether the kind of initial conditions we assume could be produced easily within a cosmological framework. As mentioned above, the problem of collapsing clouds has been wildly studied in the context of cosmological galaxy formation: however these studies were performed by taking a spherical overdensity while we used here more general shapes and a purely gravitational (and thus dissipationless) dynamics. In addition, we note that while this seems to be excluded in typical scenarios in which structure formation proceeds hierarchically from very small scales (e.g., cold dark matter type scenarios), conditions like those we assume might possibly be plausible in the case in which initial fluctuations are highly suppressed below some large scale (e.g., as occurs in warm dark matter type scenarios). A different but complementary direction would be to explore the additional effects of non-gravitational and dissipative physics, like gas dynamics, modeling the complex processes inevitably at play in galaxy formation, and how they may or may not modify formation and evolution of the structures we have focused on here.

Finally it is interesting to mention that, while it has been known for several decades that the disc of the Milky Way contains large-scale non-axisymmetric features, the full knowledge of these asymmetric structures and of their velocities fields is still lacking. The recent Gaia DR2 maps (Gaia Collaboration et al., 2018) have clearly shown that the Milky Way is not, even to a first approximation, an axisymmetric system at equilibrium, but that it is characterized by streaming motions in all three velocity components. In particular it has confirmed the coherent radial motion in the direction of the anti-center, earlier detected by López-Corredoira and González-Fernández (2016), up to 14 kpc. In addition, recent analyses of the radial velocity field in our Galaxy by using different data-sets López-Corredoira et al. (2019); López-Corredoira and Sylos Labini (2019), provide lots of new and corroborated information about the disk kinematics of our Galaxy: significant departures of circularity in the mean orbits with radial galactocentric velocities, variations of rotation speed with position, asymmetries between Northern and Southern Galactic hemisphere and others (note that the analysis of the velocity fields in external galaxies is model dependent and thus also the estimation of radial motions Sylos Labini et al. (2019)). These features of the full three-dimensional velocity field seem to be compatible with the complex velocity fields generated by the gravitational collapses we have discussed but a more detailed comparison of the models and observations is needed. Such analysis will be reported in a forthcoming work.

Acknowledgements.

This work was granted access to the HPC resources of The Institute for scientific Computing and Simulation financed by Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX- 29-01) overseen by the French National Research Agency (ANR) as part of the Investissements d’Avenir program. FSL thanks Martin López-Corredoira for very useful discussions and comments.

.