“Animals do not care about the evolution of the universe. Nor do many humans.”
— Stephen Hawking —
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.
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 lowdensity 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.
Elliptical galaxies are not flat ellipses, but they are threedimensionally elliptical (watermelon 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 supermassive 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).
Disk
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 faceon view of a disk galaxy with spiral arms, bar and bulge.
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).
Bar
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 edgeon.
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

(5.1) 
can be rewritten as
. 
(5.2) 
We will encounter Hamilton’s operator (pronounced “nabla” or “del”). It is defined as
. 
(5.3) 
We will also encounter Laplace’s operator . It is defined as
. 
(5.4) 
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

(5.5) 
The configuration of N point masses (the Nbody 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
. 
(5.6) 
This quantity is called the intensity of gravitational field.
With continuous functions, we may use Poisson’s equation (Greiner, 2004)

(5.7) 
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 selfgravitating 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 selfgravity model, we must find the global gravitational potential, in which the mass density of all parts of a galaxy is combined as
. 
(5.8) 
5.4 Initial conditions for spherical systems
As was described in Chapter 4.1, initial conditions for the Nbody 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 onedimensional 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.9) 
A uniform sphere is one of the simplest approximations of the spherical system with the uniform density distribution function (volume density) given by
(5.10) 
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).
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
(5.11) 
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. Cutoff radius must be imposed for
a finite total mass.
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
(5.12) 
Another more realistic model is Hernquist’s (1990) model given by
(5.13) 
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
, 
(5.14) 
where parameter p is a positive number. The case is Hernquist’s model and the case Plummer’s model.
Another realistic model is Jaffe’s (1983) model expressed as
(5.15) 
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 roundoff errors when we will compute with numbers of tens orders different. Therefore, we will use nondimensional 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
. 
(5.16) 
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
. 
(5.17) 
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
(5.18) 
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 acceptancerejection 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 onedimensional. 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.19) 
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

(5.20) 
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

(5.21) 
The body will stay in the system for .
Phasespace distribution function
An equilibrium spherical system depends on the phasespace coordinates (LyndenBell, 1962). The structure of a galaxy is defined by a distribution function in the imaginary space called phasespace with 6 dimensions (three spatial coordinates for positions and three corresponding velocity components) at each moment (time) t. The phasespace 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 phasespace 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 subspace of 6D phasespace) as was described in Chapter 5.6. Now, we must find the range of velocities in the phasespace 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 phasespace 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]
(5.22) 
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
. 
(5.23) 
The cumulative mass distribution relates to the gravitational potential through Poisson’s equation as follows
(5.24) 
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 Nbody 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
. 
(5.25) 
We can also define a characteristic length (Meylan and Heggie, 1996)
. 
(5.26) 
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
, 
(5.27) 
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.
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 BarnesHut Nbody 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 realworld 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 .



Start 
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.



Start 
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.



Start 
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.



Start 
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
(5.28) 
Stars in the disk have only circular velocities (so called cold disk, where no random motions are induced). Disk rotation is generated to be counterclockwise in the xy plane (right handed Cartesian coordinate system is completed with zaxis and disk rotation is chosen to by given by the right hand rule with a thumb along zaxis).
A reasonable approximation for the disk galaxy is Kuzmin’s disk. Kuzmin found in 1956 the potential of an infinitesimally razorthin 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
(5.29) 
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).
Bodies in an exponential disk are distributed in the plane with the surface mass density distribution given by
(5.30) 
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
. 
(5.31) 
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 supermassive black hole. A model realization is in the dimensionless system of units. We will evolve initial conditions with the BarnesHut Nbody 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 .



Start 
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 fallingoff 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

(5.32) 
where is the orbital velocity of a star at the distance from the center of the galaxy. The cumulative mass inside radius is
. 
(5.33) 
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 masstolight 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.
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
, 
(5.34) 
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
. 
(5.35) 
Vertical velocity dispersion (dispersion in a direction, which is perpendicular to the disk plane) is chosen as the part of tangential velocity dispersion as
. 
(5.36) 
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 Gaussiandistributed 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 BarnesHut Nbody 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.



Start 
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 BarnesHut Nbody algorithm.

Start 

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 BarnesHut Nbody algorithm. The thickening caused by the random motion of stars can be seen.
Figure 5–22: The rotation curve of the exponential disk.



Start 
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 BarnesHut Nbody algorithm.

Start 

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 BarnesHut Nbody 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.
Largescale 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 Nbody researchers (e.g. Hohl, 1971) already recognized the violation of axial symmetry in disks.
With a continuous evolution of the Nbody 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. Multicomponent model might be constructed, with components coexisting 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 BarnesHut Nbody 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 .



Start 
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.



Start 
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 darkmatter 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 luminoustodark 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 .
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 phasespace distribution function is calculated on the fly using a numerical formulation of Eddington’s formula. Once the phasespace 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 selfgravitating 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.