“Animals do not care about the evolution of the universe. Nor do many humans.”

— Stephen Hawking —

Galaxy Models

In this chapter, we will show how to create a computer model of a galaxy in order to study galaxy dynamics. The creation of initial conditions for galaxies is some kind of black magic. A century ago, nobody even knew that galaxies exist. The Milky Way galaxy and the universe were synonymous. Now, astronomers believe that galaxies mainly contain the matter of form other than we know from our everyday life.

5.1      Galaxies

A galaxy is a massive ensemble of stars and other material orbiting about a common center and its constituents held together by the mutual gravitational interaction. Galaxies come in the variety of global shapes and internal morphologies. They can be, however, broadly classified into two major types: elliptical and disk galaxies.

To our understanding, before the formation of first galaxies were in the universe immense shapeless low-density diffuse clouds of light elements. Due to gravity, they started to contract and collapsed into smaller fragments of the size of galaxies. During runaway gravitational collapse, where stars formed quickly in large numbers and started to shine at once before the gravitational collapse finished, an elliptical protogalaxy was born. The energy of gravitational collapse was not lost (dissipationless collapse) and was converted to the chaotic motion of stars. If the primordial gas was shrunken by the gravity more slowly, the gas have had enough time to start rotation and settled into a regular disk galaxy where it formed stars. Galaxy’s appearance therefore mirrors its formation. This is known as the gravitational collapse theory of a galaxy formation. In Chapter 6 will be shown that a subsequent evolution is important as well.

5.1.1     Elliptical galaxies

Elliptical galaxies are not flat ellipses, but they are three-dimensionally elliptical (water-melon like). Orbits of stars do not show a systematic rotation in such galaxies and are largely chaotic. Elliptical galaxies are smooth, featureless, almost spherical. Elliptical galaxies are slowly rotating objects (Bertola and Capaccioli, 1975) supported by pressure, i.e. their dynamics is dominated by the irregular motion of stars (elliptical galaxies are kinematically hot).

Some of the elliptical galaxies are thought to be very old, because they contain red giant stars that are known to be old and very little or no gas and dust. The lack of blue young stars shows that the elliptical galaxies do not form stars currently. These elliptical galaxies reached their shape after the collapse of initial protogalaxy cloud. Some ellipticals might be newer and result from a more recent evolution (see Chapter 6).

Elliptical galaxies are both the most massive galaxies containing up to a few trillion of stars (central dominant or cD galaxies at the centers of galaxy clusters) and the least massive galaxies with a few million of stars.

Hubble space telescope (HST) imaging of massive elliptical galaxies revealed super-massive black holes in their centers. Giant elliptical galaxy M87, for example, contains central supermassive black hole of 2 to 3 billion solar masses ().

Black holes

Black hole is an extremely dense object, where the gravitational field is so powerful that nothing can escape. We can categorize black holes through their mass as stellar black holes (BHs) with masses , intermediate black holes (IMBHs) with masses , and supermassive black holes (SMBHs) with masses  (Miller, 2006).

5.1.2     Disk galaxies


A thin disk consists of relatively young stars, open clusters (loose clusters of stars), middle age stars like the Sun, gas and dust (interstellar matter or ISM). Observations of disk galaxies are showing that these galaxies have the very thin disk whose radius is of order 10 kiloparsecs and thickness is of order 100 parsecs.

The most of the stars in a disk galaxy travel on nearly circular orbits. Disk galaxies are rotationally supported, i.e. the regular circular motion of stars (disk galaxies are kinematically cool) dominates them. Stars in disk galaxies usually rotate with a constant velocity in the range of  to  with only a low dispersion of the velocity of order . Therefore, the angular momentum must govern the structure of disk galaxies. The rotation of the disk prevents its gravitational collapse in a radial direction. We all live in the Milky Way galaxy, where stars usually rotate with a velocity of  and with a low velocity dispersion of about .

It is estimated that there are approximately  stars in the disk of the Milky Way galaxy in total, most of these stars have a slightly lesser mass than the Sun have. The Sun lies approximately 8.5 kiloparsecs from the center of the Milky Way. Disk galaxies contain mainly blue stars that are very massive and are known to live shortly. Because these young stars still exist in disk galaxies, these galaxies are referred to as young galaxies. These stars are born mainly in spiral arms. The spiral arms look like that they contain a more mass than their surroundings. However, there is at most 5 percent more of the mass in the spiral arms than outside of them. The formation of young massive bright blue and ultraviolet stars makes the spiral arms look so bright. The gas is compressed in the location of spiral arms and therefore the star formation takes a place mainly in arms. Mature disk galaxies now contain of order  solar masses of gas ().

