Galactic dynamics
| George Contopoulos and Christos Efthymiopoulos (2011), Scholarpedia, 6(5):10670. | doi:10.4249/scholarpedia.10670 | revision #91294 [link to/cite this article] |
Galactic dynamics is the study of the motions of the stars, gas and dark matter in order to explain the main morphological and kinematical features of galaxies.
In the present article we focus on topics of galactic dynamics which are most relevant to the general theory and methods of dynamical astronomy (see also scholarpedia section "Extragalactic Astronomy").
Contents |
Introduction and basic concepts
The main components of the Universe are the galaxies, which are composed of billions of stars, but also of gas, dust and dark matter. The most well known classification of galaxies is due to Hubble. The main types are elliptical galaxies (E), normal spiral galaxies (S), barred spiral galaxies (SB) and irregular galaxies (I)( Figure 1 ).
The elliptical galaxies have various degrees of ellipticity from zero (E0) to \(0.7\) (E7), and they are slowly rotating. On the other hand the spiral and barred spiral galaxies are rotating fast and they contain spiral arms. There are tight spirals (Sa, SBa), intermediate (Sb, SBb) and open spirals (Sc, SBc). They are flat systems, containing also gas and dust, out of which new stars are formed continuously. The spiral and barred spiral galaxies are also called `disc galaxies', since the spiral arms are embedded in a thin disc (the thickness is of order less than 10%). Finally, the irregular galaxies are relatively small systems that accompany large spiral galaxies.
The shapes of the galaxies are governed mainly by the gravitational interactions between their individual mass components, i.e. the stars, dark matter, and gas. The stars form the main body of the luminous part of a galaxy. The gas forms a relatively small proportion of matter (up to 10% in disc galaxies). The stars and the gas together are called `baryonic matter'. The dark matter, on the other hand, contains more mass than the stars and gas, and it extends to large distances (one order of magnitude larger than the baryonic matter).
The motions of the stars and of the dark matter elements are governed purely by their gravitational forces. The study of these motions and their combinations to form self-consistent statistical mechanical configurations, is called stellar dynamics. This constitutes the central approach to galactic dynamics. Gas dynamics, on the other hand, is governed also by dissipative forces due to pressure, radiation, magnetic fields etc. The latter also influence, to some extent, particular morphological and kinematical features of the galaxies.
The gravitational potential of a galaxy is composed of a mean field, due to the average distribution of the galactic matter, and fluctuations due to the approaches (encounters) between individual stars. An estimate of the effects of these encounters is provided by the "relaxation time" of the system (Chandrasakhar), which can be estimated by the time required by a star to change its average direction of motion by 90° due to the cumulative effects of encounters. This time is very long, of the order of 1012 years, while the periods of the motions of the stars around the center of the galaxy are of the order of 108 years ("dynamical time"). Thus, in a first approximation, we may consider the orbits of the stars as due to the general distribution of the galactic matter. The same applies to the orbits of dark matter elements. In a better approximation, however, various details of the relaxation process have to be taken into account. This issue is discussed in section 6 (N-body systems).
If \(V(\textbf{x},t)\) is the average galactic potential the equations of the average motion are
\[\tag{1} \ddot{\mathbf{x}}=-\frac{\partial V}{\partial\mathbf{x}} \]
where \(\textbf{x}=(x,y,z)\ .\)
The distribution function \(f\) is the density of matter in phase space \((\textbf{x},\textbf{v})=(x,y,z,\dot{x},\dot{y},\dot{z})\ .\) It is given by the "collisionless Boltzmann equation" \[\tag{2} \frac{\partial f}{\partial t} + \frac{\partial f}{\partial \textbf{x}}\textbf{v}- \frac{\partial f}{\partial \textbf{v}} \frac{\partial V}{\partial \textbf{x}}=0 \]
The total density \(\rho(\textbf{x},t)\) is found by integrating
\(f\)
over all velocities
\[\tag{3}
\rho(\textbf{x},t)=\int f(\textbf{x},\textbf{v},t)d^3\textbf{v}
\]
while the average potential \(V(\textbf{x},t)\) is given by Poisson's equation \[\tag{4} \nabla^2V=4\pi G\rho~~ \]
In the treatment of dark matter we may consider i) rigid halo models with a fixed dark matter density \(\rho_d\ ,\) where \(f\) and \(\rho\) refer to the stellar distribution function and density while (4) is replaced by \(\nabla^2V=4\pi G(\rho+\rho_d)\ ,\) or ii) live halo models in which the dark matter is responsive to, and affects the collective motions of stars. Rigid halo models have been used extensively in orbital studies of stars (or clusters). However, the rigid halo approach is not convenient in circumstances where the halo interacts collectively with a particular stellar sub-structure like a bar (see section 6). N-body simulations have shown that in such cases only a `live' halo model can capture correctly the effects of such interactions.
By solving the equations of motion in a given potential \(V(\textbf{x},t)\) we define smooth orbits of the stars or dark matter elements. The superposition of many orbits with appropriate weights yields a response density. The system must be self -consistent (self-gravitating) i.e. the response density must be equal to the imposed density.
The self-consistency condition is satisfied accurately in N-body simulations (section 6). However in many cases we assume a fixed potential that represents a model galaxy. In this case the self-consistency condition can be checked a posteriori, i.e. after the orbits are calculated in the fixed potential.
We frequently consider stationary models, i.e. \(V=V(\textbf{x})\) is considered independent of time. In N-body simulations, we often consider successive snapshots and study the forms of the orbits at these times.
The orbits in a given galaxy can be ordered (periodic or quasiperiodic) or chaotic. This classification is based on the number of integrals of motion obeyed by an orbit. An integral of motion is a function of the canonical variables (position and momenta) that remains constant along an orbit. The importance of integrals in galactic dynamics stems from Jeans' theorem which states that a stationary distribution function \(f(\mathbf{x},\mathbf{v})\) can depend on its arguments only through the integrals \(I_i\ ,\) i.e. \(f\equiv f(I_1(\mathbf{x},\mathbf{v}), I_2(\mathbf{x},\mathbf{v}),...)\ .\)
In stationary galactic systems of \(n\) degrees of freedom there can be at most \(n\) independent integrals which appear as arguments of \(f\ .\) The energy \(E\) is an obvious integral in all stationary systems, while the angular momentum \(L\) along the axis of symmetry is a second integral in axisymmetric systems. Additional `third integrals', of exact or approximative form, can be found in various cases and they play a key role in dynamics.
A main property of galactic kinematics is the form of the velocity ellipsoids. A velocity ellipsoid is defined at any point of space \(\mathbf{x}\) by the three principal axes arising by the diagonalization of the \(3\times 3\) velocity dispersion tensor \(\mathbf{\sigma}\) with elements \(\sigma_{ij}\) where \[\tag{5} \sigma_{ij}^2(\textbf{x})={1\over\rho(\textbf{x})} \int (v_i-V_i)(v_j-V_j)f(\textbf{x},\textbf{v})d^3\textbf{v} \]
where \(V_i\) and \(V_j\) are the mean velocities along the axes \(i\ ,\) \(j\) of an orthogonal coordinate system (\(i\ ,\) \(j\) run from 1 to 3).
If the distribution function depends only on the energy, \(f=f(E)\ ,\) the velocity ellipsoid at any point of space is a sphere, i.e. the velocity dispersion is the same in all directions. If, on the other hand, \(f\) depends both on \(E\) and \(L\ ,\) i.e. \(f=f(E,L)\ ,\) the velocity ellipsoid is a spheroid (two equal axes). Finally, if \(f\) depends on three integrals, i.e. \(f\equiv f(E,L,I_3)\ ,\) all three axes of the velocity ellipsoid are unequal. This distinction is important, because it allows one to link the kinematic observations available for a particular system to dynamical features of the same system. In fact, near the Sun the velocity ellipsoid has three unequal axes.
An important theorem in galactic dynamics is the virial theorem. This theorem states that in a stellar system in equilibrium, the following equation holds: \[\tag{6} 2T_{ij}+\Pi_{ij}+W_{ij}=0 \]
where \(T_{ij}={1\over 2}\int \rho(\textbf{x})V_iV_jd^3\textbf{x}\) and \(\Pi_{ij}=\int\rho(\textbf{x})\sigma^2_{ij}d^3\textbf{x}\ .\) The quantities \(K_{ij}=T_{ij}+\Pi_{ij}/2\) and \(\Pi_{ij}\) are called kinetic energy tensor and pressure tensor respectively, while \(W_{ij}=-{G\over 2}\int\int \rho(\textbf{x},\textbf{x}') {x_i(x_j-x_j') \over|\textbf{x}-\textbf{x}'|} d^3\textbf{x}d^3\textbf{x}'\) is the potential energy tensor. The scalar virial theorem, obtained by taking the trace of Eq.(6) reads \[\tag{7} 2K+W=0 \]
where \(K\) and \(W\) are the total kinetic and potential energies of the stellar system in equilibrium. Depending on the shape of a system, different states of virial equilibrium may exist with different degrees of kinetic energy that go to rotation or to velocity anisotropy. Accordingly, the elliptical galaxies as well as the spheroidal bulges of disc galaxies are divided into those being rotationally supported or pressure supported.
Orbits and Integrals
The study of individual orbits in fixed potential models with different degrees of symmetry is a central subject of galactic dynamics, since such models offer idealizations of the true gravitational potential of galaxies. A generic feature of such models is the co-existence of ordered and chaotic orbits. In fact, the properties of ordered orbits can be unraveled using various forms of the canonical perturbation theory. Such is the theory of the third integral as well as the Kolmogorov-Arnold-Moser (KAM) theory. On the other hand, the properties of chaotic orbits can be explored mainly by numerical means such as the Poincaré surface of section or quantities such as the Lyapunov characteristic number. Readers are deferred to (Contopoulos 2004) for an instructive introduction to these topics. Below, we review only the most basic facts relevant to the orbits in galaxies. We examine first the orbits and forms of integrals in simple (non-rotating) systems like elliptical galaxies. Then we examine the orbits and integrals in fast rotating disc galaxies. We return to the issue of self-consistency in sections 5 (spiral structure), and in section 6 (N-body simulations).
Orbits in axisymmetric galaxies
Within the approximation of collisionless stellar dynamics, the orbits in an axisymmetric galaxy are governed by a smooth time-independent axisymmetric gravitational potential \(V(R,z)\ ,\) which corresponds to the solution of Poisson's equation for an axisymmetric matter distribution \(\rho(R,z)\ .\) Then, the orbit of a star (point particle) becomes independent of its mass, and it is given by a Hamiltonian of the form \[\tag{8} H\equiv{p_R^2\over 2}+{p_\vartheta^2\over 2R^2}+{p_z^2\over 2}+V(R,z)=E \]
where \(R,\vartheta,z\) are cylindrical coordinates and \(p_R\ ,\) \(p_\theta\ ,\) \(p_z\) are the corresponding canonical momenta.
The orbits in the equatorial plane are rosettes ( Figure 2 ). Their angular momentum is constant \(p_\theta=J_0\ .\) The radial oscillations are given by \[\tag{9} \ddot{R}=-V'_0(R)+\frac{J^2_0}{R^3} \]
where \(V_0(R)=V(R,z\!=\!0)\ .\) For given energy and angular momentum the radius and angular velocity of the circular orbit \(R_0\) are given by \[\tag{10} J^2_0=R_0^3V'_0,~~~\Omega(R_0)=\sqrt{V'_0/R_0}~~. \]
For orbits close to circular, the frequency of radial oscillation is called epicyclic frequency, given by \[\tag{11} \kappa(R_0)=\left(V''_0+\frac{3V'_0}{R_0}\right)^{1/2} \]
All the orbits on the equatorial plane are ordered. On the other
hand, orbits with vertical oscillations with respect to the
equatorial plane can be ordered or chaotic.
Ordered orbits are found by calculating a `third integral' of motion in the form of a series (Contopoulos). Expanding the Hamiltonian with respect to the radius \(R_0\) of a circular orbit, the Hamiltonian describing the motion in a meridian plane takes the form \[\tag{12} H=\frac{1}{2}(\dot{x}^2+\dot{y}^2+\omega^2_1~x^2+\omega^2_2y^2)+ \varepsilon xy^2+\varepsilon'x^3+\mbox{higher order terms} \]
where \(x\) is the difference of the radial coordinate \(R\) from the reference radius \(R_0\ ,\) while \(y\) is parallel to the axis of symmetry \(z\ .\)
By definition, a formal integral of motion \(\Phi(x,y,\dot{x},\dot{y})\) is a function of the phase space coordinates which has a vanishing Poisson bracket with the Hamiltonian function, namely \[\tag{13} [\Phi,H]=\frac{\partial\Phi}{\partial x}\frac{\partial H}{\partial\dot{x}}+\frac{\partial\Phi}{\partial y}\frac{\partial H}{\partial\dot{y}}-\frac{\partial\Phi}{\partial\dot{x}}\frac{\partial H}{\partial x}-\frac{\partial\Phi}{\partial\dot{y}}\frac{\partial H}{\partial y}= 0 \]
If we develop \(\Phi\) and H in power series \(\Phi=\Phi_2+\Phi_3+...,~H=H_2+H_3+...\ ,\) we can find the successive terms of \(\Phi\) by equations of the form \[\tag{14} [\Phi_N,H_2]=-[\Phi_{N-1},H_3]-...-[\Phi_2,H_N] \]
provided that \(\Phi_2\) is a properly chosen function. The successive terms in the power series are of increasing order in a properly chosen small parameter. The latter gives the order of the distance of orbits in phase space from the equatorial circular orbit. More generally, in `third integral' expansions the small parameter can be chosen so as to reflect the deviations of the model considered from an integrable model.
A non-resonant third integral can be found if \(\omega_1\) and \(\omega_2\) satisfy no commensurability condition. We then choose, for example, \(\Phi_2=\frac{1}{2}(\dot{x}^2+\omega^2_1~x^2)\ .\) The higher order terms \(\Phi_3,\ldots,\Phi_N\) can be calculated by computer algebra.
The third integral is in general a non-convergent series. However a proper truncation provides useful approximations for the forms of galactic orbits. In fact, it has been found that a truncated third integral \(\overline{\Phi}_N=\Phi_2+...+\Phi_N\) is approximately conserved with greater accuracy as the order N increases beyond N=2. However if the expansion goes beyond a critical value \(N_{crit}\) the variations of \(\overline{\Phi}_N\) increase instead of decreasing. The order of \(N_{crit}\) is smaller when the perturbations \(\varepsilon,\varepsilon'\) etc are larger. By Nekhoroshev's theorem, the error of the approximation is exponentially small in \(1/\epsilon\) (or \(1/\epsilon'\)) at the order \(N_{crit}\ .\)
In calculating the various terms of the third integral we find that \(\Phi_N\) contains denominators of the form \(|\omega_1k_1+\omega_2k_2|\) where \(k_1\) and \(k_2\) are positive or negative integers. Thus if \( \left|\frac{\omega_1}{\omega_2}\right|\) is close to a rational number \( \left|\frac{k_2}{k_1} \right|\ ,\) the quantity \(|\omega_1k_1+\omega_2k_2|\) may be small. Then we say that we have a small divisor. In such cases we can find a resonant form of the third integral that is valid near the particular resonance \(\left|\frac{\omega_1}{\omega_2}\right|= \left|\frac{k_2}{k_1}\right|\ .\) In particular we find the corresponding resonant periodic orbits and the tube orbits that surround the stable resonant periodic orbits. Of special importance are the low order resonances \(\omega_1/\omega_2=\pm 1/1,~\pm 2/1\ ,\) etc.
The ordered orbits obeying a non-resonant third integral are called `tube' orbits if \(J_0\neq 0\ ,\) and `box' orbits if \(J_0=0\ .\) Box orbits pass arbitrarily close to the center, while tube orbits leave a hole around the axis of symmetry. On the other hand, the orbits obeying a resonant third integral are either periodic or they form thin tubes around their corresponding periodic orbits.
In the case of a galaxy with a smooth core the orbits near the center are ordered. The boxes are deformed Lissajous figures ( Figure 3 a ), while the orbits around particular resonant periodic orbits form elongated tubes ( Figure 3 b ). But there also many chaotic orbits ( Figure 3 c ), mainly in the region separating the box orbits from the main tube orbits.
If now the core of a galaxy is cuspy (e.g. the density rises as a power law as we approach the center), or if we add a central mass (e.g. a black hole) at the center of the galaxy, all the orbits near the center become chaotic. At the same time the number of tube orbits increases. In order to distinguish between ordered and chaotic orbits we use a Poincaré surface of section, i.e. a surface in phase space that intersects all the orbits, and find the distribution of the successive intersections (iterates) of each orbit by this section. In the case of orbits in the meridian plane \((R,z)\) of a galaxy with a plane of symmetry \(z=0\ ,\) the plane \(z=0\) is a Poincaré surface of section ( Figure 4 ). The energy \[\tag{15} E=\frac{1}{2}(\dot{R}^2+\dot{z}^2)+V(R,z) \]
is then one integral of motion. If we also use the `third integral' \[\tag{16} \Phi(\dot{R},\dot{z},R,z)=c \]
we eliminate \(\dot{z}\) between (15) and (16) and find a toroidal surface \[\tag{17} F(R,z,\dot{R})=q \]
on which lies the orbit. The successive intersections of an orbit by the plane \(z=0\) lie on an invariant curve \[\tag{18} F(R,0,\dot{R})=q \]
The ordered orbits thus define invariant curves of the form of (18) . On the other hand the successive iterates of a chaotic orbit are scattered irregularly ( Figure 4 ).
In a generic dynamical system the ordered and chaotic orbits coexist. In fact, if an integrable system is perturbed slightly, there is a large set of integral surfaces containing regular orbits. This result is based on the famous Kolmogorov, Arnold, Moser (KAM) theorem that can be stated as follows. If an autonomous Hamiltonian system of N degrees of freedom is close to an integrable system expressed in action-angle variables, there are N-dimensional invariant tori containing quasi-periodic motions with frequencies \(\omega_j\) satisfying a Diophantine condition \[\tag{19} |\mathbf{\omega\cdot k}|=|\omega_1 k_1+\omega_2k_2+... +\omega_Nk_N|>\frac{\gamma}{[|k_1|+|k_2|+...|k_N|]^\tau} \]
where \(\tau>N-1\ .\) The set of invariant tori is of order \([1-O(\gamma)]\ .\) thus if the perturbation \(\gamma\) is small most initial conditions generate ordered orbits. However near every resonance where \(|\mathbf{\omega\cdot k}|\) is small there are sets of chaotic orbits. On the other hand if \(\gamma\) is large a large degree of chaos appears. The transition from order to chaos is produced by an overlapping of resonances (Rosenbluth,Sagdeev, Taylor & Zaslavski, Contopoulos, Chirikov). In fact at every resonance of a nonintegrable system there is one or two unstable periodic orbits and close to every unstable orbit there is some degree of chaos ( Figure 5 a ). If the perturbation is small the chaotic zones of the various resonances are separated by invariant curves around the center. But if the perturbation increases the intervening invariant curves are destroyed and the various chaotic domains overlap and produce a large degree of chaos ( Figure 5 b ).
A general method to distinguish between ordered and chaotic orbits is by means of the Lyapunov characteristic number (LCN) defined as follows. We calculate two orbits \(x(t,x_0)\) and \(x(t,x_0+\xi_0)=x(t,x_0)+\xi\ ,\) where \(\xi\) is an infinitesimal deviation, found by solving the variational equations or by approximate numerical methods. The Lyapunov characteristic number is \[\tag{20} LCN=\lim\limits_{ t\rightarrow\infty}\frac{ln|\xi/\xi_0|}{t} \]
where \(\xi_0\) and \(\xi\) are the deviations at times \(0\) and \(t\ .\) If \(LCN=0\) the orbit is ordered and if \(LCN>0\) the orbit is chaotic. This method is applicable for any number of degrees of freedom. The inverse of the Lyapunov characteristic number is called "Lyapunov time". Beyond the Lyapunov time the orbits are unpredictable.
There are several practical methods to improve the method of the Lyapunov characteristic number, like the fast Lyapunov indicator (FLI, Froeschlé), the stretching numbers and helicity angles (Voglis), the mean exponential growth of nearby orbits (MEGNO, Cincotta, Simo) and the smaller alignment index (SALI, Skokos). On the other hand, a global analysis of phase space dynamics can be done very efficiently with the frequency analysis of Laskar.
Orbits in triaxial galaxies
In the study of triaxial galaxies, a basic starting model is the integrable model given by the Staeckel potential \[\tag{21} V=-\frac{F_1(\lambda)}{(\lambda-\mu)(\lambda-\nu)} -\frac{F_2(\mu)}{(\mu-\nu)(\mu-\lambda)} -\frac{F_3(\nu)}{(\nu-\lambda)(\nu-\mu)} \]
in elliptical coordinates \(\lambda,\mu,\nu\ .\) The most important types of orbits are i) box, (ii) and (iii) inner and outer long axis tube, and (iv) short axis tube orbits ( Figure 6 ).
The ILAT, OLAT and SAT orbits characterize the form of regular quasi-periodic orbits in generic triaxial potentials. On the other hand, the box orbits are present only in models of triaxial galaxies with a smooth core, i.e. one characterized by a flat central density profile. If, instead, the central density profile is cuspy, most box orbits disappear and are replaced by chaotic orbits which may pass arbitrarily close to the center.
If we now consider periodic galactic orbits in 3 dimensions we have a number of new phenomena.
The type of stability or instability of the periodic orbits depends on the eigenvalues of the monodromy matrix \(A\ ,\) that gives the infinitesimal deviations from the periodic orbit after one period \(T\ :\) \[\tag{22} \xi(T)=A\xi_0 \]
The eigenvalues \(\lambda\) of A satisfy an equation of the form \[\tag{23} (\lambda^2+b_1\lambda+1)(\lambda^2+b_2\lambda+1)=0 \]
Thus there are two couples of inverse eigenvalues \[\tag{24} \begin{array}{c} \lambda _1 \\ \lambda _2 \end{array} = {1 \over 2} (-b_1\pm \sqrt{ b^2_1-4}),\begin{array}{c} \lambda _3 \\ \lambda _4 \end{array} = {1 \over 2} (-b_2\pm \sqrt{ b^2_2-4}) \]
If \(|b_1|<2\ ,\) \(|b_2|<2\) the orbit is stable (S), if \(|b_1|<2\ ,\) \(|b_2|>2\ ,\) or \(|b_1|>2\ ,\) \(|b_2|<2\) the orbit is simply unstable (U), if \(|b_1|>2\ ,\) \(|b_2|>2\) the orbit is doubly unstable (DU), and if \(b_1\ ,\) \(b_2\) are complex, the orbit is complex unstable (\(\Delta\)). In a 2-D system there is only one factor of (23) and the orbits are either stable or unstable.
As an example of a 3-D system we consider the Hamiltonian \[\tag{25} H=\frac{1}{2} (\dot{x}^2+\dot{y}^2+\dot{z}^2+Ax^2+By^2+C z^2)-\varepsilon xz^2-\eta yz^2=h \]
for fixed values of \(h,A,B,C\) and varying parameters \(\varepsilon\) and \(\eta\ .\) A simple orbit, called 1a, has various types of stability, as \(\varepsilon\) and \(\eta\) vary. Figure 7 is a 'stability diagram' that contains all four types of stability-instability.
In the case of 3-D galaxies the 3-D orbits undergo an oscillation along the third dimension, whatever the form of their projection on the galactic plane. Thus we may have instability along the third dimension even if the projection of an orbit on the plane of symmetry is stable.
Another new phenomenon that appears in three or more degrees of freedom is Arnold diffusion. This is a slow diffusion that allows orbits to go very far from their initial conditions by following various chaotic layers in phase space. In the case of 2 degrees of freedom the phase space for a fixed energy is 3 dimensional. Then if there is a 2-dimensional invariant toroidal surface (and many such surfaces exist in slightly perturbed systems) the orbits inside the torus cannot cross the torus and go outside. Thus if the inner orbits are chaotic they cannot diffuse very far.
However in a system of 3 degrees of freedom the phase space (for constant energy) is 5-dimensional, but the invariant surfaces provided by the KAM theorem are 3-dimensional and cannot separate the phase space into an interior and an exterior part. Thus diffusion can carry an orbit very far.
A picture of such a diffusion is provided if we consider a 3-D phase space and 1-D invariant manifolds, like strings from the floor to the ceiling of a room. If the system is non-integrable these strings leave between them chaotic domains, and the diffusion along these domains is Arnold diffusion.
Arnold diffusion can change considerably the galactic orbits, however the time scale of this diffusion is very long and it cannot affect appreciably galaxies that are not very irregular. A consequence is that despite Arnold diffusion the galaxies do not change appreciably over a Hubble time.
Orbits in disc galaxies
The patterns of disc galaxies (e.g. bars or spiral arms) rotate with an angular velocity \(\Omega_s\) which is large compared to the slow figure rotation of elliptical galaxies. Thus, in order to study the orbits in disc galaxies it is customary to choose a rotating frame of reference.
In an axisymmetric disc, the equations of motion in polar coordinates \((R,\vartheta)\) in the rotating frame are: Eq. (6) for the radial coordinate and \[\tag{26} \dot{\vartheta}+\Omega_s=J_0/R^2 \]
for the angular coordinate.
The radius of the circular orbit \(R_0\) as well as the angular and epicyclic frequencies \(\Omega,\kappa\) are defined by (10) , (11) . The combination \[\tag{27} H=E_0-\Omega_sJ_0=h \]
is called Jacobi constant (or energy in the rotating frame).
From a given potential \(V_0(R)\) we find three main functions (curves) ( Figure 8 ), namely \(\Omega\ ,\) and \(\Omega\pm\)\(\kappa/2\ .\) The intersections of these curves by the line \(\Omega_s\) are corotation, where \(\Omega=\Omega_s\ ,\) and outer or inner Lindblad resonances, where \[\tag{28} \frac{\kappa}{\Omega-\Omega_s}=\mp\frac{2}{1} \]
In generic potentials we have one corotation, one outer Lindbald resonance (OLR) and one, zero, or two inner Lindblad resonances (ILR).
A convenient approach is to analyze the gravitational potential in terms of Fourier components \[\tag{29} V(R,\theta,z)= \sum_{m=0}^{\infty}V_m(R,z)\cos\big[m\theta-\phi_m(r)\big] \]
In most `grand design' spiral or barred galaxies, the m=2 component dominates over all other components except \(V_0\) (axisymmetric component). The main families of orbits in such galaxies are examined below. The remaining components \(V_m\) (\(m\neq\)0 or 2) describe deformations of the form of the spiral or bar structures from a purely bi-symmetric shape. It has been found observationally that these components may play a significant dynamical role in particular galaxies.
The periodic orbits at the Lindblad resonances are approximately ellipses \((2/1)\) around the center of the galaxy ( Figure 9 a ), but there are also higher order resonances, e.g. triple \((3/1)\ ,\) quadruple \((4/1)\) ( Figure 9 a ), etc. The regular non-periodic orbits are close to periodic orbits. They have two basic frequencies, one nearly equal to that of the corresponding periodic orbit, and one representing oscillations around this periodic orbit. Such orbits are called quasiperiodic ( Figure 9 b,c ).
The main families of periodic orbits inside corotation are the families \(x_1\) (stable, Figure 10 a), \(x_2\) (stable, oriented perpendicularly to \(x_1\ ,\) ( Figure 10 a) and \(x_3\) (unstable).
We call "characteristic" of a family of periodic orbits the curve that gives the position of an orbit (say the coordinate \(x\) at the point of intersection of the orbit with the x-axis (\(y=0\))) as a function of the energy or the Jacobi constant ( Figure 11 ). From the stable families bifurcate higher order resonant families, e.g. 2/1, 3/1, 4/1, etc. It is known that in perturbed (i.e. non-axisymmetric) systems at the even resonances (2/1, 4/1, etc) the family \(x_1\) forms gaps. Three main types of gaps are shown in Figure 11, but more complicated types of gaps also exist. Similar phenomena appear also beyond corotation.
Near corotation we have four equilibrium points \(L_1\ ,\) \(L_2\ ,\) \(L_4\ ,\) \(L_5\) (\(L_3\) represents the center of the galaxy). For small perturbations \(L_1\ ,\) \(L_2\) are unstable and \(L_4\ ,\) \(L_5\) are stable.
Near \(L_4\ ,\) \(L_5\) we have two families of periodic orbits, the short (SPO) and the long period orbits (LPO) ( Figure 10 a ). These are connected by a complicated set of bifurcations. The ordered nonperiodic orbits near \(L_4\ ,\) \(L_5\) are either rings, close to SPO, or bananas, close to LPO ( Figure 10 b ).
Many effects of the orbits can be understood by considering in detail the dynamics close to resonances. The most important resonances are the inner and outer Lindblad resonances, where \(\frac{\omega_1}{\omega_2}=\frac{\kappa}{\Omega-\Omega_s}= \pm{2\over 1}\ .\) Close to these resonances the Hamiltonian can be written in action angle variables in the form \[\tag{30} H=\omega_1I_1+\omega_2I_2+aI^2_1+2bI_1I_2+cI^2_2+...+V_1 \]
where \(V_1\) represents a spiral or a bar \[\tag{31} V_1=\varepsilon_0\cos 2\vartheta_2~-\left(\frac{2I_1}{\omega_1}\right)^{1/2}~[\varepsilon_+ \cos(\vartheta_1-2\vartheta_2)+\varepsilon_- ~\cos(\vartheta_1+2\vartheta_2)]+... \]
The action \(I_1\) corresponds to the radial deviations, \(I_2\) corresponds to the angular momentum, \(\vartheta_1\) represents the angle along an epicyclic oscillation around the circular orbit, and \(\vartheta_2\) represents the azimuth around the center.
Away from all resonances we can use a canonical transformation to eliminate the terms containing angles and find \(H\) in the form \(H=\omega_1J_1+\omega_2J_2+...\ .\) Then \(J_1\) and \(J_2\) are integrals of motion. When \(\frac{\omega_1}{\omega_2}\) is close to 2 (inner Lindblad resonance) we can use a canonical change of variables and find an approximate Hamiltonian of the form \[\tag{32} \overline{H}=\omega_2J_2+\gamma J_1+...-\varepsilon_+\left(\frac{2J_1}{\omega_1}\right)^{1/2}\cos\psi_1=0 \]
where \(J_1,J_2\) are the new actions, \(\gamma=\omega_1-2\omega_2\) and \(\psi_1=\vartheta^\ast_1-2\vartheta^\ast_2\ ,\) with \(\vartheta^\ast_1, \vartheta^\ast_2\) the new angles. As \(\psi_2\) does not appear in this Hamiltonian the corresponding action \(J_2\) is a second integral of motion.
The integral \(\overline{H}-\omega_2J_2\) is a function of \(J_1\) and \(\psi_1\) which gives the forms of the orbits near the inner Lindblad resonance. These orbits are close to two perpendicular deformed ellipses. The nonperiodic orbits in this region form rings around these ellipses. Thus the resonant integrals explain the main forms of the orbits inside corotation in a galaxy.
Similar results appear near the outer Lindblad resonance and near corotation.
If the disc has non-negligible thickness, we may also consider 3D orbits which have a vertical oscillation with respect to the disc. Such orbits can explain particular forms of galaxies seen edge-on, like peanut galaxies and box galaxies.
Finally, a small proportion of stars continuously escape from a galaxy to infinity. In fact after a time interval equal to the time of relaxation \(t_{relax}\ ,\) the encounters between stars (collisional effects) produce an approximately Maxwellian distribution of velocities. In a Maxwellian distribution a proportion \(0.0074\) of stars have velocities greater than the escape velocity. Thus after a time t the proportion of escaping stars is equal to \[\tag{33} \frac{dN}{N}=0.0074 t/ t_{relax} \]
As a consequence a galaxy is completely dissolved after a time of the order of 1016 years.
In the case of a given potential \(V(\textbf{x})\) the stars inside corotation cannot escape from a galaxy, but the stars outside corotation may escape to infinity.
If the galaxy is not axisymmetric only the Jacobi constant ( (27) ) is kept constant, while the energy and the angular momentum undergo variations, that appear as random. Then after a sufficient time the energy may become positive and the stars escape.
If a system is close to integrable and there are invariant curves beyond corotation surrounding the galaxy the orbits inside them cannot ever escape. However in strongly perturbed systems there are no such invariant curves and the orbits may escape after a sufficiently long time. But even if the orbits are chaotic, nevertheless they are `sticky' and only after very long times such orbits can escape. Stickiness means that the chaotic orbits spend transient but long intervals of time in a restricted domain of the phase space around islands of stability or other invariant sets (like invariant manifolds; section 5.5), before filling a much larger chaotic domain, or escaping to infinity.
Construction of equilibria
The study of individual orbits in static potentials aims primarily to identify which types of orbits support the main observed morphological and kinematical features of galaxies. In order, however, to construct fully self-consistent stellar equilibria, one has to create an appropriate statistical mixture of many orbits ensuring that the resulting distribution function f is a solution of the collisionless Boltzmann equation ( (2) ). There are various methods to accomplish this goal. Some important methods are i) the hierarchy of Jeans' equations, ii) inversion formulae, and iii) the numerical construction of self-consistent equilibria (Schwarzschild). Finally (iv) one can make an `ad hoc' choice of a distribution function model.
Hierarchy of Jeans' equations
The collisionless Boltzmann equation ( (2) ) is, in general, very hard to solve in terms of all six arguments (positions and velocities). In practice, however, only the lowest moments of the distribution function can be compared to observations. In Cartesian coordinates, the moments are defined by: \[\tag{34} \mu_{i,j,k}(\mathbf{r})=\rho(\mathbf{r})\overline{v_x^i v_y^j v_z^k}= \int v_x^iv_y^jv_z^k f dv_x dv_y dv_z \]
The zero-th moment \(i=j=k=0\) is the density itself. The first order moments \(i=1,j=0,k=0\) (and cyclic permutations) depend on the mean streaming velocities in the three main directions, etc. Multiplying all terms of Boltzmann's equation with power terms \(v_x^iv_y^jv_z^k\) we obtain the so-called hierarchy of Jeans' equations, which relate the various moments defined by ( (34) ). The most important Jeans' equations are the continuity and the momentum equations which relate the zero, first, and second order moments of \(f\ .\) These equations must be supplemented with a closure condition, which is equivalent to the equation of state in hydrodynamics. Then, this system of equations can be solved and yields specific models that can be compared to observations. Special attention, however, has to be paid on aposteriori controls that the resulting models do not exhibit unphysical properties (for example, the density may turn to be negative in some region of the modelled galaxy, the velocity dispersion may show anomalous peaks etc.).
Inversion formulae
If we make some specific assumptions about the geometry and kinematics of an observed galaxy, we can deduce models of the distribution function by an inversion procedure starting from some functions determined by observations. The simplest such example, referring to spherical systems, is Eddington's inversion formula. Starting from the density \(\rho(r)\ ,\) we derive the potential through Poisson's equation \(\nabla^2 V=4\pi G\rho\) and then use the solution \(V(r)\) to find \(r(V)\ ,\) and finally \(\rho(V)\ .\) If we now assume that the distribution function is isotropic in the velocities, we find its form via \[\tag{35} f(E)={1\over\sqrt{8}\pi^2} {d\over dE}\int_{E}^0 {d\rho\over dV} {dV\over\sqrt{V-E}} \]
Various generalizations of Eddington's formula exist. For example,
the Ossipkov - Merritt formula for anisotropic spherical systems
reads:
\[\tag{36}
f(Q(E,L))={1\over\sqrt{8}\pi^2}{d\over dQ}\int_{Q}^0 {d\rho_Q\over dV}
{dV\over\sqrt{V-Q}}
\]
where \(Q=E+L^2/2r_e^2\ ,\) \(L\) is the angular momentum, \(\rho_Q=\rho(1+r^2/r_a^2)\) and \(r_a\) is a characteristic radius beyond which the velocity distribution becomes anisotropic.
In the case of axisymmetric systems, the inversion method yields the `even' part of the distribution function, i.e. the part depending on the energy and on the square of the angular momentum along the axis of symmetry (Lynden-Bell, Hunter, Dejonghe). The odd part, on the other hand, can be found by an inversion algorithm only if the azimuthal velocity distribution is known (Merritt).
Schwarzschild's method
A numerical method to examine the relative contribution of various types of orbits in models of stellar equilibria is Schwarzschild's method. This method tries to find whether there is an appropriate statistical combination of orbits in a model with fixed spatial density \(\rho(\mathbf{r})\) (called 'imposed density') such that the superposition of the orbits yields a 'response density' equal to the imposed density. To this end, a grid of initial conditions is specified in a properly chosen subset of the phase space (e.g. on equipotential surfaces). The orbits with these initial conditions are integrated for sufficiently long time intervals. This creates a ‘library of orbits’. By now dividing the space into a large number \(N_c\) of small cells, we calculate the time spent along each orbit inside every cell. We furthermore assign statistical weights \(w_o (o = 1, . . .,N_o)\ ,\) where \(N_o\) is the total number of orbits considered. Finally, we check whether there is a solution for these weights that satisfies the self-consistency condition \[\tag{37} \sum_{o=1}^{N_o}w_o t_{oc}=m_c~~~~c=1,\ldots,N_c \]
where \(t_{oc}\) is the time spent by the orbit labeled by o in the cell labeled by c, and \(m_c\) is the total mass in the cell c as derived from the imposed density function. We furthermore impose the constraint \(w_o\geq 0\ ,\) i.e. the weights can only be positive quantities. If a solution under the imposed constraints is found, we call it a self-consistent model of a galaxy. In practice, instead of a perfect solution we look for a minimizer of the difference between imposed and response density, which can be found by various algorithmic techniques such as linear or quadratic programming, non-negative least squares, the Lucy algorithm, or ‘entropy’ functional methods. Usually, the solutions found are non-unique, and further constraints can be introduced, for example, when kinematic data are available for a particular galaxy.
The method of Schwarzschild has been applied with great success in various models of axisymmetric or triaxial galaxies. An important question regards the relative contribution of regular and chaotic orbits to the so-constructed stellar equilibria. In early models of triaxial galaxies it was found that most orbits are regular. However, in models with central cusps or central black holes, it was found that the chaotic orbits contribute significantly to the self-consistent equilibria (by percentages reaching 50%). Finally, applications of variants of the same method in disc galaxies (Contopoulos) have revealed that chaotic orbits play an important role in the self-consistency condition near and beyond the corrotation region of barred galaxies.
Gas Dynamics
Up to about 10% of the total baryonic matter in disc galaxies can be in the form of gas or dust. The gas is present in the inner disc, but also at distances exceeding many optical lengths. The formation of galactic discs from cosmological initial conditions is itself attributed partly to dissipative processes taking place in the gas, during an early epoch of galaxy formation.
The main condition for a mass component of a galaxy to be described by gas dynamics is \[\tag{38} \lambda<<l \]
where \(l\) is the mean-free path, while \(\lambda\) is the scalelength over which the distribution function exhibits significant variations. In terms of the temporal evolution of the distribution function \(f(\mathbf{x},\mathbf{v},t)\ ,\) such a condition can be accounted for by including `collisional' (dissipative) terms in Boltzmann's equation, namely \[\tag{39} \frac{\partial f}{\partial t} + \frac{\partial f}{\partial \textbf{x}}\textbf{v}- \frac{\partial f}{\partial \textbf{v}} \frac{\partial V}{\partial \textbf{x}}= \frac{df}{dt_c} \]
where the term \(df/dt_c\) describes non-zero variations of the distribution function's convective derivative (the latter are zero only in pure stellar dynamics). The index c in the r.h.s. means `collisions', since collisional effects mainly account for the time change in the distribution function. For example, dynamical friction phenomena may become important when massive objects (e.g. black holes, giant molecular clouds or clusters) move in the background of stars and gas. The timescale in which such phenomena affect the global form of the distribution function is of the order of the collisional relaxation time (see section 6), but the local form can be affected in significantly shorter times.
By considering various velocity moments of the distribution function,
( (39) ) gives rise to specific forms of
hydrodynamical equations (i.e. the continuity and force equations).
The main role of the gas is that it produces secular evolution of the galaxies. Gas-dynamical phenomena of particular interest are:
i) Jeans instability: if we consider a volume of gas with homogeneous density \(\rho\ ,\) the time evolution of a density perturbation \(\delta\rho\) generated by any mechanism within this volume depends on the scale \(L\) in which this perturbation extends. Namely, if \(L>L_J=v_s(G\rho)^{-1/2}\ ,\) where \(v_s\) is the velocity of sound, the perturbation becomes unstable, i.e. its pressure can no longer sustain it against its self-gravity. The length \(L_J\) is called Jeans' length. A stellar analog of this phenomenon occurs with the velocity of sound replaced by the velocity dispersion of the stars. Jeans instability is important because it leads to the growth of even very small initial density perturbations, and various applications of it have been studied, ranging from the early epoch of galaxy formation to the later epochs, where e.g. the spiral structure is formed.
ii) Gas infall: giant gas clouds with masses comparable to those of dwarf galaxies, or giant gas streams, may cause accretion of mass in the centers of disc galaxies. The origin of the infalling gas lies mainly inside the disc, but a considerable fraction may come from the surrounding environment, in which case gas accretion is characterised as a remnant of cosmological infall. The accretion rate could be as high as 1 solar mass per year, resulting in significant secular evolution of the mass distribution in galaxies within one Hubble time. This, in turn, affects both the morphology and kinematics of a galaxy, but it also induces enhanced star formation and chemical (metalicity) evolution of galaxies.
In numerical simulations, the motions of infalling gas components often appear to be highly non-circular. However, the gas may settle to particular locations in the disc, close to the main disc resonances. In this way, rings can be formed near the inner and outer Lindblad resonances.
iii) Gas shocks: The spiral arms in disc galaxies are described as density waves embedded in the galactic disc. As a density wave propagates, large pressure gradients are formed in the wavefront, which may induce gas shocks. Such shocks are also accompanied by rapid star formation. In fact, observation of the phase lag between the loci of the maximum of the spiral density and of the density of suitable nearby tracers with relatively short lifetime (like open clusters of young stars), has been proposed as an observational method to determine the pattern speed of spiral arms.
iv) Hot coronae: elliptical galaxies are often surrounded by hot gas coronae (of temperature \(10^6\)K). The total mass in such coronae can rise up to \(10^{10}\) solar masses.
Numerical studies of gas dynamics offer insights regarding the role of dissipative processes in the secular evolution of galaxies. Besides general mesh hydrodynamical schemes, some numerical techniques well suited to the study of galactic gas dynamics are the flux-split (FS2, van Albada) and the sticky particle (Calberg) techniques. A very widely used technique is smooth particle hydrodynamics (SPH). This is a Lagrangian scheme where the equations of motion of individual `gas particles' (which correspond, physically, to large masses of gas) are integrated numerically. The SPH particles are considered to occupy a volume (e.g. a sphere) in space. Quantities like pressure, temperature and viscosity can be defined by taking weighted-averages inside the volume. The particles may lose energy due to cooling (caused e.g. by radiation of the gas). This is taken into account by considering the physical processes which cause the gas to radiate. The SPH method is very flexible for many purposes, compared to the less flexible but more accurate mesh hydrodynamical methods. SPH can be combined with N-body stellar dynamical methods (section 6) to explore the evolution of galaxies. It has been applied to studies of gas dynamics in isolated galaxies and in cosmological simulations of galaxy formation and the formation of large scale structures.
Finally, the comparison of photometric data in various wavelengths reveals to what extent gas dynamics follows or is differentiated from stellar dynamics in particular galaxies. This issue is important in the case of spiral structure, since we often find structures in the gas (e.g. spiral arms extending to large distances) that do not have a stellar counterpart, and vice versa.
Spiral Structure
The current paradigm of spiral structure in disc galaxies is based on the description of the spiral arms as density waves rather than material structures (i.e. always composed by the same stars). The origin of density wave theory lies in the observation that, if all matter is in nearly circular motion around the center of a disc galaxy, material non-axisymmetric structures embedded in the disc wind very quickly. This is due to the differential rotation of a galaxy, namely the fact that the angular rotation speed decreases from the center outwards. The average decrease with distance \(\Delta\Omega/\Delta R\ ,\) of order \(10 Gyr^{-1}Kpc^{-1}\ ,\) causes such structures to become tightly wound in a small time period (of order 1Gyr), acquiring pitch angles which are only a fraction of one degree. This is opposed to the fact that the pitch angles of observed spiral arms are substantially larger (a few degrees), leading to the conclusion that the observed spiral arms cannot be material.
The most natural solution to the above winding dilemma is found by adopting that the spiral arms of the spiral galaxies are density waves, i.e. the stars pass through the spiral arms but stay longer on the average close to them. Thus the spiral arms are not composed of the same stars always, but they are waves, i.e. the stars move with different angular velocities \(\Omega\) around the center, while the spiral arms (density maxima) rotate with a particular angular velocity \(\Omega_s\ .\) Stars inside corotation have \(\Omega>\Omega_s\) and overtake the spiral arms, while stars outside corotation have \(\Omega<\Omega_s\) and move in a retrograde direction with respect to the spiral arms.
The existence of spiral density waves has been established both by observations and by numerical N-body simulations.
Linear density wave theory
The density wave theory of spiral arms was developed first by Lindblad, and later by Lin, Shu, Kalnajs, Toomre, Lynden-Bell etc.
The linear theory of density waves starts with a given axisymmetric model that consists of the functions \(f_0,\varrho_0\) and \(V_0\) satisfying equations (2),(3),(4), and finds the first order perturbations \(f_1,\varrho_1\ ,\) and \(\!V_1\) that have a spiral form (or a bar form in special cases). In the simplest case instead of \(\varrho_0, \varrho_1\) we have the surface densities \(\sigma_0, \sigma_1\ .\)
The linear problem consists in finding self-consistent solutions of Eqs. (2),(3),(4) of the form \(f=f_0+f_1,\ \sigma=\sigma_0+\sigma_1\) (where \(\sigma\) is the surface density corresponding to a spatial density \(\rho=\sigma\delta(z)\ ,\) with \(\delta(z)\) the delta function) and \(V=V_0+V_1\) with \(\!V_1\) of the form \[\tag{40} V_1=Ae^{i(\varphi+\omega t-2\vartheta)} \]
where \(A=A(r)\) is the amplitude, and \(\varphi=\varphi(r)\) is the phase of the spiral, while \(m\) is the number of spiral arms (usually \(m\!=\!2\)). This problem leads to an integral equation that has eigenfrequencies \(\omega=2\Omega_s\) (for \(m\!=\!2\)), where \(\Omega_s\) is the pattern velocity of the spiral.
In the case of tight spiral arms the self-consistency condition leads to the Lin dispersion relation that relates the frequencies \(\frac{\omega_1}{\omega_2}=\frac{\kappa}{\Omega-\Omega_s}\) at a distance \(r\) with the inclination of the spiral arms (or local radial wavenumber \(k=\frac{d\varphi}{dr}\)), the local surface density \(\Sigma\) and the radial dispersion of velocities \(\sigma_R\ .\) The wavenumber is found as a function of the radius \(r\ .\) Thus we can derive the value of \(\varphi\ ,\) i.e. the form of the spiral arms all the way from the inner Lindblad resonance, past corotation, to the outer Lindblad resonance.
Assuming \(\phi\) to increase in the same sense as pattern rotation, the wave is called leading if \(k>0\ ,\) and trailing if \(k<0\ .\)
The linear density wave theory relies on an important local approximation, namely that if both the disc potential and surface density are analyzed in modes (like (40) ), the potential perturbation induced by one mode at a point of the disc depends only on the density perturbation for the same mode at the same point. This (WKB) approximation leads to a local relation between the two quantities which ignores how the density perturbation in a global scale affects the potential locally. At a distance \(R\ ,\) this approximation is valid to order \((kR)^{-1}\ ,\) thus, it holds better in the short wavelength limit (\(|kR|>>1\)).
Nonlinear theory - Termination of spiral arms
The linear theory is inadequate to deal with the regions near the inner and the outer Lindblad resonances, because there the perturbation \(\!f_1\) tends to infinity.
In fact if we expand \(f_1\) as in the case of the third integral, we find that \(f_1\) contains denominators of the form \((\omega_1 \mp 2\omega_2)\) (- at the ILR and + at the OLR) and this tends to zero when we approach the resonance. Therefore \(f_1\) is larger than \(f_0\) near the resonances and the basic assumption of the linear theory that \(f_1\) is a small perturbation of \(\!f_0\) is not satisfied.
In these cases we need a different perturbation analysis (Contopoulos), namely one has to use the resonant integrals of motion of the collisionless Boltzman equation (2). Away from all resonances \(f\) is a function of the generalized actions \(J_1,J_2\) as seen in section 2.3. Close to the inner Lindblad resonance \(f\) is a function of \(H\) and \(J_2\) (where \(H\) contains not only \(J_1\) and \(J_2\ ,\) but also the resonant angle \(\psi=\theta_1-2\theta_2\!\)). Thus a different expansion is valid in this case.
The use of a nonlinear theory of density waves allows us to describe the resonant phenomena near the Lindblad resonances and find the form of the spiral arms near these resonances.
The nonlinear effects are important also near corotation. In this case the Lin dispersion relation has no infinities but the amplitude of the wave tends to infinity as we approach corotation (Shu). Thus in this case also we have to use resonant integrals of motion in \(f\) and not the actions \(J_1,J_2\ .\)
In normal spirals the non axi-symmetric perturbation is weak, of the order \(2-10\%\ ,\) The linear density wave theory is applicable all the way to corotation if the amplitude of the spirals is of the order of \(2-5\%\ .\) But in the range \(5-10\%\) nonlinear effects due to higher order resonances are important. In particular the periodic orbits near the \(4/1\) resonance are similar to squares with 4 round corners or 4 loops. But while inside the \(4/1\) resonance the periodic orbits are oriented in a way supporting the spirals ( Figure 12 ), beyond the \(4/1\) resonance the orientation of the orbits is quite different and does not support the spiral arms. Thus the stronger normal spirals tend to terminate near the \(4/1\) resonance. This has been verified by many numerical experiments.
Preference of trailing waves
The observations indicate that most spiral galaxies host trailing spiral arms. On the other hand, the anti-spiral theorem (Lynden-Bell and Ostriker) states that for any collisionless steady-state solution of Boltzmann's equation representing trailing spiral arms, there is a symmetric solution (which is found be reversing the sings of all velocities in the distribution function) representing leading spiral arms. Thus, there are no steady-state spiral solutions, and the preference of trailing waves must be based on a time-dependence of the amplitude of the density waves such that trailing waves have a growing amplitude while leading waves have decreasing amplitude in time.
The variation of amplitude in time can be caused by both dissipationless or dissipational mechanisms. Some proposed stellar-dynamical mechanisms are:
Angular momentum transfer (Lynden-Bell and Kalnajs). The torques exerted between the spiral arms cause angular momentum to be transferred across the disc. In the case of trailing spiral arms, the transfer is from the center outwards, while in the case of leading spiral arms it is inwards, which leads to a number of non-physical effects.
Swing amplification (Julian and Toomre etc). If a small trailing wave-packet is formed, e.g. by shear instability, this moves towards the disc center with a group velocity \(v_g=d\omega/dk\ .\) After passing through the center it emerges as a leading wave, which propagates outwards and is transformed again to a trailing wave near corotation. A careful examination of the local epicyclic motions inside the wave shows that upon the transition to a trailing wave the perturbation is substantially amplified.
Resonance effects in the ILR (Contopoulos). When the density wave theory is implemented in the inner Lindblad resonance, it is found that, starting with a slightly growing imposed density of trailing waves both inside and outside the ILR, the response density forms a short bridge which smoothly connects these arms. If, however, we start by imposing leading spiral arms inside and outside ILR, the bridge formed by the response density is trailing.
Modal theory - theories of multiple pattern speeds
Two important open problems regarding spiral structure refer to i) the longevity and degree of quasi-stationarity of the spiral arms, and ii) whether all structures in a galaxy rotate with a common, or with multiple pattern speeds.
Regarding the assumption of quasi-stationary spiral arms, an important development going beyond the local assumptions of the Lin-Shu linear density wave theory arose by considering the problem of global modes in the galactic disc. This problem is hardly tractable in pure stellar dynamics, but numerical and semi-analytical solutions have been proposed where the effects of a cold gas component are also taken into account. This so-called modal theory (Lin and Bertin) describes evolving spiral or bar-spiral global modes. Their maintenance relies on the fact that the gas causes self-regulation, i.e. cooling of the disc, that compensates the stellar disc heating (the temporal increase of the random motions of the stars). The so-produced patterns are characterized as `quasi-stationary spiral structure' (QSSS), because the underlying global modes are long-lived, despite the fact that their superposition may lead to rather rapidly varying observed patterns.
On the other hand, despite the fact that the basic form of density wave theory examines the case where the spirals rotate with a unique pattern speed, there are indications, both from observations and from dynamical or N-body models, of the existence of multiple pattern speeds (Sellwood). In a well studied scenario, an inner bar exhibits fast rotation, while the spiral arms rotate more slowly. In fact, the spiral arms may be generated by recurrent instabilities due to the dynamical coupling of the bar with disc modes beyond corotation. The spiral arms in such models recurrently appear and disappear. In N-body simulations it is found that, when spiral arms are present, their pattern speed is close to a resonant relation with the bar's pattern speed.
Termination of bars - chaotic spiral arms
The non-axisymmetric perturbation in barred galaxies is much higher than in normal spiral galaxies (of order \(50\%\)). In fact, near corotation we can find ordered orbits near the Lagrangian points \(L_4,L_5\ ,\) while near the points \(L_1, L_2\) we have the interaction of many resonances and large chaos is generated in this region. The existence of large chaos near corotation tends to destroy the bar of a barred galaxy. In fact the orbits near corotation spread in an irregular way and do not support a continuation of the bar beyond corotation. The chaotic domain is larger if the bar is stronger.
Another reason for the termination of bars near corotation is that even if there are ordered orbits beyond corotation, these tend to produce density maxima perpendicularly to the bar, and not along the bar.
Numerical experiments and observations of real galaxies show that, in fact, the bars terminate a little before corotation.
However, in strong bars we observe spiral arms beyond the ends of the bar, starting close to the Lagrangian points \(L_1\),\(L_2\ .\) These spiral arms are rather different from the spiral arms of normal galaxies in the inner parts of the galaxies. In fact a systematic study has shown that the orbits of the outer spiral arms are chaotic. These spiral arms are still density waves, because they are not composed always of the same matter. The chaotic orbits produce maxima of density near their apocentra and pericentra, where \(\dot{r}=0\ .\) The apocentra and pericentra are close to the unstable asymptotic curves (invariant manifolds) of the unstable periodic orbits \(PL_1, PL_2\) around \(L_1, L_2\ .\) These manifolds have a trailing spiral form (U,U\('\) in Figure 13 a), and on the leading side they form an envelope of the bar (UU,UU\('\!\)). After longer integration times the unstable manifolds are seen to support a considerable part of the trailing spiral arms ( Figure 13 b ). In fact it has been found in many numerical experiments that the chaotic orbits tend to populate both the trailing spiral arms and an envelope of the bar while the main body of the bar is composed of ordered orbits (Voglis). A particular chaotic orbit supporting the bar and the spiral arms is shown in Figure 13 c.
These chaotic spiral arms may last for a considerable fraction of the Hubble time. Later on the orbits become more complicated and some orbits escape to infinity. But many chaotic orbits are `sticky', i.e. they remain close to the spiral arms and to the envelope of the bar ( Figure 13 c ) for very long times before escaping from the galaxy. Thus `stickiness', i.e. the tendency of chaotic orbits to remain confined close to invariant structures (islands of stability or invariant manifolds) appears to play a key role in understanding the structures formed in barred galaxies beyond corotation.
N-body Systems
N-body simulations
A particular method to explore the dynamics and evolution of a stellar system is by N-body simulations. These simulations are by definition self-consistent because the motions are due to forces generated by the masses themselves.
In the case of small stellar systems, like a cluster, we consider a system of N point masses, where N is of the order of 102 - 104 and we follow the orbits of all these points. In the case of large systems, like galaxies, we take N of the order of 105 - 108 but not the actual number of stars in a galaxy, which is of the order of \(10^{12}\ .\) In large systems we use various approximations to speed the calculations. E.g. we calculate in detail the forces due to nearby stars, but we consider the remote stars as forming a relatively small number of large masses affecting a particular star. The accurate N-body calculations for a given time T require \(N^2\) calculations of forces, but the approximate methods (like tree or mesh methods) require O(N log N) calculations. Implementation of the so-called `self-consistent field' method (Clutton-Brock) in isolated systems with simple geometry has, finally, led to the production of O(N) algorithms.
Similar calculations can be used for the gas dynamics of galaxies, that take into account also the effects of pressure. Such is the SPH (smooth particle hydrodynamics) method, in which the motions of fluid particles, acted upon by both gravitational and non-gravitational forces, are integrated.
These studies are mainly aimed at finding the evolution of stellar systems. E.g. we may study the collapse or merging of stellar systems until they form an almost stationary configuration. We may also consider the differences between the evolution of the stellar component and the gas component of a galaxy. Other problems of interest refer to the tidal effects between neighboring galaxies, or the collisions of galaxies and the formation of tails or supermassive central cores in some galaxies.
In such studies one considers two different types of relaxation, collisional and collisionless. Collisional relaxation is due to the close approaches of stars and their time scale is the time of relaxation. In a spherical system of stars of mass \(m\) with average radius \(R\) and density \(n_0\) it is \[\tag{41} t_{relax}=\sqrt{\frac{n_0}{Gm}}\frac{R^3}{\ln(N/2^{3/2})} \]
where N is the total number of stars. The relaxation time is much longer than the "dynamical time" \[\tag{42} t_0=\frac{R}{\sqrt{\langle v^2 \rangle }} \]
where \(\langle v^2 \rangle\) is the average of the squares of the velocities. On the other hand the collisionless relaxation takes place in only a few dynamical times. The rate of relaxation increases when the Lyapunov time \(t_L\) of chaotic orbits (section 2.1) is small.
In general, the evolution of an N-body simulation is characterized by the following phases which evolve in different timescales:
(a) Violent relaxation and/or growth of collective instabilities . These phenomena evolve in a timescale which is of the order of the dynamical time \(\!t_0\ .\) An abrupt collapse of a protogalaxy or a merging event (e.g. two colliding spiral galaxies forming an elliptical galaxy) is characterized by rapid fluctuations of the gravitational potential which result in a fast dynamical relaxation. The resulting stationary state gives a distribution of velocities with \[\tag{43} \langle v^2 \rangle=\frac{GM}{2R} \]
where \(M\) is the total mass of the system. On the other hand, the initial conditions (positions and velocities of particles) may be chosen close to an unstable equilibrium state, in which case the simulation checks the expected fast growth of instabilities and the form of endstates to which the N-body system is led.
(b) Virial equilibrium and/or secular evolution. Simulations of slowly or non-rotating 3D N-body systems may form states which are close to virial equilibrium, i.e. systems satisfying the virial equations (6) and (7). The main interest, in this case, is to examine the types of orbits and of the phase-space structures which support self-consistency all along the simulation. On the other hand, an N-body system may exhibit secular evolution. Two main cases are: i) addition of a central mass in simulations of elliptical galaxies (section 2), ii) bar-halo interaction in barred-spiral galaxies. The latter is an example of a live halo which exchanges angular momentum with the bar due e.g. to dynamical friction. The exchange results in a slowing down of the bar, gradually shifting the positions of resonances and affecting also the spiral structure beyond the bar.
(c) Mixing of chaotic orbits. This takes place over a Lyapunov time \(t_L\ ,\) and may be driven by small fluctuations of the N-body potential around its average form.
(d) Collisional relaxation. This type of relaxation becomes effective in a long timescale, which gives also an estimate of the time of dissolution of the system (section 2.3).
N-body simulations have become a very important tool for investigations in Galactic Dynamics (see specializing scholarpedia article on N-body simulations).
Violent relaxation
A common result found in many N-body simulations is that during the collapse or merging of stellar systems we have a fast approach to a quasi-stationary equilibrium state. This approach takes place in a timescale of the order of a few dynamical times. Thus, it cannot be attributed to two-body relaxation, whose timescale is much longer (of the order of the collisional relaxation time).
An interpretation of this phenomenon was developed in the 60s by the theory of violent relaxation (Lynden-Bell). This theory considers the statistical distribution of the so-called `phase-space elements', i.e. smooth elements of mass whose motion in phase space is governed by the collisionless Boltzmann equation. Then, it can be shown that the rapid time variations of the potential lead to a fast mixing of the phase-space elements. According to the `Lynden-Bell statistics' (which is similar to a Fermi-Dirac statistics but without mass segregation), violent relaxation should lead to an equilibrium coarse-grained distribution function which has the form \[\tag{44} \overline{f}=\frac{f_0~e^{-\beta[E-E_0]}}{1+e^{\beta(E-E_0)}} \]
where E is the energy and \(f _0\ ,\) \(\beta\) and \(E_0\) are constants. If \(\overline{f}\) is small this distribution becomes \[\tag{45} \overline{f}=f_0~e^{-\beta(E-E_0)} \]
and it looks like a Boltzmann distribution. However in the Boltzmann distribution, which is due to collisions (encounters), \(\beta\) is proportional to the masses of the stars \(m\ ,\) while the Lynden-Bell distribution is collisionless and \(\beta\) is independent of \(m\ .\)
The theory of violent relaxation has provided a paradigm of application of statistical mechanics in collisionless self-gravitating stellar systems. On the other hand, N-body simulations show that the quantitative predictions of the Lynden-Bell statistics represent well the endstates of only idealized systems with `smooth' cores (i.e. ones with a nearly flat density profile in the central parts). In fact, in real galaxies there is a variety of physical mechanisms leading to strong deviations from Lynden-Bell statistics. Examples are: (i) Partial mixing, i.e. the mixing of particles (or of the `phase space elements') is not complete. Obstructions to mixing are posed by the existence of additional constraints in the form of local approximately conserved quantities which restrict the allowable motions in phase space. In addition, the variations of the potential may decay well before there was sufficient time for mixing to become complete (this turns out to be true in particular in the outer parts of galaxies). (ii) Large local variations of the initial phase space density before relaxation. For example, when a massive galaxy merges with a smaller satellite galaxy, a substantial fragment of the smaller galaxy may be driven to the center of the more massive galaxy without any effective mixing taking place.
Phenomena like the above cause an effect of memory of the initial conditions (the conditions before relaxation) in the final state. The theory of partially relaxed systems represents an important open challenge for statistical mechanics in general. Various attempts towards such a theory have been proposed over the years, but their usefulness in comparison with N-body experiments is still unclear (see Efthymiopoulos et al. 2007 for a review).
Collective instabilities
Another topic studied via N-body simulations is collective instabilities in stellar systems. These are caused when a system satisfies some instability criterion, and lead to abrupt redistribution of the mass and/or velocities within a system. Important cases are:
(a) Axisymmetric instabilities: an equilibrium distribution of matter in a galactic disc is prone to axisymmetric instabilities when the velocity dispersion in the radial direction is small compared to the circular velocity. The threshold to instability is provided by Toomre's Q-parameter: \[\tag{46} Q={\kappa\sigma_R\over 3.36G\Sigma} \]
where \(\kappa\) is the epicyclic frequency, \(\sigma_R\) is the velocity dispersion in the radial direction, and \(\Sigma\) is the surface density. Instability occurs if \(Q<1\ .\)
(b) Non-axisymmetric instabilities in disc galaxies: depending on the form of the radial profile of the density distribution, disc galaxies may develop a non-axisymmetric instability leading to the formation of bars. This effect was observed in N-body simulations (Hohl). The threshold of instability is given by the so-called Ostriker-Peebles criterion: \[\tag{47} {T\over W}> 0.14 \]
where T is the rotational kinetic energy and W the potential energy of the disc. Transient spiral arms and/or ejection of material from the galaxy may occur after the onset of this type of instability.
(c) Radial instability in spherical systems: these may lead to a collapse or a number of radial pulsations of a stellar system. The main criteria for the stability of spherical systems have been given by Antonov. Briefly, a condition for stability with respect to both radial and non-radial perturbations is that the distribution function \(f\) should be a decreasing function of the energy E. However, a system develops a radial instability if \[\tag{48} {d^3\rho\over dV^3}>0 \]
over a radial domain, where \(\rho\) and \(V\) are the density and potential of the spherical system respectively.
(d) Radial orbit instability. This is manifested in spherical galaxies, when the velocity dispersion in the radial direction of motion is much smaller than in the transverse direction. The radial orbit instability leads to the loss of spherical symmetry and the formation of triaxial galaxies. Various thresholds for its onset have been proposed in the literature (Polyachenko and Shukhman, Merritt, Palmer and Papaloizou). One such criterion (Polyachenko and Shukhman) is \[\tag{49} {2T_r\over T_t}>1.75 \]
where \(T_r\) and \(T_t\) are the kinetic energy in the radial and transverse directions respectively. The radial orbit instability is one of the earliest types of instability observed in N-body simulations (Henon), and it is believed to affect mainly the dark haloes of galaxies by turning their shape from spherical to triaxial.
(e) Warp or bending instabilities. These are analyzed in terms of the growth of bending modes, i.e. collective vertical oscillations of a thin disc.
N-body simulations have also been used in the study of live haloes (bar-halo interaction) and of the formation and secular evolution of galaxies.
Finally the N-body simulations can be viewed as providing an experimental method to check the theoretical considerations in galactic dynamics.
Further reading
Bertin, G. and Lin, C.C.: 1996, Spiral structure in galaxies: a density wave theory, MIT Press, Cambridge, Massachusetts.
Bertin, G.: 2000, Dynamics of Galaxies, Cambridge University Press, Cambridge.
Binney, J. and Tremaine, S.: 2008, Galactic Dynamics, second edition, Princeton University Press, New Jersey.
Boccaletti, D. and Pucacco, G.: 1996, Theory of Orbits, Springer, Berlin.
Combes, F., Boissé, P., Mazure, A., Blanchard, A.: 2001, Galaxies and Cosmology, Springer, Berlin.
Contopoulos, G., 2004, Order and Chaos in Dynamical Astronomy, Springer, New York.
Contopoulos, G. and Patsis, P.A. (Eds): 2007, Chaos in Astronomy, Astrophysics and Space Science Proceedings, Springer, Berlin.
Efthymiopoulos, C., Voglis N., and Kalapotharakos, C., 2007, Special Features of Galactic Dynamics, Lecture Notes in Physics 729, 297.
Fridman A.M. and Polyachenko V.L.: 1984, Physics of Gravitating Systems, Springer, Berlin.
Palmer, P.: 1995, Instabilities in Collisionless Stellar Systems, Cambridge University Press, Cambridge.
