Transient spiral arms from far out of equilibrium gravitational evolution

Transient Spiral Arms from Far Out-of-equilibrium Gravitational Evolution

David Benhaiem Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, Via dei Taurini 19, I-00185 Roma, Italia Michael Joyce Laboratoire de Physique Nucléaire et de Hautes Énergies, UPMC IN2P3 CNRS UMR 7585, Sorbonne Universités, 4, place Jussieu, 75252 Paris Cedex 05, France Francesco Sylos Labini Centro Studi e Ricerche Enrico Fermi, Via Panisperna 00184, Roma, Italia Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, Via dei Taurini 19, I-00185 Roma, Italia INFN Unit Rome 1, Dipartimento di Fisica, Universitá di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italia
July 24, 2019yyyyzzzzJuly 24, 2019
July 24, 2019yyyyzzzzJuly 24, 2019
July 24, 2019yyyyzzzzJuly 24, 2019
July 24, 2019yyyyzzzzJuly 24, 2019

We describe how a simple class out-of-equilibrium, rotating, and asymmetrical mass distributions evolve under their self-gravity to produce a quasi-planar spiral structure surrounding a virialized core, qualitatively resembling a spiral galaxy. The spiral structure is transient, but can survive tens of dynamical times, and further reproduces qualitatively noted features of spiral galaxies such as the predominance of trailing two-armed spirals and large pitch angles. As our models are highly idealized, a detailed comparison with observations is not appropriate, but generic features of the velocity distributions can be identified to be the potential observational signatures of such a mechanism. Indeed, the mechanism leads generically to a characteristic transition from predominantly rotational motion, in a region outside the core, to radial ballistic motion in the outermost parts. Such radial motions are excluded in our Galaxy up to 15 kpc, but could be detected at larger scales in the future by GAIA. We explore the apparent motions seen by external observers of the velocity distributions of our toy galaxies, and find that it is difficult to distinguish them from those of a rotating disk with sub-dominant radial motions at levels typically inferred from observations. These simple models illustrate the possibility that the observed apparent motions of spiral galaxies might be explained by non-trivial non-stationary mass and velocity distributions without invoking a dark matter halo or modification of Newtonian gravity. In this scenario the observed phenomenological relation between the centripetal and gravitational acceleration of the visible baryonic mass could have a simple explanation.

galaxies: formation — Galaxy: formation —methods: numerical — Galaxy: kinematics and dynamics — galaxies: spiral
journal: ApJ

1 Introduction

The arms of spiral galaxies are one of the most striking and remarkable features of the visible universe revealed by astronomy. They have been the subject of much study, both observational and theoretical, over many decades. Several competing theories have been advanced to explain their physical origin, but no single one has emerged definitively as the correct framework (see, e.g., Dobbs & Baba (2014)). Understanding of their motions is of particular importance because it is the observed apparent (i.e., on the line of sight — LOS) motions in the outer parts of spiral galaxies that have led to the supposition that much of the gravitating matter in them is not visible (Rubin, 1983). These same motions have led also to alternative scenarios involving strong modifications of Newtonian gravity (Milgrom, 1983). In this paper we show how mass distributions qualitatively resembling those of the visible components of spiral galaxies can result from the far out-of-equilibrium dynamics of purely self-gravitating systems, starting from a class of very simple idealized initial conditions. We study in particular the generic features of the velocity distributions of the structures produced by this mechanism, and consider their qualitative compatibility with observations of motions in spiral galaxies.

Our approach is different from standard theoretical ones, in which spiral structure arises by perturbation (internal or external) of an equilibrium system, and the large-scale motions are modeled assuming a stationary mass distribution. Indeed, our study illustrates how, for intrinsically non-stationary models, the relation between apparent motions and the associated mass distribution can be completely different from that in stationary models. In particular, we show that the observation of a non-Keplerian rotation curve in the outer part of such a structure does not necessarily require the existence of an extended dark matter halo or modification of Newtonian gravity, and could instead be consistent with non-axisymmetric radial motion of weakly bound and unbound mass.

We note that, because our models involve only Newtonian gravity, the physics we describe could potentially be applicable to astrophysical systems of very different natures and sizes — to dwarf galaxies that are inferred from their motions to be even more dark-matter dominated than spirals (see, e.g., Combes (2002)); to protoplanetary disks, which have been revealed in observations in the last couple of years to have spiral-like structure (see, e.g., Christiaens et al. (2014)); or even possibly to circumplanetary disks, whose existence is still inconclusive (see, e.g., Ward & Canup (2010)). In a forthcoming work (Benhaiem et al., 2017) that is complementary to this paper, we will describe the physical mechanism in much greater detail, using both for a broad range of initial conditions and also numerical simulations with larger particle numbers.