Figure 5–1: A face-on view of a disk galaxy with spiral arms, bar and bulge.


A central bulge is a small spherical component dominated by old red stars. This galactic component of a disk galaxy resembles a miniature elliptical galaxy. The bulge surrounds the central super massive black hole. The Milky Way galaxy harbors the black hole with the mass of about 3.6 million solar masses (Eisenhauer et al., 2005). Supermassive black holes lurk at the centers of most, if not all, galaxies (Silk and Rees, 1998).


Bars are a common feature among disk galaxies. These bars have various morphologies in the sense of length, thickness, and shape. From surveys of spiral galaxies is estimated that about 75 percents of spiral galaxies have bars or ovals (Grosbøl et al., 2002). The Milky Way galaxy also has a bar (Binney, 1995).

Stellar halo

A roughly spherical halo of galaxy contains globular clusters (GCs), i.e. isolated dense stellar clusters of millions of stars, and other ISM. Like in elliptical galaxies, the motion of stars in globular clusters is chaotic. Globular clusters usually contain old red (giant) stars and are therefore very old. The stellar halo has a mass in the range of 15 to 30 percent of the mass of disk. The diameter of the halo is approximately the same as the diameter of the disk.

Figure 5–2: The disk galaxy as on Figure 5–1, but seen edge-on.

5.2        Modeling galaxies

Galaxies did not appeared suddenly out of nothing as we will show in Chapter 6. However, it is necessary for the sake of clarity to make some simplifying hypotheses. In this chapter, we will study the evolution of galaxies by performing simulations of isolated systems that are in the steady state (equilibrium[42]).

Elliptical galaxies contain mainly stars and very little gas. They can be readily approximated with mass bodies representing stars. Disk galaxies contain larger quantities of gas, but still there is a significantly larger proportion of stars than of the gas. Galaxy is in the first approximation a large collection of stars or with even greater idealization, the ensemble of mass bodies hold together by the gravity. We can suppose that many of the properties we study in galaxy dynamics can be understood using a very simple approximation composed of following models.

Mathematical review

First, I will review some useful mathematical notions used throughout the text.

Vector quantities like  are usually expressed as , where  are positions (components of a vector) along unit base vectors (directions) .

The summation convention is used during repetitive additions. For example, the equation


can be rewritten as



We will encounter Hamilton’s operator  (pronounced “nabla” or “del”). It is defined as



We will also encounter Laplace’s operator . It is defined as



5.3        Gravity

In secondary school physics we are taught that the gravity is a force acting remotely and instantaneously between all bodies as was stated by Sir Isaac Newton[43], based on his observations of nature and reasoning. According to Newton’s theory “every particle of matter in the universe attracts every other particle with a force which is directly proportional to the product of their masses and inversely proportional to the square of their distances apart”.

In the mathematical description of gravity, we will start with a scalar gravitational potential field rather than a vector force field. Suppose we have some arrangement of mass points fixed in space and we are interested in their cumulative effect on some given point in space. Total gravitational potential of N bodies with masses  in a certain point in space is additive and is given by the total sum of all potentials from individual bodies  at distances  from the given point in space as


The configuration of N point masses (the N-body system) can be generalized to a continuous mass density distribution function , which is a function of position . If we write  as a function of coordinates  then  is the mass contained in a little box of volume , located at the point . A mass density  in a specific place is given by the sum of mass sources (point masses) in an infinitely small volume element. The mass density distribution function  is the function describing the distribution of matter (stars) in the system (galaxy). A transition from a discrete to continuous form can be approximately expressed as .

In a similar way, the gravitational potential generated by the continuous mass distribution instead of individual bodies as sources of the potential can be expressed continuously. Summation (5.5) is turning into integration in places (macroscopically) continuously filled with matter.

The gradient of the potential is then the gravitational force acting on a body with a unit mass, apart from a minus sign. The exploring body of the unit mass will experience at any point given by a position  force



This quantity  is called the intensity of gravitational field.

With continuous functions, we may use Poisson’s equation (Greiner, 2004)


