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

David Benhaiem, Michael Joyce, Francesco Sylos Labini

► To cite this version:

David Benhaiem, Michael Joyce, Francesco Sylos Labini. Transient Spiral Arms from Far Out-of-equilibrium Gravitational Evolution. The Astrophysical Journal, 2017, 851 (1), pp.19. 10.3847/1538-4357/aa96a7. hal-01727320

HAL Id: hal-01727320

https://hal.sorbonne-universite.fr/hal-01727320

Submitted on 9 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

David Benhaiem¹, Michael Joyce², and Francesco Sylos Labini³,¹,⁴

¹Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, Via dei Taurini 19, I-00185 Roma, Italia

²Laboratoire de Physique Nucléaire et de Hautes Énergies, UPMC IN2P3 CNRS UMR 7585, Sorbonne Universités, 4, place Jussieu, F-75252 Paris Cedex 05, France

³Centro Studi e Ricerche Enrico Fermi, Via Panisperna I-00184, Roma, Italia

⁴INFN Unit Rome 1, Dipartimento di Fisica, Università di Roma Sapienza, Piazzale Aldo Moro 2, I-00185 Roma, Italia

Received 2017 August 7; revised 2017 October 7; accepted 2017 October 25; published 2017 December 6

## Abstract

We describe how a simple class of 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.

**Key words:** galaxies: formation – galaxies: spiral – Galaxy: formation – Galaxy: kinematics and dynamics – methods: numerical

## 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 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 also led 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 (D. Benhaiem et al. 2017, in preparation) 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 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 (Hénon 1973; van Albada 1982; Aarseth et al. 1988; David & Theuns 1989;

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 3.0 licence](https://creativecommons.org/licenses/by/3.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
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 timescales 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 $b > -1$) 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 $b \ge -0.5$) 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 kinds 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), 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 Section 2 we present the details of our numerical simulations. Section 3 is devoted to a discussion of the three-dimensional and two-dimensional results of our simulations and their relations with some key observational results on spiral galaxies. Finally, in Section 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: N particles distributed randomly, with uniform mean density, inside an ellipsoidal region, and velocities that 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 kind 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 characterized 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 $a_1 \ge a_2 \ge a_3$: the ellipsoids can be prolate, oblate, or triaxial and they are specified by the flatness parameter $\iota = (a_1/a_3) - 1$; (ii) the initial virial ratio $b_{\text{rot}} = 2K_{\text{rot}}/W_0$, where $K_{\text{rot}}$ 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 $W_0$ is the initial gravitational potential energy.<fnref n="5" /> We have explored a large parameter range in this family of initial conditions, extending down to $b_{\text{rot}} = -1$, which, although strictly “virial,” is well out of equilibrium for the chosen velocity distribution. We follow the evolution under self-gravity until a time $t \approx (50 \div 100)\tau_d$ where $\tau_d$ is the characteristic timescale for their mean-field evolution defined as

$$ \tau_d = \sqrt{\frac{\pi^2 a_3^3}{8GM}}, \qquad (1) $$

where M is the initial mass and G is Newton’s constant. All simulations<fnref n="6" /> are performed for $N = 10^5$ particles, using the gadget-2 code (Springel et al. 2001), adopting a force-smoothing that 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 the results for a broad representative range of these initial conditions, and also for a range of particle numbers extending up to an order of magnitude larger, will be provided in a separate paper (D. Benhaiem et al. 2017, in preparation).

## 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 timescale $t \sim \tau_d$. 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 longest semi-principal 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 semi-principal 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 Figure 1—ranging from some qualitatively resembling more grand design spirals, and others resembling barred spirals,

<footnote n="5">⁵ Because the force-smoothing at small scales is a factor of 10 smaller than the initial interparticle distance, the difference between $W_0$ (computed using the Newtonian potential) and the Clausius virial term (computed using the exact forces acting on particles) is negligible.</footnote>