The class of models we consider as initial conditions consists of asymmetrical and isolated self-gravitating clouds with some angular momentum. The dynamics of isolated self-gravitating systems from out-of-equilibrium initial conditions has been extensively studied for several decades (Henon, 1973; van Albada, 1982; Aarseth et al., 1988; David & Theuns, 1989; Aguilar & Merritt, 1990; Theuns & David, 1990; Boily et al., 2002; Barnes et al., 2009). Broadly speaking, such systems relax quite efficiently to virial equilibrium, i.e., on time scales of the order of a few times the characteristic dynamical time. Early studies showed that spherical configurations with little isotropic velocity dispersion (i.e. sub-virial, with an initial virial ratio ) could produce equilibrated structures resembling elliptical galaxies, with a surface brightness notably close to the observed de Vaucouleurs law (van Albada, 1982). One generic feature of such sub-virial collapses (for ) is that they lead to the ejection of some of the initial mass (see, e.g., Joyce et al. (2009); Sylos Labini (2012, 2013)): the strong contraction of the initial configuration leads to a rapidly varying mean-field, which causes particle energies to also rapidly vary, leaving some of them weakly bound and others with positive energy. Two of us have recently studied (Benhaiem & Sylos Labini, 2015, 2017) the evolution from configurations that are initially ellipsoidal or of an irregular shape and found them to give rise to a virialized central core surrounded by very flattened configurations made by both weakly bound and ejected particles.

These results, combined, have led us to the idea that, with some initial rotational motion, it might be possible to generate a spiral structure from these kind of initial conditions. Indeed, the large radial velocities are generated in a small region (of the order of the minimal size reached) in a very short time (much less than one dynamical time), and thus the radial distance these particles subsequently travel once they are outside the core can be expected, given approximate conservation of angular momentum, to be correlated with the integrated angle they move through.

The paper is organized as follows. In Sect.2 we present the details of our numerical simulations. Sect.3 is then devoted to a discussion of the three-dimensional and two-dimensional results of our simulations and their relation with some key observational results on spiral galaxies. Finally, in Sect. 4 we draw our main conclusions. In the Appendix we detail how we constructed the projected velocity maps from our simulated mass distributions.

2 Simulations

We have considered a very simple set of initial conditions that combines the characteristics described above: breaking of the spherical symmetry of the initial mass distribution, a velocity distribution that

is sub-virial (or, more generally, out of equilibrium), and some coherent rotation. More precisely, we consider the following: particles distributed randomly, with uniform mean density, inside an ellipsoidal region, and velocities which correspond to a coherent rigid body-like rotational motion about the shortest semi-principal axis. Although these are ad hoc and clearly too idealized to describe a physically realistic situation, in the context of the theory of galaxy formation these kinds of initial conditions have often been argued to be reasonable (see, e.g., Eggen et al. (1962)). Nevertheless, they are very different from those described in current scenarios for galaxy formation in the context of cold dark-matter-dominated cosmological models, which are characterixed by hierarchical collapse. We note, however, that in cosmological scenarios with very suppressed initial fluctuations at very small scales (e.g. in models with warm dark matter), a monolithic collapse from a quasi-uniform initial state may be a more reasonable approximation. In any case, our goal here is to identify and study a physical mechanism and its possible observational signatures, and not to provide a realistic modeling of great complexity.

The parameters we choose to characterize our initial configurations are then (i) the ratios of the semi-axes of lengths : the ellipsoids can be prolate, oblate, or triaxial and they are specified by the flatness parameter ; (ii) the initial virial ratio , where is the kinetic energy of the rotational component of the motion, which has an angular velocity independent of radius (i.e. solid body rotation) and parallel to the shortest semi-principal axis, and is the initial gravitational potential energy 111Because the force-smoothing at small scales is a factor of 10 smaller than the initial interparticle distance, the difference between (computed using the Newtonian potential) and the Clausius virial term (computed using the exact forces acting on particles) is negligible. We have explored a large parameter range in this family of initial conditions, extending down to which, although strictly “virial,” is well out of equilibrium for the chosen velocity distribution. We follow the evolution under self-gravity until a time where is the characteristic time scale for their mean-field evolution defined as


where is the initial mass and is the Newton’s constant. All simulations222See Joyce et al. (2009); Sylos Labini (2012); Benhaiem & Sylos Labini (2015, 2017) for details. are performed for particles, using the gadget-2 code (Springel et al., 2001), adopting a force-smoothing which is approximately one-tenth of the initial mean interparticle separation. In this paper we report in detail results for just one chosen simulation, whose features are representative of this class of models. Greater details on numerical issues and analyses of results for a broad representative range of these initial conditions, and also for a range of particle numbers extending up an order of magnitude larger, will be provided in a separate paper (Benhaiem et al., 2017).

3 Results

3.1 Three-dimensional properties