Here we get a connection between the gravitational potential  and the mass density . Poisson’s equation serves as one of the principal governing equations in our examination of stability, structure and dynamics of self-gravitating systems. The mass density required for the generation of the potential can be found by solving Poisson’s equation. On the other side, from the known mass density, we can determine the gravitational potential through Poisson’s equation.

Different components of galaxies (disk, bulge… and halo) have their own mass density function within the global gravitational potential. Poisson’s equation gives us a relation between them. For a self-gravity model, we must find the global gravitational potential, in which the mass density of all parts of a galaxy is combined as



5.4        Initial conditions for spherical systems

As was described in Chapter 4.1, initial conditions for the N-body system must be set up. The most convenient is to start with elliptical galaxies. Systems with a spherical symmetry can be described by functions depending only on a one-dimensional distance  from the center of the system. In the centrally symmetric mass distribution with the spherical symmetry are the mass density  and gravitational potential  same in each direction with the origin at the center of the system. The spherical symmetry therefore implies  and . Poisson’s equation in one dimension is expressed as



5.4.1     Uniform sphere

A uniform sphere is one of the simplest approximations of the spherical system with the uniform density distribution function (volume density) given by



where  is the mass of the sphere and  is its radius. As we would expect from the mass density distribution function , outside the sphere potential falls off as in Keplerian case. This is known as Newton’s theorem. It states that “a body inside a spherical shell of matter experiences no net gravitational force from the shell” (the mass exterior to the shell has no effect) and that “the gravitational force on a body lying outside a closed spherical shell experiences a gravitational force which is the same as if all the matter in the shell was concentrated at a point at its center” (Greiner, 2004).

5.4.2     Isothermal sphere

Isothermal sphere has interesting property since its rotational velocity is unchanging with radius (it is constant). Therefore the isothermal sphere is useful for modeling flat rotational curves (more on Chapter 5.11). The mass density distribution function of isothermal sphere is given by



where  is a central density and a is a scale length. This distribution function (see
Figure 5–3) gives the system with an infinite mass and must be therefore truncated at some distance from the center. Cut-off radius must be imposed for a finite total mass.

5.4.3     Plummer’s and Hernquist’s model

We should now turn to more realistic density distribution functions that mimic spherical components of galaxy. Plummer’s (1911) model is such a simple approximation with the mass density distribution function given by



Another more realistic model is Hernquist’s (1990) model given by



Plummer’s and Hernquist’s models are members of the same family of models (Evans and An, 2005). They can be both described by the mass density distribution function as



where parameter p is a positive number. The case  is Hernquist’s model and the case  Plummer’s model.

5.4.4     Jaffe’s model

Another realistic model is Jaffe’s (1983) model expressed as



5.5         Units and scales

It is often useful to reduce units and equations describing a physical system to a dimensionless form, both for physical insight and numerical convenience. Let us imagine that we are using the gravitational constant in SI units that is of order . However, typical distances in the world of galaxies are of order  in SI units. A finite numerical representation of real numbers in a digital computer will most likely yield to large numerical round-off errors when we will compute with numbers of tens orders different. Therefore, we will use non-dimensional units throughout our simulations. We adopt  for Newton’s gravity constant, a total mass of system  and a scale length  is also set equal to 1. As a general rule, all quantities in this text are given in these intrinsically scale free units unless otherwise noted.

It is useful to summarize some commonly used units:


Let us consider a practical choice of units on measures of galaxies. A typical disk galaxy has a size of orders of kiloparsecs, contains hundreds of billions of stars, and a rotational velocity in such galaxy is usually up to . The models may be therefore compared with real galaxies using, for example, the following scaling


where we employed a convention to write  for the unit of quantity .  is a unit length,  is a unit mass, and  is a unit velocity. With these, the unit time can be derived as



Figure 5–3: Density profiles of spherical components for modeling galaxies. All have .

Figure 5–4: Density profiles of spherical components for modeling galaxies on logarithmical scale axes. All have .

5.6        Positions of bodies in galaxy

A probability  that we will find a body (i.e. star) in the volume of spherical shell limited by radii  relates to the radial mass density distribution function  as



First, we will introduce the amount of mass in the spherical system bounded by the sphere with a radius . This quantity is called a cumulative mass distribution  and is given by



In numerical formulation is the density profile  numerically integrated with the radius  changing between a minimal radius equal to zero and a maximum radius . The cumulative mass distribution  of a model is tabulated on a grid with points spaced logarithmically in . To check and compare numerical results with analytic solutions, we have plotted the analytic and numeric  of Plummer’s model to see the difference (Figure 5–5).