<footnote n="6">⁶ See Joyce et al. (2009), Sylos Labini (2012), and Benhaiem & Sylos Labini (2015, 2017) for details.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Configurations resulting from four different initial conditions.</figcaption></figure>

and even flocculent spirals in some cases.<fnref n="7" /> 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 $t = 1$ and $a_2 = a_3$ (i.e., a prolate initial ellipsoid), and $b_{rot} = -1.0$. We choose this case because, even if it corresponds to a case that is not so far out of equilibrium and characterized by a less violent contraction and expansion, it produces structure that is fairly typical of all cases. Shown in Figure 2 are configurations of the evolved configuration at different times<fnref n="8" /> 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 $t \approx 3$, while the core is triaxial with a flatness parameter $t \approx 1$ and corresponds to a triaxial ellipsoid. We note that, once formed, the spiral-like 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 Figure 3).

Figure 4 shows the density profile $n(r)$, the velocity profile $v(r)$, and the energy profile $\epsilon(r)$ computed as averages in radial bins of constant logarithmic width. During the time evolution, the outer tail of $n(r)$ 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

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Configurations resulting from four different times (see the labels). The solid line in the upper left panel corresponds to the initial ellipsoid.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Configuration at $t = 25$: arrows are proportional to velocities.</figcaption></figure>