We observe, as expected given the chosen initial velocity distribution and normalization, a significant contraction and a subsequent re-expansion of the system on a time scale . Associated with this behavior is, as anticipated, also a strong injection of energy into a significant fraction of the particles, which are those initially located furthest from the center (i.e. close to the semi-major axis) and which pass through the center of the structure latest during the collapse. Correspondingly, we observe an amplification of the spatial asymmetry during this phase (with, in particular, a more rapid contraction along the shortest axis). In addition to these features, which have been studied extensively in previous works (Benhaiem & Sylos Labini, 2015, 2017), we find that these systems are qualitatively characterized in their outer parts by spiral-like structure, with a rich variety of forms — see Fig. 1 — ranging from some qualitatively resembling more grand design spirals, and others qualitatively resembling barred spirals and even flocculent spirals in some cases333Here, and in the following figures, we use units of length in which , and units of time in which ; energies are given in units in which .. Detailed analysis of the evolving configurations confirms that the emergence of this spatial organization— associated with a velocity distribution with very specific characteristic properties that we will describe below — is indeed the result of the injection of energy into some of the mass around the time of maximal contraction, which gives it large radial velocities in addition to the initial rotational motion.

Figure 1: Configurations resulting from four different initial conditions.

As anticipated above, we focus here, for simplicity, on the detailed analysis of just one specific initial condition, with and (i.e. a prolate initial ellipsoid), and . We choose this case because, even if it corresponds to a case that is not so far out of equilibrium characterised by a less violent contraction and expansion, it produces structure which is fairly typical of all cases. Shown in Fig. 2 are configurations of the evolved configuration at different times 444see for the full movies of the time evolution. projected on the plane orthogonal to the initial shortest semi-principal axis, along which the structure is (as expected) very flattened in extent compared to the observed projection: diagonalizing the inertia tensor to determine the principal axes and eigenvalues, we find a typical offset of a couple of degrees from the initial axes, but a much larger ratio for the eigenvalues, corresponding to a flatness parameter , while the core is triaxial with a flatness parameter and corresponds to a triaxial ellipsoid. We note that, once formed, the spiral-like arms expand radially, slowly changing shape.

Figure 2: Configurations resulting from four different times (see the labels). The solid line in the upper left panel corresponds to the initial ellipsoid.

Indeed the velocity field of the particles in the outer part of the object is almost radial and directed outward (see Fig. 3).

Figure 3: Configuration at : arrows are proportional to velocities

Fig. 4 shows the density profile , the velocity profile and the energy profile computed as averages in radial bins of constant logarithmic width. During the time evolution, the outer tail of is stretched to larger and larger distances. In general, when the system contraction during the collapse is strong enough to produce a large change of the particle energy distribution, the tail of the density profile is well fit by a power-law behavior with (Sylos Labini, 2013). Correspondingly the velocity and the energy profiles also extend to larger and larger scales. At the largest radii, as indicated by the average value particles are unbound (with ), while in the core region particles are strongly bound (i.e. below ); there is then an extended intermediate region in which many particles are marginally bound (i.e ).

Figure 4: From top to bottom: (i) density profile, (ii) velocity profile, and (iii) energy profile, at two different times: (black dots) and at (red dots).

The energy distribution at two different times () together with that of the initial conditions, is shown in Fig. 5: we note that a large change of has occurred during the gravitational collapse of the cloud at while at later times the shape of the distribution remains approximately the same.

Figure 5: Energy distribution at (black), at (red), and at (green).

The upper panel of Fig. 6 shows the average, in spherical shells of radius , denoted , of the radial component

of the velocity, and of the “transverse” velocity

defined parallel to the angular momentum relative to the origin (at the center of the structure). Thus, in particular, a coherent rotation of the shell in a plane corresponds to = .

The middle panel of Fig. 6 shows the velocity anisotropy

Finally, the lower panel of Fig. 6 shows , the mass that would be enclosed inside this radius if the motions were purely circular and the mass distribution spherically symmetric, and the mass actually enclosed inside the radius .

Figure 6: Configuration at . Upper panel: components of particle velocities averaged in spherical shells as a function of radius. Middle panel: anisotropy parameter . Lower panel: mass estimated from the velocity assuming stationary circular orbits, and the actual enclosed mass.

According to the behaviors observed, we can divide the structure into three regions: (i) an inner part (R1) in which, as

there is no significant net rotation, and given that , the velocity distribution is close to isotropic; (ii) an intermediate range of radii (R2), extending over about a decade, in which deviates strongly from zero as a net coherent rotational motion develops and dominates at larger radii, i.e.

correspondingly (lower panel of Fig. 6), there is a good agreement between the estimated and actual enclosed mass in this region; and (iii) an outer region (R3) in which the rotational motion of the particles is still coherent, but radial motions, with almost negligible dispersion, are now predominant, i.e.

Region R3 is also characterized clearly by the behavior of the estimated enclosed mass, which greatly overestimates the actual enclosed mass. This reflects the fact that the mass is weakly bound or even unbound rather than bound on circular orbits.

Measurement of the particle energies (see Fig. 4) shows that the transition from R2 to R3 is indeed approximately that from unbound to bound particles, and that in the outer part of R3 all particles are unbound. Indeed, asymptotically the knee between the two regions is precisely the transition from bound to unbound orbits, with the shell with corresponding to particles with zero energy. At large distance in R3 we have, correspondingly, a linear growth with distance of the radial velocity that is simply a reflection of the ballistic radial motion. Thus, when we study these curves as a function of time, region R1 and the inner part of R2 are stationary to a very good approximation, while the boundary with R3 propagates out progressively and the size of R3 itself grows linearly in time, with the maximal velocity remaining fixed.