Figure 5–5: The cumulative mass of Plummer’s model computed in an analytic and numeric form. The numeric representation is so precise that no deviation from the analytic expression can be seen. The numeric representation (green curve) completely covers analytic representation (red curve).

The cumulative mass distribution  is essentially a probability distribution function (PDF). Numerical libraries contain a variety of random number generators, but how can we realize an arbitrary PDF?

Given the known PDF, it can be realized by a random (Monte Carlo[44]) sampling from the PDF. We can find the maximum  of the function  (overall mass of the system of bodies); we also know the maximum radius . Then, we generate a couple of random numbers with a uniform distribution. We generate a random number R with the uniform distribution in an interval  and expand it (multiply by ) to  range. Then we generate a second random number P with the uniform distribution in interval  and expand it (multiply by ) to  range. If the value of function  is equal or smaller than the random number P (so that the number P is under the graph of point ) then we will accept the distance R. Otherwise, we reject this number, generate again two new random numbers R and P and perform the test again. The method is called acceptance-rejection method. It was used by John von Neumann[45] and is therefore sometimes called “von Neumann’s method”.

We may discover also other method for sampling from the arbitrary PDF. We may invert the function  to a function  and thus obtain . We can find the maximum value  of the function  for . Then we generate the random number M from the uniform distribution and expand it to interval . Finally, we will receive the position of the body by evaluating .

Up to now, all functions were one-dimensional. Now, we must scatter bodies to all three dimensions. The initial radius  of a body is randomly determined according to the radial surface density distribution function . The initial angles  and  are drawn uniformly between  and . Particles are then positioned in three dimensions using


5.7        Velocities of bodies

To have a body on a stable orbit, its kinetic energy  must be exactly the same as a potential energy , i.e.,  or [46]. The kinetic energy is defined as , where m, v are the mass and velocity of the body, respectively. The potential energy is defined as , where  is the distance of the body from the center of the system. The potential energy characterizes a mutual force between bodies – the potential energy is shared at least between two bodies. The gravitational potential  is always negative. The potential energy  is resulting from the interaction of the body and the companion system.

The velocity for the stable orbit can be derived as


The body will stay in the system until its kinetic energy  will not be greater than the potential energy  of gravity that binds this body to the system. When , the body will leave the system with so called escape velocity . This velocity can be easily derived for the spherically symmetric system, when we take into account Newton’s theorems. Then the whole system is looking like the mass point of a mass M with the gravitational potential  and the escape velocity required for the body to escape infinitely far from the gravitational field will be


The body will stay in the system for .

Phase-space distribution function

An equilibrium spherical system depends on the phase-space coordinates  (Lynden-Bell, 1962). The structure of a galaxy is defined by a distribution function  in the imaginary space called phase-space with 6 dimensions (three spatial coordinates for positions and three corresponding velocity components) at each moment (time) t. The phase-space is given by the product of configuration space (positional space, represents generalized coordinates) and velocity space (represents generalized momentum).

Such a system of bodies is subjected to the potential  that is completely described by its phase-space distribution function . It is a probability density of finding some mean number  of stars  in a small box  near a position , i.e. between  and ,  and ,  and , and velocities in the range  about , i.e. between  and ,  and ,  and . The probability density is never negative.

We have already found the position of the star in the system (a 3D position in the configuration sub-space of 6D phase-space) as was described in Chapter 5.6. Now, we must find the range of velocities  in the phase-space for the concrete position and also a correct statistical weight that tells us what probability distribution function is for . As was shown in Equation (5.21), the escape velocity is given by the gravitational potential energy that comes from the gravitational potential of the system at the given point of space. The velocity distribution function at the position of each particle can be found from the distribution function . The distribution function f is assumed to be zero  for energies  greater than zero (while assuming a unit mass), because the stars with the positive energy will leave the system.

We must find the concrete phase-space distribution function  for the desired models. The probability isotropic distribution function f of any spherical model depends only on the overall specific energy E (energy per unit mass) of the system of stars and is available as (Eddington, 1916)[47]



Abel integral transform is required to relate density profile  of models and the distribution function  to the potential generated by . However, analytical distribution functions are known for only a handful of models, such as e.g. Plummer’s model. In order to generate any model that is more realistic, one has to find the steady state distribution function  that reproduces the desired density  numerically.