distribution, the tail of the density profile is well fit by a power-law behavior with $n(r) \sim r^{-4}$ (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 $\epsilon(r)$, particles are unbound (with $\epsilon > 0$), while in the core region particles are strongly bound (i.e., $\epsilon$ below -1.5); there is then an extended intermediate region in which many particles are marginally bound (i.e., $0 > \epsilon > -0.5$).

The energy distribution $P(\epsilon)$ at two different times ($t = 25, 50$), together with that of the initial conditions, is shown in Figure 5: we note that a large change of $P(\epsilon)$ has occurred during the gravitational collapse of the cloud at $\approx \tau_d$ while at later times the shape of the distribution remains approximately the same.

<footnote n="7">⁷ Here, and in the following figures, we use units of length in which $a_3 = 1$, and units of time in which $\tau_d = 1$; energies are given in units in which $GNm^2/a_3 = 1$.</footnote>

<footnote n="8">⁸ See *goo.gl/L1fRzZ* for the full movie of the time evolution.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4. From top to bottom: (i) density profile, (ii) velocity profile, and (iii) energy profile, at two different times: $t = 25$ (black dots) and at $t = 45$ (red dots).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. Energy distribution at $t = 0$ (green), at $t = 25$ (black) and at $t = 50$ (red).</figcaption></figure>

The upper panel of Figure 6 shows the average, in spherical shells of radius $r$, denoted $\langle \cdots \rangle$, of the radial component

$$v_r = \frac{\mathbf{r} \cdot \mathbf{v}}{r}$$

of the velocity, and of the “transverse” velocity

$$v_t = \frac{\mathbf{r} \times \mathbf{v}}{r},$$

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 $\langle |\mathbf{v}| \rangle = |\langle \mathbf{v}_r \rangle|$.

<figure><img src="image_9.png" /><figcaption>Figure 6. Configuration at $t = 45$. Upper panel: components of particle velocities averaged in spherical shells as a function of radius. Middle panel: anisotropy parameter $\beta(r)$. Lower panel: mass estimated from the velocity assuming stationary circular orbits, and the actual enclosed mass.</figcaption></figure>

The middle panel of Figure 6 shows the velocity anisotropy

$$\beta(r) = 1 - \frac{\langle v_r^2 \rangle}{2 \langle v_r^2 \rangle}.$$

Finally, the lower panel of Figure 6 shows $v^2r/G$, the mass that would be enclosed inside this radius if the motions were purely circular and the mass distribution spherically symmetric, and the mass $M(<r)$ actually enclosed inside the radius $r$.

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

$$\langle |\mathbf{v}_t| \rangle \gg |\langle \mathbf{v}_t \rangle|,$$

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

$$\langle |\mathbf{v}_r| \rangle \approx |\langle \mathbf{v}_r \rangle| \gg \langle |v_r| \rangle;$$

correspondingly (lower panel of Figure 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.,

$$\langle |\mathbf{v}| \rangle \approx \langle |v_r| \rangle \approx |\langle v_r \rangle|.$$

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 Figure 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 $\langle|v_r|\rangle \approx 0$ 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 Figure 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 the 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 $b_{rot}$). 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 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 $200 \text{ km s}^{-1}$, we define the dimensionless parameters as follows:

$$v_{200} = \frac{a_3 / \tau_d}{200} \text{ km s}^{-1}.$$

We can then write

$$a_3 \approx \left( \frac{200 v_{200}}{n} \times t_{\text{Gyr}} \right) \text{kpc},$$

where n is the number of dynamical times in our simulation and $t_{\text{Gyr}}$ is its duration given in billion of years. Thus, for $n \approx 50$, as in Figure 2, and taking $t_{\text{Gyr}} \sim 1$, which corresponds to a mass (by using Equation (1))

$$M = \frac{\pi^2 a_3^3}{8G\tau_d} \approx 10^{11} M_{\odot},$$

we have that region R1 extends to $\sim 2$ kpc, and region R2 extends to $\sim 50$ kpc. 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 that

generated the disk and arms occurred much more recently than the formation of the oldest stars in these galaxies (with an age $\sim 10$ 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; 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 the compatibility of the large-scale motions of our generated mass distributions with the 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 15 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., $> 15$ 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—along which there is maximal variation of the projected velocities—with the projected 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 30% 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 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 $v_{los}$ of random observers, characterized by two angles: the inclination angle $i$, defined as the angle between the vector $u_o$ giving the orientation of the observer's LOS and the vector $u_g$ in the direction of the shortest semi-principal axis of the model galaxy, and an angle $j$ 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, 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 Figure 7 are the projected velocity maps for the same

simulation analyzed above, for an observer with $i = 30^\circ$ and $j = 30^\circ$. The maps have been averaged on a grid of size $64^2$ (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 $40^\circ$;

(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 using the one-dimensional LOS velocity profile along the estimated kinematic axis);

(vi) the two-dimensional LOS velocity 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., using Equation (7)).

We then explored (see Figure 9) the full range of $i$ and $j$. Only for $j$ very close to $\pi/2$ (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 10%–30%, except for $j \to 90$, in which they can be as high as $\approx 50\%–70\%$. In these images one can discern clearly that our model galaxies are non-axisymmetric at larger radii, and as we have noted, there is a strong correlation between the direction of the outer radial velocities with the intrinsic longest semi-principal 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 longest semi-principal axis is typically very close to the major axis of the projected image, and the large radial velocities project out their components 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.

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
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Projection for $i = 30^{\circ}$, $j = 30^{\circ}$ (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.</figcaption></figure>

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 Figure 6 and of Figure 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, the observed phenomenological relation, $a_c \propto \sqrt{g(r)}$, where $a_c$ is the centripetal acceleration inferred from the apparent motions, and $g(r)$ is the gravitational acceleration of the visible baryonic mass (see e.g., McGaugh et al. 2016), 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 $a_c \approx \frac{v_{\max}^2}{r}$ where $v_{\max}$ is the inferred constant velocity of rotation, while $g(r) \approx GM_c/r^2$, where $M_c$ is the mass in the virialized core. Thus,

$$a_c \approx \frac{v_{\max}^2}{r} \approx \sqrt{a_0 g(r)}$$

where

$$a_0 = \frac{v_{\max}^4}{GM_c}$$

The observed approximate constancy of $a_0$ for different systems then corresponds to $v_{\max}^4 \propto M_c$, i.e., the Tully-Fisher relation.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Projection for *i* = 30°, *j* = 30°. 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.</figcaption></figure>

## 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 (D. Benhaiem et al. 2017, in preparation), 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 velocity fields that 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 that 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 detail 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 timescale for collapse, as we have discussed, must be assumed to be 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 author of this work were granted access to the HPC resources of The Institute for Scientific Computing and Simulation financed by Région 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.

## Appendix

We detail here how we construct the projected velocity maps reported in Section 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 $\lambda_i$, where $\lambda_1 \le \lambda_2 \le \lambda_3$, and corresponding eigenvectors $\lambda_i$. The longest semi-principal axis is then designated by a unit vector $\mathbf{u}_1$ parallel to $\lambda_1$, the intermediate semi-principal axis by a unit vector $\mathbf{u}_2$ parallel to $\lambda_2$, and the shortest semi-principal axis by $\lambda_3$. The plane of the galaxy is then orthogonal to $\lambda_3$. We then
<figure><img src="image_12.png" /><figcaption>Figure 9. Residuals for i = 60° and different values of j from 0° to 80° in steps of 10°.</figcaption></figure>

rotate from our original Cartesian axes (of the simulation) to determine the components of the particle positions, $x_i$, and their velocities, $v_i$, in the new basis {$u_i$}.

## A.2. Orientation of the Observer

Following standard conventions (see, e.g., Beckman et al. 2004) we define the inclination angle $i$ of the observer as the angle between his LOS and a vector orthogonal to the plane of the galaxy, which we take to be $u_3$. Furthermore, as the galaxy is non-axisymmetric about this axis, we define an azimuthal angle $j$ 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

$$u_o = \sin(i)\cos(j)u_1 + \sin(i)\sin(j)u_2 + \cos(i)u_3, \quad (2)$$

## 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

$$\begin{aligned} u_x &= \sin(j)u_1 - \cos(j)u_2, \\ u_y &= \cos(j)u_1 + \sin(j)u_2 \end{aligned} \quad (3)$$

in the plane of the galaxy. The vector $u_y$ is thus parallel to the axis of the projection in the plane of the galaxy of the observer LOS, while $u_x$ 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

$$u_x'' = u_x, \quad u_y' = u_o \times u_x. \quad (4)$$

Using the expressions above, a little algebra gives

$$\begin{aligned} u_1 &= \sin(j)u_x' + \cos(j)\cos(i)u_y' + \cos(j)\sin(i)u_o \\ u_2 &= -\cos(j)u_x' + \sin(j)\cos(i)u_y' + \sin(j)\sin(i)u_o \\ u_3 &= -\sin(i)u_y' + \cos(i)u_o. \end{aligned} \quad (5)$$

The position coordinates ($x', y'$) of the particles in the plane of projection, and projected velocity $v_{los} = v \cdot u_0$, can then be calculated, for any given observer ($i, j$), as

$$\begin{aligned} x' &= x_1 \sin(j) - x_2 \cos(j) \\ y' &= x_1 \cos(i)\cos(j) + x_2 \cos(i)\sin(j) - x_3 \sin(i) \\ v_{los} &= v_1 \sin(i)\cos(j) + v_2 \sin(i)\sin(j) + v_3 \cos(i). \end{aligned} \quad (6)$$
### 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 spectrographic 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 $\Delta$ which is a small fraction of the minor axis. From these LOS velocity profiles along the kinematic axis $v_{los}(R)$, we estimate the mass $M_c(R)$ enclosed in the radius $R$ assuming that particles are in stationary circular orbits as

$$M_c(R) = \frac{v_{los}^2(R)R}{\sin(i)G}, \quad (7)$$

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

$$v_{los}(r, \phi) = v_{\theta} \sin(i)\cos(\theta) + v_{R} \sin(i)\sin(\theta), \quad (8)$$

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

$$\tan(\theta) = \tan(\phi) / \cos(i) \\ R = r \cos(\phi) / \cos(\theta). \quad (9)$$

For a purely rotating axisymmetric model, $v_R = 0$ and $v_\theta = v_\phi(R)$. The kinematic axis is that along which there is maximal variation of the projected velocity, i.e., $\theta = \phi = 0$.

### 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 axis, assumed to correspond to 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 $i$: we do this by minimizing the

residuals between the rotational model, computed for a generic $i$, and the actual data on each grid cell. To do so we compute first, for each grid cell, labeled by $\alpha$ and centered on projected coordinates $x'_\alpha, y'_\alpha$, the polar coordinates as defined above:

$$r_{\alpha} = \sqrt{(x'_{\alpha})^2 + (y'_{\alpha})^2} \\ \phi_{\alpha} = \arccos(x'_{\alpha}/r_{\alpha}) \\ R_{\alpha} = r_{\alpha} \sqrt{\cos(\phi'_{\alpha})^2 + \sin(\phi'_{\alpha})^2 / \cos(i)^2} \\ \theta_{\alpha} = \arctan(\tan(\phi'_{\alpha}) / \cos(i)). \quad (10)$$

Then, for the given value of the inclination angle $i$, we use Equation (8) (with $v_R = 0$) to compute the LOS velocity of the rotational model, denoted $v_{los, model}^\alpha$. 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 $v_{los}(R)$ 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$, i.e.,

$$\text{Residuals} = \sum_{\alpha} |v_{\text{los}}^{\alpha} - v_{\text{los, model}}^{\alpha}|. \quad (11)$$

## ORCID iDs

Francesco Sylos Labini https://orcid.org/0000-0003-2236-4537

## References

Aarseth, S., Lin, D., & Papaloizou, J. 1988, ApJ, 324, 288

Aguilar, J., & Merritt, D. 1990, ApJ, 354, 33

Barnes, E. I., Lanzel, P. A., & Williams, L. L. R. 2009, ApJ, 704, 372

Beckman, J. E., Zurita, A., & Vega Beltrán, J. C. 2004, LNEA, 1, 43

Benhaiem, D., & Sylos Labini, F. 2015, MNRAS, 448, 2634

Benhaiem, D., & Sylos Labini, F. 2017, A&A, 598, A95

Boily, C., Athanassoula, E., & Kroupa, P. 2002, MNRAS, 332, 971

Brown, W. 2015, ARA&A, 52, 15

Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJL, 785, L12

Combes, F. 2002, NewAR, 46, 755

David, M., & Theuns, T. 1989, MNRAS, 240, 957

Dobbs, C., & Baba, J. 2014, PASA, 31, e035

Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748

Elmegreen, B. G., Seiden, P. E., & Elmegreen, D. M. 1989, ApJ, 343, 602

Erroz-Ferrer, S., Knapen, J. H., Font, J., & Beckman, J. E. 2015, HiA, 16, 328

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Hénon, M. 1973, A&A, 24, 229

Henry, A. L., Quillen, A. C., & Guttermuth, R. 2003, AJ, 126, 2831

Joyce, M., Marcos, B., & Sylos Labini, F. 2009, MNRAS, 397, 775

Kalberla, P. M. W., & Dedes, L. 2008, A&A, 487, 951

López-Corredoira, M., & González-Fernández, C. 2016, AJ, 151, 165

McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, PhRvL, 117, 201101

Milgrom, M. 1983, ApJ, 270, 365

Milgrom, M. 2016, arXiv:1609.06642

Rubin, V. C. 1983, Sci, 220, 1339

Sofue, Y. 2017, PASJ, 69, R1

Sofue, Y., & Rubin, V. 2001, ARA&A, 39, 137

Springel, V., Yoshida, N., & White, S. D. M. 2001, NewA, 6, 79

Sylos Labini, F. 2012, MNRAS, 423, 1610

Sylos Labini, F. 2013, MNRAS, 429, 679

Theuns, T., & David, M. 1990, Ap&SS, 170, 276

Tully, R. B., & Verheijen, M. A. W. 1997, ApJ, 484, 145

van Albada, T. 1982, MNRAS, 201, 939

Vogel, S. N., Rand, R. J., Gruendl, R. A., & Teuben, P. J. 1993, PASP, 105, 666

Ward, W. R., & Canup, R. M. 2010, AJ, 140, 5