These angle-averaged data do not give information about the angular dependence of the radial velocities in particular, which are very non-trivial: the presence of the spiral-like structure visible in Fig. 2 is a reflection of the fact that the particles’ transverse motions become correlated with their radial velocities because of the approximate conservation of angular momentum (and energy) of the ejected particles, which, once outside the core, move in an approximately stationary central potential. The particles forming the spiral structure preferentially have a radial velocity oriented along directions close to the initial longest semi-principal axis, and the structure is elongated the most along directions in which the radial velocities are maximal. Clearly, the precise form of the spiral structure depends directly on the dispersion of the energies of the high-energy particles at the time of collapse compared to their transverse velocity at this time (and thus, in particular, on the parameter ). The non-stationary nature of the structure also manifests itself in the evolution of the form of the spiral structure. In particular, it becomes more elongated (and less axisymmetric) in time.

3.2 Estimation of Typical Length/Time Scales

Even if our models are too simple and idealized to be meaningfully confronted with observations in any great detail, we can consider the qualitative compatibility of the features of the mass distributions generated with the observed properties of real astrophysical systems. In particular, we focus here on the primary astrophysical motivation for our study — spiral galaxies — although, as noted, several other applications could also be explored. For any such comparison we evidently need to relate approximately relate the scales of our toy model to physical scales. Bearing in mind that the typical scale of observed rotation velocities in disk galaxies is km/sec, we define the dimensionless parameters as follows:

We can then write

where is the number of dynamical times in our simulation and is its duration given in billion of years. Thus, for as in Fig. 2, and taking , which corresponds to a mass (by using Eq.1)

we have that region extends to , and region extends to . Thus, in order to have a structure that would possibly be compatible with the typical size of spiral galaxies, we need to assume that the collapse process which generated the disk and arms occurred much more recently than the formation of the oldest stars in these galaxies (with an age Gyr). This is very different from the usual hypothesis that the disk, and its spiral structure, are at least as old as the oldest stars. From the observational point of view, however, there is no definitive evidence establishing the age of spiral arms; rather several observational studies have suggested that spiral arms are not long-lived (Elmegreen et al., 1989; Vogel et al., 1993; Tully & Verheijen, 1997; Henry et al., 2003).

3.3 Characteristic of spiral arms

We note that the arms formed in our models are always trailing. This is a simple consequence of the approximate conservation of angular momentum for the outgoing particles, which means that the transverse components of their velocities decrease with their radial distance. Although, as mentioned, a rich variety of forms of the arms can be obtained with different initial conditions, two dominant arms as in our chosen simulation are very easily produced, with pitch angles of the order of tens of degrees. Thus, our model naturally reproduces very common features of spiral galaxies, which are very difficult to explain within the much explored framework of density wave theory (Dobbs & Baba, 2014), although density variations associated with spiral arms in our models are larger than they are in reality.

3.4 Apparent velocity maps

Let us now consider further the compatibility of large-scale motions of our generated mass distributions with observed apparent motions in disk galaxies. Depending on the initial conditions we choose, the details of the kinematic properties will change (e.g. the exact radial dependence of the velocities), but it is a generic property of this mechanism of generation of the spiral structure that there is a clear transition from predominantly rotational motion to predominantly radial motion, the latter being in the outermost parts the ballistic motion of freed particles. This is the case simply because the particles that are furthest from the center at long times are unbound or very loosely bound outgoing particles that have lost almost all their transverse velocity because of angular momentum conservation. Let us focus on this characteristic feature.

Decades of study of various different observational tracers of the velocity fields provide strong evidence for predominantly rotational motions in disk galaxies (Sofue & Rubin, 2001). For what concerns our Galaxy, in which apparent motions have been measured over four decades in scale (Sofue, 2017), the angular dependence of the projected velocities, inferred from HI emission in particular, shows convincingly that the motion of the disk is very predominantly rotational up to a scale of order kpc (Kalberla & Dedes, 2008; Sofue, 2017): as mentioned above, such a coherent rotation is also characteristic of the region R2 in our models. For this reason, the key observation thus concerns the nature of the motion at large distance, i.e. kpc, in our Galaxy. In this respect it is interesting to note that there is nevertheless also evidence for significant coherent radial motions beyond a few kpc and increasing with radius (López-Corredoira & González-Fernández, 2016). Beyond this scale the constraints are much weaker, but in the near future measurements from the GAIA satellite (Gaia Collaboration et al., 2016) will make it possible to distinguish the nature of the motions at much larger scales. In addition the GAIA satellite will be able to shed light on the nature of hyper-velocity stars that are unbound from the Milky Way and shows a surprising anisotropic distribution (Brown, 2015). A population of such stars might possibly correspond to ejected particles in our models.