Gravitational potential

In order to solve Equation (5.22), we must find the gravitational potential  connected to the mass density  (that we have for various spherical models) through Poisson’s equation (5.7). Let us remind Poisson’s equation in one dimension given by



The cumulative mass distribution  relates to the gravitational potential  through Poisson’s equation as follows


In the numerical formulation is the gravitational potential  determined by numerical integration of the cumulative mass distribution  and is tabulated on a grid between the minimal radius equal to zero and radius equal to ten times the maximum radius  with grid points spaced logarithmically in . The cumulative mass distribution  is obtained by the interpolation of the grid points.

Determination of distribution function

The integrand in Equation (5.22) depends on the derivative of the density . We must therefore express mass density  as the function of the gravitational potential . Both the mass density  and the gravitational potential  are obtained from a numerical interpolation of  and  at the logarithmically spaced points  and acquired values are placed on a grid that represents . The numerical derivation is then performed on the function . Numerical integration is applied to the integral from Equation (5.22) and the result is placed on the grid of  evaluated at logarithmically spaced points . Finally, Equation (5.22) must be differentiated; and we get the distribution function .




Figure 5–6: The gravitational potential of Plummer’s model computed in an analytic and numeric form. The numeric representation is so precise that no deviation from the analytic form can be seen.



Figure 5–7: Numerically computed gravitational potentials of various spherical models.


Figure 5–8: The mass density distribution  as the function of the gravitational potential .



Figure 5–9: Numerically computed distribution function of various spherical models.


Figure 5–10: The distribution function  of the Plummer’s model computed in the analytic and numeric form. A minute deviation can be seen as the function is approaching zero.

Once the N-body system is generated, its realization must be verified. Every closed system held together by the gravity must globally obey the virial theorem[48]. The virial theorem states that the mean values of the kinetic energy  and potential energy  are related as  or  in the closed system. The total energy  is conserved.

Given a total mass M and the mean kinetic energy  of the system of stars, we may get from the definition of the kinetic energy a characteristic velocity



We can also define a characteristic length (Meylan and Heggie, 1996)



These characteristic quantities are sometimes known as the virial velocity and radius. Their ratio is an estimate of time that a typical star takes to cross the system. The crossing time is a measure of the time that take for a star to traverse the diameter of the system. The crossing time is defined (Meylan and Heggie, 1996) as



where  is a measure of the size of the system and  a measure of the mean stellar velocity. We choose  and  to be equal to  and , respectively.

5.8        Spherical galaxies

Our first objective is a simulation of elliptical galaxies with an eccentricity of zero (spherical galaxies). We will perform the simulation of two spherical systems – Plummer’s and Hernquist’s models.

We will use a program for initial conditions generation to create a spherical galaxy consisting of stars with a position, velocity and mass each. We will evolve initial conditions with the Barnes-Hut N-body simulation code with the opening angle .

Table 5–1: Parameters of simulation for Plummer’s and Hernquist’s model.

Model realization is in the dimensionless system of units. Gravitational forces were calculated after smoothing the mass distribution using the softening length . Model is comparable to a real-world dwarf spherical galaxy with following radius, mass and Newton’s gravitation constant


The time unit is comparable to , the timestep to  and the overall simulation with  steps covers .



224 million years

448 million years




671 million years

895 million years

1.12 billion years

Figure 5–11: A time sequence of Hernquist’s model in the xy plane. Practically no evolution can be seen.


224 million years

448 million years




671 million years

895 million years

1.12 billion years

Figure 5–12: A time sequence of Hernquist’s model in the yz plane.


224 million years

448 million years




671 million years

895 million years

1.12 billion years

Figure 5–13: A time sequence of Plummer’s model in the xy plane. Practically no evolution can be seen.


224 million years

448 million years




671 million years

895 million years

1.12 billion years

Figure 5–14: A time sequence of Plummer’s model in the yz plane.

Stability test of the initial models shows a little evolution in these models. Both Hernquist’s and Plummer’s models are in the dynamic equilibrium from the beginning.

5.9        Initial conditions for flattened systems

A spiral galaxy is flat and can be represented by a flattened potential. The disk of the spiral galaxy is generally very flat and we will describe them as infinitely thin. We will suppose that disks are axisymmetrical and therefore, initial conditions will be set up in cylindrical coordinates.

An initial angle  is drawn uniformly between  and . In the xy plane are disk particles positioned using


