Galaxy Models

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

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

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 (

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

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.

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.

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

We will encounter Hamilton’s operator

(pronounced “nabla” or “del”). It is defined as

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

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

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

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

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.

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

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

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

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

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

To have a body on a stable orbit, its kinetic energy

must be exactly the same as a potential energy

, i.e.,

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

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 velocities in the range

about

, i.e. between

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.

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

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

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

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.

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

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


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


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

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.

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

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.

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


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.

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

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

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

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.

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

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.

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

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

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 Barnes-Hut N-body algorithm.

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.


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 Barnes-Hut N-body algorithm.

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.

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.

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.

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

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.

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.

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

5 Galaxy Models

5
Galaxy Models