Let us now consider constraints on velocity fields from external disk galaxies. In this case apparent (LOS) velocities are probed robustly out to scales of several tens of kiloparsecs, and in some cases even larger (Sofue, 2017), but the strength of evidence for rotation depends on the scale and weakens at larger scales. These measurements are both one-dimensional (i.e., along the major axis of the observed galaxy) and two-dimensional (mapping out the full projected velocity field). The former measurements provide a direct measure of rotational velocities, but only on the assumption that the galaxy is in fact a rotating disk: in this case the major axis of the projection (which is an ellipse) is orthogonal to the LOS, and thus motions parallel to the LOS must be rotational. In our models, the mass distribution is not a disk— indeed it is clearly non-axisymmetric at larger radii— and, furthermore, as we have noted, there is intrinsically a strong correlation of the direction of the outer radial velocities with the intrinsic longest semi-principal axis.

As a result, we show below that there is generically a contribution, which may be very large, along the projected major axis. In our models, as a result, even at length scales where the motion is purely radial, we will infer a non-trivial rotation curve from a one-dimensional measurement.

For two-dimensional data the evidence for predominant rotation (and the inferred rotation curves) is based on the quality of best fits to rotating axisymmetric disk models provided by two-dimensional data. In particular, two-dimensional velocity maps of numerous galaxies show the pattern distinctive of a rotating disk: the alignment of the kinematic axis — that along which there is maximal variation of the projected velocities — with the projected semi-major axis. Such alignment is, however, far from perfect and very significant angular offsets are frequently observed (and attributed to the breaking of axisymmetry by bars). Furthermore, very significant residuals are typically measured in such fitting procedures— typically of the order of or even larger— and these are attributed to radial motions (see e.g. Erroz-Ferrer et al. (2015)).

To qualitatively evaluate whether the radial motions that are dominant in the outer parts of the spiral structure in our models can be strongly excluded by observations, as one might naively expect, we have thus examined whether the projected motions of our toy galaxies can provide fits to rotating disk models of a comparable quality to those provided by the observed galaxies. To do so, as detailed in the Appendix, using our distributions we have generated the projected LOS velocity maps of random observers, characterised by two angles: the inclination angle , defined as the angle between the vector giving the orientation of the observer’s LOS and the vector in the direction of the shortest semi-principal axis of the model galaxy, and an angle defined as the angle between the projection of the LOS in the plane of the galaxy and its longest semi-principal axis. To fit the resulting two-dimensional map with a rotating disk model, for we determine the velocity as a function of distance along the axis of maximal variation, and use it as the input rotational velocity for a rotating disk, for which we analytically determine the projection. Shown in Fig. 7 are the projected velocity maps for the same simulation analyzed above, for an observer with and . The maps have been averaged on a grid of size (mimicking the finite resolution of measured maps); the different panels show the following:

  • (i) the two-dimensional projection of the mass distribution, with the kinematic axis and the major axis of the projection indicated; in this case the angle between the two axes is about ;

  • (ii) the two-dimensional LOS velocity dispersion map; the largest dispersion is in the core where the velocities are isotropic;

  • (iii) the two-dimensional LOS velocity map;

  • (iv) the two-dimensional LOS velocity map in which the radial velocities have been removed, illustrating that the motions are indeed very predominantly radial;

  • (v) the two-dimensional LOS velocity map of the best-fit rotating disk model (this is obtained by using the one-dimensional LOS velocity profile along the estimated kinematic axis);

  • (vi) the two-dimensional LOS velocity residual map.

Figure 7: Projection for (from top to bottom). Left panel: projection of the object on the observer’s sky; the kinematic axis (red) and the major axis (black) are shown. Middle panel: two-dimensional LOS velocity dispersion map. Right panel: two-dimensional LOS velocity map. Left panel: two-dimensional LOS velocity for the case in which the three-dimensional radial velocity has been set to zero. Middle panel: two-dimensional rotational map derived from the LOS velocity profile. Right panel: two-dimensional residual map.

Figure 8 shows, respectively, the one-dimensional LOS velocity profile along the kinematic axis and along the axis orthogonal to it (upper panel) and (lower panel) the mass estimated by assuming that the velocities are circular, and the actual mass (i.e., by using Eq.6).

Figure 8: Projection for . Upper panel: LOS velocity profile along the kinematic axis and along the axis perpendicular to it. Bottom panel: ratio between the mass estimated from the LOS velocity (assuming it to be circular and stationary) and the actual mass.