Stars in the disk have only circular velocities (so called cold disk, where no random motions are induced). Disk rotation is generated to be counter-clockwise in the xy plane (right handed Cartesian coordinate system is completed with z-axis and disk rotation is chosen to by given by the right hand rule with a thumb along z-axis).

5.9.1     Kuzmin’s disk

A reasonable approximation for the disk galaxy is Kuzmin’s disk. Kuzmin found in 1956 the potential of an infinitesimally razor-thin disk with Plummer’s potential in the plane of the disk. The initial mass density distribution (surface density)  of Kuzmin’s disk in the disk plane is



where  is a distance from the center of the disk and  is a radial scale length. This mass density belongs to a more general family of disks introduced by Toomre (1963).

5.9.2     Exponential disk

Bodies in an exponential disk are distributed in the plane with the surface mass density distribution given by



5.9.3     Kepler’s disk

Kepler’s disk has the potential of a point mass  given by . Its rotation is that obeyed by the planets in the solar system and Kepler’s three laws. A circular velocity computed after introducing Kepler’s potential is



Figure 5–15: An initial setup of Kepler’s disk with circular velocities.

5.10   Disk galaxies

Kepler’s model

As the first approximation, we can model a disk galaxy with Kepler’s model. As the planets in the solar system revolve around the Sun, also the stars (with planetary systems) in the disk galaxy are revolving on (nearly) circular orbits around the supermassive black hole at the center. The two massive bodies (SMBH and single star) are essentially unaffected by the third body (other stars).

Let us assume that the galaxy is made of two components: the stellar disk with uniformly distributed stars and the super-massive black hole. A model realization is in the dimensionless system of units. We will evolve initial conditions with the Barnes-Hut N-body simulation code.

Table 5–2: Parameters of the Kepler’s disk simulation.

The conversion to physical units is performed by scaling the Kepler’s digital galaxy to the system with known physical properties. If we scale it to the disk of the Milky Way galaxy, the unit length, unit mass, unit velocity, gravitational constant and unit time are


The timestep is therefore equal to  and the overall simulation covers .


205 million years

409 million years







614 million years

818 million years

1.02 billion years

Figure 5–16: A time sequence for Kepler’s disk model in the xy plane. The disk is ideally flat in the yz plane.

Figure 5–17: A rotation curve[49] of Kepler’s model.

5.11   Rotation curves

We will now compare a rotation curve from the simulations of Kepler’s disk with a rotation curve of a real galaxy as measured by a real experiment. Such work can be carried out with a radio telescope that traces neutral hydrogen (HI) atoms and carbon monoxide (CO) molecules. HI atoms are hard to observe, but CO molecules are copying the presence of HI atoms and their motion. Spectroscopy and the Doppler Effect permits us to measure the speed of gas orbiting around the center of the galaxy.

These observations of galactic rotation curves in disk galaxies have revealed that rather than falling-off as we can see in the plot of our computer simulation of Kepler’s disk (see Figure 5–17) and that is in an agreement with Newtonian physics and observed mass, the rotational velocities remain constant at large radii. This is valid for rotation curves of the Milky Way galaxy and other real galaxies (Fuchs et al., 1998).

Figure 5–18: Rotation curves: Theoretical (based on Newton’s law of gravity and visible matter) and experimental (based on observations).

Based on an astronomical observation of a mass distribution in disk galaxies can be inferred that the number of stars is dropping towards the edge of the stellar disk. From Newton’s law of gravity, we obtain


where  is the orbital velocity of a star at the distance  from the center of the galaxy. The cumulative mass inside radius  is



If the orbital velocity  is constant with increasing  (as is known from observations) then for
G = const. the only term in Equation (5.33) that can change with increasing  is . The cumulative mass  must grow even in vast distances from galaxy center. However, observed light or luminosity  that is believed to trace the matter distribution tends to a finite limit as we reach the edge of the luminous material in the galaxy. Therefore, the flatness of the rotation curve in large distances from the center of the disk, where number of stars is beginning to drop, cannot be explained by the visible mass alone.

Dark matter

