Transient Spiral Arms from Far Outofequilibrium Gravitational Evolution
Abstract
We describe how a simple class outofequilibrium, rotating, and asymmetrical mass distributions evolve under their selfgravity to produce a quasiplanar 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 twoarmed 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 subdominant 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 nontrivial nonstationary 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.
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 outofequilibrium dynamics of purely selfgravitating 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 largescale motions are modeled assuming a stationary mass distribution. Indeed, our study illustrates how, for intrinsically nonstationary 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 nonKeplerian 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 nonaxisymmetric 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 darkmatter 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 spirallike 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 selfgravitating clouds with some angular momentum. The dynamics of isolated selfgravitating systems from outofequilibrium 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. subvirial, 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 subvirial 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 meanfield, 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 threedimensional and twodimensional 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 subvirial (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 bodylike rotational motion about the shortest semiprincipal 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 darkmatterdominated 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 quasiuniform 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 semiaxes 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 semiprincipal axis, and is the initial gravitational potential energy ^{1}^{1}1Because the forcesmoothing 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 selfgravity until a time where is the characteristic time scale for their meanfield evolution defined as
(1) 
where is the initial mass and is the Newton’s constant. All simulations^{2}^{2}2See Joyce et al. (2009); Sylos Labini (2012); Benhaiem & Sylos Labini (2015, 2017) for details. are performed for particles, using the gadget2 code (Springel et al., 2001), adopting a forcesmoothing which is approximately onetenth 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 Threedimensional properties
We observe, as expected given the chosen initial velocity distribution and normalization, a significant contraction and a subsequent reexpansion 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 semimajor 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 spirallike 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 cases^{3}^{3}3Here, 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.
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 ^{4}^{4}4see goo.gl/L1fRzZ for the full movies of the time evolution. projected on the plane orthogonal to the initial shortest semiprincipal 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 spirallike arms expand radially, slowly changing shape.
Indeed the velocity field of the particles in the outer part of the object is almost radial and directed outward (see Fig. 3).
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 powerlaw 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 ).
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.
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 .
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 angleaveraged data do not give information about the angular dependence of the radial velocities in particular, which are very nontrivial: the presence of the spirallike 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 semiprincipal 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 highenergy particles at the time of collapse compared to their transverse velocity at this time (and thus, in particular, on the parameter ). The nonstationary 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 longlived (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 largescale 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ópezCorredoira & GonzálezFerná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 hypervelocity 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 onedimensional (i.e., along the major axis of the observed galaxy) and twodimensional (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 nonaxisymmetric 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 semiprincipal 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 nontrivial rotation curve from a onedimensional measurement.
For twodimensional 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 twodimensional data. In particular, twodimensional 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 semimajor 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. ErrozFerrer 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 semiprincipal 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 semiprincipal axis. To fit the resulting twodimensional 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 twodimensional 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 twodimensional LOS velocity dispersion map; the largest dispersion is in the core where the velocities are isotropic;

(iii) the twodimensional LOS velocity map;

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

(v) the twodimensional LOS velocity map of the bestfit rotating disk model (this is obtained by using the onedimensional LOS velocity profile along the estimated kinematic axis);

(vi) the twodimensional LOS velocity residual map.
Figure 8 shows, respectively, the onedimensional 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).
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 nonaxisymmetric 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 threedimensional 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 subdominant, radial motions and nonaxisymmetric 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 threedimensional 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.
The reason for this surprising result is the strong correlation in the alignment of the kinematic axis and the longest semiprincipal axis of the projected distribution, which is a characteristic of our outofequilibrium 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. nonaxisymmetric distributions with predominantly radial velocities very nontrivially 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 socalled Modified Newtonian Dynamics (Milgrom, 1983, 2016). Indeed, scaleindependent radial motions would give an inferred scaleindependent rotation curve, and thus where is the inferred constant velocity of rotation, while where is the mass in the virialized core. Thus,
where
The observed approximate constancy of for different systems then corresponds to , i.e. the TullyFisher relation.
4 Discussion
In summary, we have shown, using simulations of evolution from very simple toy initial conditions, that transient spirallike structure may be generated in the far outofequilibrium evolution of a relaxing selfgravitating system. As will be detailed further in a forthcoming work (Benhaiem et al., 2017), the spatial organization in a spirallike 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 nonaxisymmetric 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 nonstationary 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 nongravitational 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.
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 semiprincipal axis is then designated by a unit vector parallel to , the intermediate semiprincipal axis by a unit vector parallel to , and the shortest semiprincipal 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 nonaxisymmetric 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
(1) 
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
(2)  
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
(3) 
Using the expressions above, a little algebra gives
(4)  
The position coordinates of the particles in the plane of projection, and projected velocity , can then be calculated, for any given observer , as
(5)  
A.4 Onedimensional Apparent Velocity Profiles
Most observations of apparent velocities are not fully twodimensional, but given along a specific axis (corresponding to the orientation of the slit used for the spectographic measurements). In order to obtain such onedimensional 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 threedimensional 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
(6) 
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
(7) 
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
(8)  
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 threedimensional major axis of the nonaxisymmetric 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., ErrozFerrer et al. (2015)).
To find the bestfitting 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:
(9)  
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 bestfitting inclination angle, we minimize the sum of the residuals in all the cells with respect to , i.e.
(10) 
References
 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., LyndenBell, 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
 ErrozFerrer et al. (2015) ErrozFerrer, 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ópezCorredoira & GonzálezFernández (2016) LópezCorredoira, M., & GonzálezFerná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 eprints, 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