We have then explored (see Fig. 9) the full range of and . Only for very close to (i.e. an observer with an LOS almost exactly orthogonal to the axis along which radial velocities are maximal) do we fail to obtain a fit to a rotating disk model with residuals compatible with the level reported in the literature for such fits applied to observational data. These residuals are small in all cases, i.e. of the order of except for in which they can be as high as . In these images one can discern clearly that our model galaxies are non-axisymmetric at larger radii and, as we have noted, that there is a strong correlation of the direction of the outer radial velocities with the intrinsic major axis: the velocities in the outer parts of the structure are radial and very preferentially oriented along the axis which is significantly elongated in the structure. As a result there is generically a contribution from these radial velocities along the projected major axis. In addition, except for very small inclination angles, the projection of the three-dimensional major axis is typically very close to the major axis of the projected image, and the large radial velocities project out their component along this latter axis. There is thus in practice a rough degeneracy between rotating disk models with significant, but sub-dominant, radial motions and non-axisymmetric models with a specific pattern of radial velocities like the one in our models. This is very clearly illustrated by comparing in each case the map in which the three-dimensional radial velocity is set to zero and the map of the best fit rotational model: despite the fact that most of the signal at large scales comes from the radial velocities they can be fit quite well by the rotational model.

Figure 9: Residuals for and different values of from to in steps of .

The reason for this surprising result is the strong correlation in the alignment of the kinematic axis and the longest semi-principal axis of the projected distribution, which is a characteristic of our out-of-equilibrium structures: as we have seen, the velocities in the outer parts of the structure, which we are resolving in these mock measurements, are radial and very preferentially oriented along the axis, which is significantly elongated in the structure. In projection the major axis typically remains very close to this axis — other than for very specific observers, looking along the axis with very small inclination angles — and the large radial velocities project out their components along this axis.

A much more detailed and sophisticated analysis of observed projected velocity maps of spiral galaxies would evidently be required to establish or exclude their possible compatibility with velocity distributions qualitatively similar to those in our models, i.e. non-axisymmetric distributions with predominantly radial velocities very non-trivially correlated with the spatial distribution. As we have illustrated with our models, the motions in the outer parts of galaxies are in fact predominantly radial; there is no need to invoke a dark matter halo to explain them. Indeed, as illustrated in the lower panels of Fig. 6 and of Fig. 8, the mass estimate using the hypothesis of rotational motions leads to an inferred mass that grows strongly with radii, while the actual enclosed mass does not grow at all.

3.5 Flat Rotation Curves and Correlation between the Centripetal and Gravitational Acceleration

We conclude by speculating on two further important observational results about velocity fields, and their possible explanation within the framework suggested by our models.

One of the noted properties of rotation curves of spiral galaxies is that they are typically flat as a function of scale at the largest scales probed by observation, although a great variety of behaviors are in fact observed in individual galaxies (see, e.g., Sofue (2017)). Our models are not predictive in this respect: we can obtain very different behaviors depending on the range of scale considered, and notably whether we assume the region observed corresponds to R2 or R3. Furthermore, the precise functional dependence on scale may be very different if we modify, for example, the radial dependence of the initial angular velocity. We note, however, that, if we consider the region R3, in which radial motions dominate, the rotation curve (inferred by supposing the projected motions to arise from a rotating disk) will progressively flatten in time: as the velocities are essentially ballistic the same velocity range extends over a range of scale, which grows monotonically with time.

In this hypothesis of purely radial velocities with an approximately constant (i.e. very slowly increasing) amplitude, we note finally that one may also obtain very trivially in models like ours, also the observed phenomenological relation, , where is the centripetal acceleration inferred from the apparent motions, and the gravitational acceleration of the visible baryonic mass (see, e.g., McGaugh et al. (2016)), and which also underlies the so-called Modified Newtonian Dynamics (Milgrom, 1983, 2016). Indeed, scale-independent radial motions would give an inferred scale-independent rotation curve, and thus where is the inferred constant velocity of rotation, while where is the mass in the virialized core. Thus,


The observed approximate constancy of for different systems then corresponds to , i.e. the Tully-Fisher relation.

4 Discussion

In summary, we have shown, using simulations of evolution from very simple toy initial conditions, that transient spiral-like structure may be generated in the far out-of-equilibrium evolution of a relaxing self-gravitating system. As will be detailed further in a forthcoming work (Benhaiem et al., 2017), the spatial organization in a spiral-like structure arises dynamically as particles that gain significant energy during an initial collective contraction and expansion of the system move outward, with the more energetic particles losing their transverse motion faster. The mechanism is completely different in nature from the usual perturbative mechanisms widely studied to explain such structure. Despite the unrealistically simplified nature of the models, we have argued that a qualitative comparison with observational data is possible: our models show that the mechanism will generate structures with velocity fields which have a very characteristic behavior. This is a transition to predominantly radial motion with very small dispersion in the outermost parts. Surprisingly, we have found that the projected motions of these regions can typically be quite compatible with a rotating disk model, up to residuals attributed to radial motion which are very significant but of the order typically found in fitting rotating disk models to observations. This suggests the possibility that these motions could be explained without invoking either dark matter or a modification of Newtonian gravity, which are unavoidable if these galaxies are modeled as stationary and rotating. Rather, these motions might be consistent with the purely Newtonian gravitational dynamics of the visible mass if the outer parts of the galaxy are far from stationary and the motions are predominantly radial and spatially correlated in a non-axisymmetric distribution, rather than rotational. Instead of providing a single predictive model, we have opened a Pandora’s box of models, a different framework — of completely non-stationary mass distributions— that must be compared in much greater details with observations. Any such model is obviously also very simplistic, not just because of the idealization of the initial conditions but also in that it neglects everything but gravitational dynamics. Any detailed quantitative model will of course necessarily need to consider more complex initial conditions and also incorporate non-gravitational physics. There are other obvious apparent shortcomings of the toy model. For example, (i) spiral arms correspond to modest variations in mass density, and (ii) the time scale for collapse, as we have discussed, must be assumed short compared to the ages of old stars. The former may plausibly be related to the low mass resolution we have used, while the latter constraint may change in more complex initial conditions. Nevertheless we believe it is remarkable and tantalizing that the simple framework we have discussed produces structures bearing so much qualitative resemblance to astrophysical objects, and suggesting the possibility of a different and simple explanation for their observed projected motions.