It is supposed that invisible matter supplies this missing gravity. Jan H. Oort (1932) studying the motion of stars in the Milky Way  galaxy was the first one, who noticed that some matter in our galaxy is “missing”. The flat rotation curves of spiral galaxies at large galactocentric distance suggest that the halo of the invisible matter might exist around spiral galaxies; this dark halo might continue to contribute to  out to very large radii. Fritz Zwicky[50] studied the motion of galaxies within Coma cluster of galaxies in 1933. He discovered that galaxy velocities are so large that they should fly out of the galaxy cluster long ago. Since the cluster shows no signs of flying apart, some force of invisible source must hold it together. It is generally supposed that some additional matter invisible to us must gravitationally bound these galaxies together (Zwicky, 1937).

However, idea of dark matter (DM) was not accepted till extensive observational results of Rubin and Ford (1970). Visible parts of the spiral galaxy and the gas layers cannot account for the observed high rotation velocities in the outer parts of spirals, if the stellar mass-to-light ratio is constant with the radius. This is an indirect evidence that luminous portions of galaxies represent only a small fraction of the total galactic mass. Pioneering simulations of galaxies performed by Ostriker and Peebles (1973) also revealed that the inclusion of the massive halo will stabilize a galaxy disk against strong axisymmetric instability. Otherwise, the galaxy simply flies apart. The Milky Way and other disk galaxies may be therefore embedded in halos of unseen matter.

5.12   Activation of the disk

Orbits of stars in real disk galaxies are not perfectly circular. We will create more realistic models of the disk galaxy with a hot disk, where a velocity dispersion (fluctuation, deviation) is introduced. A way to quantify velocity dispersion in the disk is the Toomre’s (1964) parameter  expressed as



where  is the radial velocity dispersion,  is the frequency of small oscillation about the circular radius (local Lindblad epicyclic frequency) and  is the stellar surface density. The stability criterion against axisymmetric perturbations requires for the stellar disk to set . For , the disk is locally unstable and  corresponds to the cold disk. The epicyclic frequency is opt as , where  is an angular velocity .

Tangential (azimuthal) velocity dispersion  is expressed as



Vertical velocity dispersion (dispersion in a direction, which is perpendicular to the disk plane) is chosen as the part of tangential velocity dispersion as



From Equation (5.36) can be seen that even if the initial mass density of the disk is planar, vertical dispersion will lead to the thickening of the ideally flat disk.

Hot disk realizations

We will create a disk with kinematically hot (supported by random motions, rather than the rotational motion) stars. The thickness of the disk is set such that Toomre’s stability criterion is satisfied. The radial velocity dispersion  corresponds to the value of the Solar neighborhood, i.e. . The Gaussian-distributed random velocity dispersions are equal to ,  and .

Kuzmin’s and exponential model

We will use a program for initial conditions generation to create a disk galaxy. A model realization is in the system of units in which the gravitational constant , the total mass , and the scale radius , are all equal to unity. We will evolve initial conditions with the Barnes-Hut N-body simulation code with the opening angle  and the softening length .

Table 5–3: Parameters of the Kuzmin’s and exponential model simulations.

The timestep was set equal to  and the overall simulation covers .


Figure 5–19: An initial setup for the hot disk with the velocity dispersion added to circular velocities.





205 million years

409 million years




614 million years

818 million years

1.02 billion years

Figure 5–20: A time sequence for the hot exponential disk with  in the xy plane evolved with the Barnes-Hut N-body algorithm.



205 million years

409 million years

614 million years

818 million years

1.02 billion years


Figure 5–21: A time sequence for the hot exponential disk with  in the yz plane evolved with the Barnes-Hut N-body algorithm. The thickening caused by the random motion of stars can be seen.


Figure 5–22: The rotation curve of the exponential disk.




205 million years

409 million years




614 million years

818 million years

1.02 billion years


Figure 5–23: A time sequence for the hot Kuzmin’s disk with  in the xy plane evolved with the Barnes-Hut N-body algorithm.



205 million years

409 million years

614 million years

818 million years

1.02 billion years


Figure 5–24: A time sequence for the hot Kuzmin’s disk with  in the yz plane evolved with the Barnes-Hut N-body algorithm. The thickening caused by the random motion of stars can be seen.


Figure 5–25: The rotation curve of Kuzmin’s disk.

First, we will check the stability of the model. The whole simulation is showing a continuous evolution and changes of both disks. After the formation of a strong axisymmetric instability, the system relaxes to a steady state. After relaxation, the system forgets its initial conditions.

Large-scale deviations from the axial symmetry can be recognized in the simulations – the rise of a strong central bar and spiral arms. Initial conditions of the disk galaxy are developing into the bar and two or four spiral arms. Early N-body researchers (e.g. Hohl, 1971) already recognized the violation of axial symmetry in disks.