The authors of this work were 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. We detail here how we construct the projected velocity maps reported in in Sect.3 from our simulated mass distributions. This projection is defined for a random observer at infinity. It is convenient, in order to understand the dependence on the orientation of the observer’s LOS, to define this orientation with respect to the principal axes of the mass distribution. Having done so, it then straightforward to determine the projected velocities as a function of this orientation and the components of the position and velocity in the principal axes.

A.1 Principal Axes

We compute the inertia matrix of the mass distribution relative to an origin taken at the minimum of the gravitational potential. We then determine its eigenvalues , where , and corresponding eigenvectors . (The longest semi-principal axis is then designated by a unit vector parallel to , the intermediate semi-principal axis by a unit vector parallel to , and the shortest semi-principal axis by . The plane of the galaxy is then orthogonal to . We then rotate from our original Cartesian axes (of the simulation) to determine the components of the particle positions, , and their velocities, , in the new basis .

A.2 Orientation of the Observer

Following standard conventions (see, e.g., Beckman et al. (2004)) we define the inclination angle of the observer as the angle between his LOS and a vector orthogonal to the plane of the galaxy, which we take to be . Furthermore, as the galaxy is non-axisymmetric about this axis, we define an azimuthal angle as the angle between the projection into the galaxy plane of the LOS and the major axis. Thus, we write the unit vector parallel to the LOS as


A.3 Determination of Projected Velocities

To define the axes giving the observer’s plane of projection it is convenient first to define the set of axes


in the plane of the galaxy. The vector is thus parallel to the axis of the projection in the plane of the galaxy of the observer LOS, while is the axis in the plane of the galaxy orthogonal to the observer LOS.

The projected plane, orthogonal to the LOS, is then spanned by the unit vectors


Using the expressions above, a little algebra gives


The position coordinates of the particles in the plane of projection, and projected velocity , can then be calculated, for any given observer , as


A.4 One-dimensional Apparent Velocity Profiles

Most observations of apparent velocities are not fully two-dimensional, but given along a specific axis (corresponding to the orientation of the slit used for the spectographic measurements). In order to obtain such one-dimensional velocity profiles we define two such slits: one aligned parallel to the kinematic axis, i.e. the axis along which there is the maximum gradient of the LOS velocity (details below), and one orthogonal to this direction. We have also considered projections along the major axis and minor axis of the projected distribution (defined following a procedure analogous to that described above for the three-dimensional case). The slit is assumed to be rectangular, of a width which is a small fraction of the minor axis. From these LOS velocity profiles along the kinematic axis , we estimate the mass enclosed in the radius assuming that particles are in stationary circular orbits as


where the inclination angle is estimated from the projection as described below.

A.5 Velocities for Rotating Disk Model

If one models a galaxy as a disk, the projected LOS velocities can be written (see, e.g., Beckman et al. (2004)) as


where are polar coordinates in the plane of the projection, with defined relative to the axis orthogonal to the observer LOS (i.e. parallel to the axis defined above), and and are the components of the velocity field given in polar coordinates in the plane of the galaxy (with defined relative to the same axis , which is also in the plane of the galaxy). The polar coordinates are related by the transformation


For a purely rotating axisymmetric model, and . The kinematic axis is that along which there is maximal variation of the projected velocity, i.e. .

A.6 Fitting to a Rotating Disk Model

To fit our projected velocity data to a rotating axisymmetric disk we first estimate from our data the orientation of the kinematic axis. We determine the kinematic axis as the axis passing through the center of mass of the distribution and along which the difference of the velocities at the two extreme points is maximal.

While this axis must strictly be the major axis of the projection if the underlying distribution is really a disk, this is generally not the case for our distributions that are not axisymmetric. However, because in our models the directions of the radial velocities are strongly correlated with the real three-dimensional major axis of the non-axisymmetric distribution (see below), the offset between the kinematic axis and the projected major axis is, in fact, typically (i.e. for a large fraction of random observers) quite small. Such offsets are, indeed, typically seen in observations (see, e.g., Erroz-Ferrer et al. (2015)).

To find the best-fitting rotating disk model, we need to determine the inclination angle : we do this by minimizing the residuals between the rotational model, computed for a generic , and the actual data on each grid cell. To do so we compute first, for each grid cell, labeled by and centered on projected coordinates , the polar coordinates as defined above:


Then, for the given value of the inclination angle , we use Eq.7 (with ) to compute the LOS velocity of the rotational model, denoted . Note that in the case where the unprojected size of the galaxy is larger than the maximum distance at which the LOS velocity profile extends, we perform a linear fit over the last five points of and then extrapolate using this fit to a higher radius. Finally, in order to get the best-fitting inclination angle, we minimize the sum of the residuals in all the cells with respect to , i.e.



  • Aarseth et al. (1988) Aarseth, S., Lin, D., & Papaloizou, J. 1988, Astrophys. J., 324, 288
  • Aguilar & Merritt (1990) Aguilar, L., & Merritt, D. 1990, Astrophys. J., 354, 73
  • Barnes et al. (2009) Barnes, E. I., Lanzel, P. A., & Williams, L. L. R. 2009, Astrophys. J., 704, 372
  • Beckman et al. (2004) Beckman, J. E., Zurita, A., & Vega Beltrán, J. C. 2004, Lecture Notes and Essays in Astrophysics, 1, 43
  • Benhaiem et al. (2017) Benhaiem, D., Joyce, M., & Sylos Labini, F. 2017, submitted
  • Benhaiem & Sylos Labini (2015) Benhaiem, D., & Sylos Labini, F. 2015, Mon.Not.R.Astron.Soc., 448, 2634
  • Benhaiem & Sylos Labini (2017) —. 2017, Astron.Astrophys., 598, A95
  • Boily et al. (2002) Boily, C., Athanassoula, E., & Kroupa, P. 2002, Mon. Not. R. Astr. Soc., 332, 971
  • Brown (2015) Brown, W. 2015, Ann.Rev.Astron.Astrophys., 52, 15
  • Christiaens et al. (2014) Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, Astrophys.J.Lett., 785, L12
  • Combes (2002) Combes, F. 2002, New Astronomy Review, 46, 755
  • David & Theuns (1989) David, M., & Theuns, T. 1989, Mon. Not. R. Astr. Soc., 240, 957
  • Dobbs & Baba (2014) Dobbs, C., & Baba, J. 2014, Pub.Astron.Soc.Austr., 31, e035
  • Eggen et al. (1962) Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, Astrophys.J., 136, 748
  • Elmegreen et al. (1989) Elmegreen, B. G., Seiden, P. E., & Elmegreen, D. M. 1989, Astrophys.J., 343, 602
  • Erroz-Ferrer et al. (2015) Erroz-Ferrer, S., Knapen, J. H., Font, J., & Beckman, J. E. 2015, Highlights of Astronomy, 16, 328
  • Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, Astron.Astrophys., 595, A1
  • Henon (1973) Hénon, M. 1973, Ann. Astrophys., 24, 229
  • Henry et al. (2003) Henry, A. L., Quillen, A. C., & Gutermuth, R. 2003, Astronom.J., 126, 2831
  • Joyce et al. (2009) Joyce, M., Marcos, B., & Sylos Labini, F. 2009, Mon. Not. R. Astron. Soc., 397, 775
  • Kalberla & Dedes (2008) Kalberla, P. M. W., & Dedes, L. 2008, Astron.Astrophys., 487, 951
  • López-Corredoira & González-Fernández (2016) López-Corredoira, M., & González-Fernández, C. 2016, Astron.J., 151, 165
  • McGaugh et al. (2016) McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Physical Review Letters, 117, 201101
  • Milgrom (1983) Milgrom, M. 1983, Astrophys.J., 270, 365
  • Milgrom (2016) —. 2016, ArXiv e-prints, arXiv:1609.06642
  • Rubin (1983) Rubin, V. C. 1983, Science, 220, 1339
  • Sofue (2017) Sofue, Y. 2017, Pub.Astron.Soc.Jap., 69, R1
  • Sofue & Rubin (2001) Sofue, Y., & Rubin, V. 2001, Ann.Rev.Astron.Astrophys., 39, 137
  • Springel et al. (2001) Springel, V., Yoshida, N., & White, S. D. M. 2001, New Astronomy, 6, 79
  • Sylos Labini (2012) Sylos Labini, F. 2012, Mon. Not. R. Astron. Soc., 423, 1610
  • Sylos Labini (2013) —. 2013, Mon. Not. R. Astron. Soc., 429, 679
  • Theuns & David (1990) Theuns, T., & David, M. 1990, Astrophys. Sp. Sci., 170, 276
  • Tully & Verheijen (1997) Tully, R. B., & Verheijen, M. A. W. 1997, Astrophys.J., 484, 145
  • van Albada (1982) van Albada, T. 1982, Mon. Not. R. Astr. Soc., 201, 939
  • Vogel et al. (1993) Vogel, S. N., Rand, R. J., Gruendl, R. A., & Teuben, P. J. 1993, Publ. Astron. Soc. Pac., 105, 666
  • Ward & Canup (2010) Ward, W. R., & Canup, R. M. 2010, Astronom.J., 140, 5
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description