With a continuous evolution of the N-body system, the bar becomes weaker, but seems to last forever (exceeding the age of the universe) (Sellwood, 1981). The variation of mass breaks the axisymmetry and gravitation amplifies it by creating the spiral arms and bars. It is inevitable that a spiral galaxy will develop the axisymmetric instability at some point in its evolution.

5.13   Composite models

Galaxies can be better approximated with more than one component. Galaxies are composed of various parts that can be viewed in isolation as spherical or disk distributions. When we add these components together, we obtain a composite model of galaxy.

It is interesting that  does not have to be related to  through Poisson’s equation. Multi-component model might be constructed, with components co-existing in the mutual gravitational potential  (disk, bulge, halo). Our numerical spiral galaxy will be such a composite model of a disk, spherical bulge and halo. For each component of the combined model, the distribution function can be found separately.

Initial conditions are generated from the isotropic distribution functions , , and , which represents the bulge, disk, and halo, respectively. We assume that the total potential is spherical even when the potential is aspherical because of the disk. Influence to current investigations can be neglected.

More techniques for generating initial conditions were introduced. Barnes (1988) constructed isolated components and then slowly induced them into the presence of each other. Hernquist (1993) used local Maxwellian approximation. Kuijken and Dubinski (1995) employed a technique involving spherical harmonic expansion of the potential. Boily et al. (2001) used a method without the need of computing anisotropic velocity dispersions.

Galaxy experiments

Our objective will be a simulation of the composite model of the Milky Way galaxy. The most intuitive is to start with a total disk mass that is equal to the mass of 400 billion stars of solar mass. A bulge mass will be 10 % of the disk mass and stellar halo will contain 30 % of the disk mass (see Table 5–4). We will evolve initial conditions with the Barnes-Hut N-body simulation code with the opening angle  and the softening length . Model galaxy will be generated with the following units


Table 5‑4: The Milky Way galaxy model G1.


In this system of units is the timestep equal to  and the overall simulation with 150,000 timesteps covers more than .


1 billion years

2 billion years




3 billion years

4 billion years

5 billion years

Figure 5–26: The evolution of the Milky Way galaxy G1 composite model in the xy plane.



1 billion years

2 billion years




3 billion years

4 billion years

5 billion years

Figure 5–27: The evolution of the Milky Way galaxy G1 composite model in yz plane.

Milky Way with dark matter halo

In the Milky Way galaxy model G2, we will add the isothermal dark-matter halo component. Our galaxy is probably contained in the halo of dark matter particles with radius nearly 250 kpc that is about 20 times more massive than visible parts of Galaxy, leading to luminous-to-dark ratios of 1:20.

Table 5‑5: Milky Way galaxy model G2.

Simulation parameters and units are the same as in the simulation of G1 model. The timestep is equal to  and the overall simulation with 150,000 timesteps covers more than .

5.14   Concluding remarks

In Chapter 5, we have shown how to create a computer model of a galaxy in order to study galaxy dynamics. We found that the construction of galaxy in a controllable way is difficult. Due to immense complexity, all models are very artificial in comparison to real galaxies. In spite of that, the initial density distribution function of models are in good agreement with observations of real galaxies. The galaxy models created were single or multicomponent systems in stable dynamical equilibrium.

Initial conditions were generated as follows. Specifying the mass density distribution function, we first calculate the model’s cumulative mass distribution function and corresponding gravitational potential. Then the mass density distribution function is expressed as a function of gravitational potential. The phase-space distribution function is calculated on the fly using a numerical formulation of Eddington’s formula. Once the phase-space distribution function has been calculated, one can start to randomly sample particles from the distribution function. If the use of another mass density profile is requested, all that is necessary, is to override a virtual function “rho” according to the chosen mass density profile. Through this approach, many kinds of models may be constructed. Models created in this way are quickly getting into equilibrium.

We created realizations of an elliptical galaxy from various spherical models. The spherical models were in equilibrium from the beginning. We have created disk models that showed continuous evolution. We saw that dynamically cold disk without a dark halo spontaneously formed features resembling galactic bar and spiral arms. It has been shown that a self-gravitating disk system is unstable unless a certain velocity dispersion and dark halo were included. All models were evolved for more than 1 billion years and movies from all computations were produced.