.

.

.

Inhaltsverzeichnis

Foreword

The integration of differential equations and the numerical solution of integrals has gained much attention since the flowering of computer technology in the natural sciences. This book, which is called “Simply Integrate!’, contains descriptions of a collection of different integrators such as Taylor methods, Leapfrog methods, Hermite schemes, 2nd, 4th and 8th-order Runge-Kutta methods, or fancy ones, like a time-transformed 8th-order composition scheme where the highest accuracy is needed. An example are highly unstable periodic solutions of the three-body problem. In extreme cases one uses a “shooting method” from both sides to obtain closed orbits. A measure for the degree of instability is the dominant eigenvalue of the monodromy matrix in Floquet theory, which is not covered in this book, but for example in the book by Contopoulos (2002). The biggest part of the present book is concerned with “science that is hoary with age” (Platon, Timaios) but parts of it are also concerned with novel inventions and methods of modern physics and astrophysics.

The first main chapter of the present book is dedicated to the few-body problem. Starting with an introduction to the two-body problem, which is analytically solvable the author presents initial conditions for three- and four-body problems. These can be solved with the hereafter presented integration methods mentioned above. If heuristically means “concerning the inventions”, one may say: The chapter proceed heuristically, from simple to complex algorithms, from the two-body problem to the few-body problem.

The second main chapter of this book is about galactic dynamics and orbital mechanics in analytic gravitational potentials. At first spherically symmetric models of stellar systems are discussed and later multi compenent models of galaxies. The aim is to integrate orbits in these potentials with a self-written integrator, against the “ignava ratio”.

The third main chapter of this book is dedicated to linear stability analysis. The linear stability analysis is presented with many examples. Every example uses a general classification of equilibrium points, which is presented in the beginning. In this sense, the individual case is explained by summing it under a general classification or mathematical law. This is the paradigmatic method of the nomothetical sciences (Hempel, 1977; Gadamer, 1990).

Further chapters are concerned with the stellar initial mass function and astrophysical Poisson equations.

The field of work of the author of this book is “experimental stellar dynamics”. The word “experimental” means in the context of astrophysics, that computer experiments are being used to study the dynamical evolution of stellar systems, where stars interact through gravity. Since the author’s N-body simulations require a huge amount of computer power, he has made use of Beowulf PC clusters as well as the largest supercomputers in Europe, field programmable gate arrays (FPGAs), graphics processing units (GPUs) and the grid technology. For the plots in this book, GNUPLOT, IDL, R, OMNIGRAFFLE and GIMP have been used. A website of this book with more material can be found at the internet address http://www.simplyintegrate.de. On the one hand, the readers of this book should get an impression what it means to do computational physics, on the other hand, they should see how much fun it is!

And the scientist and truth seeker, who has once read the book Burning Daylight by Jack London, knows well about the comparableness between the scientific research and the utterly exhausting search for gold nuggets at the Klondike River. It remains, to exclaim from time to time with Jack London’s hero: “Burning daylight!”

1 The few-body problem

1.1 Introduction to the two-body problem

We follow the lecture notes of Kühn & Spurzem (2005) first. A similar approach to the problem can be found in the book by Binney & Tremaine (2008). Newton’s equation of motion for the relative motion of two bodies under their mutual gravitational force is

(1.1)

where G is the gravitational constant, M = m₁ + m₂ is the total mass, r = (x₂ - x₁,y₂ - y₁,z₂ - z₁) the relative vector, r = |r| the relative distance and the two dots denote a second derivative with respect to time t in the notation of Newton. The bound solutions of Eqn. (1.1) are ellipses, the general solutions are conic sections. Kühn & Spurzem (2005) note that the Hamiltonian of the relative motion leads to μ = -Gm₁m₂r∕r³ with the reduced mass μ = m₁m₂∕M and the above equation is that divided by μ. Therefore we can recover the sum M = 4π²a³∕(GT²) of the two masses, e. g. of a double star, using Kepler’s third law from the knowledge of the relative motion (semimajor axis a and orbital period T of the relative motion), provided that we can deproject the orbit and learn to know about its inclination with respect to the celestial sphere (see Figure 1, left-hand side). To obtain the individual masses, the knowledge of the ratio m₁∕m₂ is needed in addition. It can be obtained from the knowledge of the center of mass and that of the motion of the double star components around it (see Figure 1, right-hand side).

Numerical solution: For a numerical integration we need to replace Eqn. (1.1) by two coupled first-order differential equations,

Note that these cannot be solved independently. A dimensionless formulation can be introduced (e. g. Kühn & Spurzem, 2005). We define

Using this we obtain

where the prime denotes a derivative with respect to the scaled time τ in the notation of Lagrange.

The simplest Taylor integrator, the Euler integrator, can be written as

Euler-Methode / Euler method:

where a = v′ = dv∕dτ is the dimensionless acceleration, Δt is the time step and the subscripts ‘0’ (and ‘1’) stand for the current (next) time step. The local error (in one time step) of the Euler method is of order a₀(Δt)² in position and ₀(Δt)² in velocity.

Although it is possible to implement the Euler scheme straightforwardly in the form of a computer program, the author recommends to begin the integration exercises with the initial conditions for three- and four-body problems given in the next section. The reason is that the computer program should be written generally enough, such that a few-body problem with N > 2 can be solved, where the variable N stands for the particle number.

Integrals of motion: We prove the temporal constancy of the integrals of the relative motion,

where r = |r| and v = |v|. We use the Graßmann identity

We use Eqn. (1.14) to calculate

(1.15)

and find

Analytical solution: For an analytical solution of Eqn. (1.1) we consider the expression

(1.19)

This can be written as

which is the well-known solution of a conic section with the orbital plane (e. g. Kühn & Spurzem, 2005) with the phase angle ψ = ψ(t).

Hamiltonian formulation: The conserved Hamiltonian is given by

(1.22)

where T and V are kinetic and potential energy, respectively, and H₀ a constant. With the dimensionless expressions for the relative motion of the two-body problem

(1.23)

we obtain

where we have used the separability property of the Hamiltonian, according to which T = T(v) depends only on v and V = V (r) depends only on r. The notations of the gradient operators are ∇_r  ⋅ =   ⋅ and ∇_v  ⋅ = ⋅, where ⋅ denotes a placeholder for a scalar quantity. We recognize the equations of motion (1.2) and (1.3) of the two-body problem.

Lagrangian formulation: The Lagrange function is given by

(1.26)

in dimensionless phase space coordinates. The Euler-Lagrange equations of motion are given by

(1.27)

or

(1.28)

We recognize the dimensionless formulation of Eqn. (1.1).

1.2 Initial conditions for three- and four-body problems

The initial conditions are given in the form m x y z v_x v_y v_z.

Acht / Eight

Pytagoreisches Problem / Pythagorean problem

Herz / Heart

Schwimmer / Swimmer

Ducati / Ducati

Stern / Star

Special solutions of the three- or four-body problem are described in:

Szebehely & Peters (1976); Hénon (1976); Moore (1993); Simó (2000); Chencincer & Montgomery (2000); Broucke (2002); Vanderbei (2004); Martynova & Orlov (2009); Šuvakov & Dmitrašinović (2013); Ouyang & Xie (2013).

After a discussion of the unit scaling in the next section many different integration algorithms are described in following sections. The author recommends, after writing your own integration program according to one of these algorithms, to begin the integration exercises with the figure “eight”, since it is easy to integrate due to its stability and the fact, that close encounters do not occur at all. Afterwards, you may try to solve the Pythagorean problem. The Pythagorean problem is difficult to integrate. A very accurate integrator must be used which is able to cope with the numerous close encounters.

In addition to the initial conditions for the figure “heart” mentioned above we also state the initial conditions for two other choreographic solutions (Simó, 2000):

Choreographie 1 / choreography 1

Choreographie 2 / choreography 2

Particularly the second solution is highly unstable. A “shooting method” from both sides in combination with a highly accurate integrator is recommended to obtain a closed solution. For this purpose, the sign of the velocity vectors must be reversed. Then one integrates (“shoots”) from both sides, once with the standard initial conditions and once with the initial conditions with reversed signs of the velocity vectors. If the accuracy is high enough, one may hope that the solutions meet each other in the middle after half an orbital period. For further choreographies of the three-body problem we refer to the works of Carles Simó.

1.3 Few-body units

We follow the procedure from Sverre Aarseth’s direct N-body program NBODY6 (Aarseth, 1999, 2003), but extend the method for the purpose of scaling to arbitrary units. We assume that the N-body system has units in which the gravitational constant has the value G the total mass has the value M and the total energy has the value E. We would like to scale to units in which these quantities have the values G′, M′ and E′.

For example, we might be interested in scaling to standard N-body units according to Heggie & Mathieu (1986), which are defined by G′ = M′ = -4E′ = 1.

Or, for example, we might want to scale a double or triple star system to physical units, in which G′ = (222.3)^-1 pc³M_⊙^-1Myr^-2, M′ = 1 M_⊙ and E′ = 1 M_⊙pc²Myr^-2.

We must proceed as follows:

Calculation of the total mass: We calculate the total mass firstly,

(1.29)

where N is the number of particles.

Translation into the center of mass system: We calculate the center of mass and its velocity,

If necessary (x_cm ⁄= 0 or v_cm ⁄= 0), we apply a translation to the center of mass frame,

where the prime denotes the quantities in the center of mass frame.

Calculation of the energy: Now we calculate the total energy in the center of mass frame,

(1.35)

where T and V are kinetic and potential energy. It can be seen that the positions and velocities are coupled through the masses, which enter quadratically in the potential energy, and the coupling constant G.

Scaling the masses, positions and velocities: We define three scaling factors

(1.36)

We may set

(1.37)

for standard N-body units (Heggie & Mathieu, 1986). Then we scale the individual masses, positions and velocities as follows:

That’s it. The system has now units in which the gravitational constant has the value G′, the total energy is E′ and the total mass is M′.

1.4 Taylor integrators

Taylor integrators are based on the evaluation of the Taylor series

with Δt = (t-t₀), where x₁,v₁, x₀ and v₁ are the positions and velocities at the next and current time. The simplest Taylor integrator is the Euler integrator which is first-order accurate.

The acceleration and its time derivatives for the Newtonian gravitational force are given by

where x_ij = x_j - x_i and x_ij = |x_ij| are the relative position of particles j and i and its modulus. Note that for a_i⁽²⁾ and a_i⁽³⁾ more than one summation have to be carried out which is numerically expensive.

1.5 Leapfrog/Verlet integrators

1.5.1 Standard Leapfrog methods

Hamiltonian flow: As a first “tool” we will introduce the Hamitonian flow. The Poisson brackets are given by

(1.47)

where X = X(x,v,t) and Y = Y (x,v,t) are two “observables”, i.e. differentiable functions of the phase space coordinates x,v and time t (e. g. Kuypers, 1997). We define operators with a hat as

(1.48)

where ⋅ denotes a placeholder. Furthermore, we define the 6-component phase space vector as

(1.49)

We note that is the time evolution “operator”. The time evolution of each component of Q is governed by

(1.50)

The above equations (1.50) are called “flow equations” for the phase space coordinate Q_i. They are equivalent to the Hamiltonian equations of motion, which are written in standard notation as

(1.51)

The solution of the flow equations (1.50) can be written formally as

(1.52)

where exp is called the “Hamiltonian flow”.

Baker-Campbell-Hausdorff identity: As a second “tool”, we introduce the Baker-Campbell-Hausdorff (BCH) identity: The product of the exponentials of any non-commutative operators X and Y can be expanded as

(1.53)

with

(1.54)

to second order. It we apply the BCH identity twice we can find

(1.55)

with

(1.56)

which is also second-order accurate.

Leapfrog schemes: The Leapfrog methods are based on the splitting of the Hamiltonian,

(1.57)

where the kinetic energy T = T(v) depends only on the velocities and the potential energy V = V (x) depends only on positions. Let us call this property the “separability property”.

We can write

where the last approximate equality is valid because of the expansion (1.56).

The Hamiltonian equations of motion (1.50) can be written second-order accurate as

(1.59)

or alternatively as

(1.60)

Switching to the standard notation we can obtain from (1.60)

(1.61)

Therefore the integration can be performed in three time-ordered steps. Here we need to consider the time-ordering of the operators with a hat which is dictated by the fact that they are non-commutative. Note that the time ordering can also be written formally by applying the Wick-Dyson time ordering operator . A typical non-commutative operation from every-day life is, for example, washing and drying clothes.

Algorithms: The DKD (“Drift-Kick-Drift”) Leapfrog, also called “Position Verlet” is given by (e. g. Sutmann, 2006)

DKD Leapfrog:

1. “Drift”:

2. “Kick”:

The KDK (“Kick-Drift-Kick”) Leapfrog, also called “Velocity Verlet” is given by (e. g. Sutmann, 2006)

KDK Leapfrog:

1. “Kick”:

2. “Drift”:

Properties of leapfrogs: The leapfrog methods have the following good properties (c. f. McLachlan & Quispel, 2001):

1.

Symplecticity: The leapfrog is symplectic, i.e. it respects the structure and properties of Hamiltonian systems. For example, it does not show any secular drift in energy error.

2.

Symmetry conservation: The leapfrog ist time-symmetric or symmetric under time reversal. Trajactories, which are integrated forwards or backwards in time, respectively, are exactly identical. This property also implies that there is no secular drift in the energy error.

3.

Volume conservation: The phase space volume is conserved. For example, a circular orbit stays circular or the amplitude of a rosette orbit as is shown in Figure 9 stays constant during a long-term integration.

4.

Integral conservation: The leapfrog conserves angular momentum exactly, as a simple calculation shows (noticed by Patrick Glaschke):

We note that in practical computations every numerical scheme is prone to roundoff errors. For the connection to symplectic integration we refer to the work of Skeel (1999).

1.5.2 Logarithmic Hamiltonian method (LogH)

The LogH method, a leapfrog with adaptive time steps is in detail described in Preto & Tremaine (1999). It is based on a so-called “functional Hamiltonian” (cf. Mikkola & Tanikawa, 1999; Preto & Tremaine, 1999; Mikkola, 2008).

We introduce now a logarithmic Hamiltonian (LogH)

(1.75)

where 0 < ϵ ≪ 1 is an accuracy parameter and P₀ is an integral of motion which equals minus the initial energy. We remark that a different functional form for the Hamiltonian can be used as well. However, it should obey the separability property with respect to x and v.

In the special case of the relative motion of the two-body problem, which was treated in chapter 1, we have

(1.76)

with the dimensionless relative velocity v = |v| and the dimensionless distance r = |r|.

We obtain the equations of motion for the logarithmic Hamiltonian method,

where τ is a transformed time. The time step function is given by

(1.79)

Preto & Tremaine (1999) remark already the nice property that the Equations (1.77), (1.78) and (1.79) do not contain any square root operations which are computationally expensive.

The DKD (“Drift-Kick-Drift”) leapfrog algorithm can, according to Preto & Tremaine (1999) and Mikkola (2008) be written as follows with an accuracy constant ϵ and an integral of motion P₀ which is equal to minus the initial energy:

LogH Methode / LogH method:

1. “Drift”:

2. “Kick”:

Because of the Hamiltonian formulation the LogH method is going to preserve the two-body trajectory except for a phase error: If we integrate a two-body problem with the LogH method, the numerical solution will always stay on the corresponding conic section. The same holds for the TTL method which is discussed in the next section.

1.5.3 Time-transformed leapfrog (TTL)

The TTL method, also a leapfrog with adaptive time steps is in detail described in Mikkola & Aarseth (2002) and the Cambody lecture notes by Mikkola (2008).

The Hamiltonian equations of motion are given by

where x, v and a are position, velocity and force on a unit mass and the dot denotes a derivative with respect to time. Note that these are coupled and cannot be solved independently of each other. A time transformation

(1.88)

can be introduced, where Ω(x) > 0 can be the absolute value of the potential energy in the case of a few-body problem, or its square, or cube, or, for example,

(1.89)

We introduce an auxiliary variable W = Ω. The equations of motion in the transformed time become

where the prime denotes a derivative with respect to the transformed time. The trick is to calculate W from the differential equation

(1.93)

With this trick the leapfrog scheme can be applied, since Ω = Ω(x) similar to a = a(x) and the separability property can be respected. Switching to the transformed time we obtain

(1.94)

We obtain the differential equations

(1.95)

These can be solved, for example with a DKD (“Drift-Kick-Drift”) leapfrog method. We note, however, that the separability property is not exactly realized in the equations (1.95).

The algorithm can be written as follows: TTL-Methode / TTL method:

1. “Drift”:

2. “Kick”:

1.5.4 Higher-oder composition schemes

The composition schemes are higher-order leapfrog methods and in detail described in the review by Sutmann (2006). Valid is, for example, the

Theorem of Yoshida, Qin and Zhu: Let ∏ _i exp be a scheme of order 2n with the time step Δt. Then the time-ordered composition

(1.102)

with

(1.103)

is a scheme of order 2n + 2 (Yoshida, 1990; Qin & Zhu, 1992), where is the Wick-Dyson time ordering operator. We refer to McLachlan & Quispel (2001) for the proof.

Many other composition possibilities exist. Mathematicians try to minimize the number of computationally expensive force evaluations for a given order of the integrator. For example using the coefficients given in Table 1, Table 2 or Table 3, 4th-order, 6th-order or 8th-order composition schemes can be written in a “Position Verlet” fashion as a loop over n = 9, n = 7 or n = 15 coefficients in n = 9, n = 7 or n = 15 steps,

FOR i=1, n DO BEGIN

1.6 Hermite integrators

1.6.1 Standard Hermite method

The standard 4th-order Hermite scheme (Makino, 1991; Makino & Aarseth, 1992) is in detail described in Glaschke (2006, Section 4.4) or in Ernst (2009, Section 8.2.1). It is a predictor-corrector method, which is particularly used for direct N-body simulations with N ≫ 1, since it combines the advantages of being a fast integrator and a higher-order method. We cite the latter of the above descriptions in the following subsection:

Prediction step: “Positions and velocities of the ith particle are predicted according to

where the superscript ‘p’ stands for ‘predicted’ and the subscripts ‘0’ and ‘1’ for the current time t₀ and the following time t₁, respectively. The acceleration a_i,0 and jerk _i,0 have been calculated from positions x_i,0 and velocities v_i,0 at time t₀ according to the analytical formulas

where x_ji = x_j - x_i and v_ji = v_j - v_i are relative positions and velocities between the ith and jth particles, respectively. Acceleration a_i,1^p and jerk _i,1^p at time t₁ can be calculated as well from the predicted positions (1.107) and velocities (1.108) using the analytical formulas (1.109) and (1.110).

Interpolation of higher derivatives: A second way to write down a_i,1^p and _i,1^p approximately is by a Taylor series expansion which contains higher-order derivatives of acceleration:

Since we know the quantities Δt, a_i,0, _i,0, a_i,1^p and _i,1^p, we can solve the system of equations (1.111) – (1.112) to obtain the Hermite interpolation of the higher-order derivatives of acceleration:

Correction step: With these quantities, positions and velocities from the prediction step can be corrected to higher orders,

where the superscript ‘c’ stands for ‘corrected’.” (Ernst, 2009, Section 8.2.1)

1.6.2 Iterated time-symmetric Hermite method

An extension is the iterated time-symmetric Hermite method (Kokubo, Yoshinaga & Makino, 1998)

Prediction step: The prediction is given by

Correction step: The corrector step is given by

where the Hermite interpolations in Eqns. (1.113) and (1.114) have been used implicitly.

Iteration: The iteration of the corrector is achieved by returning to Eqns. (1.109)–(1.110) with the predicted positions replaced by the more accurate values x^c,v^c (Glaschke, 2006, Abschnitt 4.5).

1.6.3 Higher-order Hermite methods

Hermite mothods of higher order are explicitely described by Nitadori & Makino (2008).

Prediction step: For example, the predictor for a 6th-order Hermite method is given by

Correction step: The corrector for the 6th-order Hermite method is given by

An eighth-order method is also described in Nitadori & Makino (2008).

1.6.4 Adaptive time steps for the Hermite scheme

The standard 4th-order Hermite scheme is based on the Taylor expansions

for the ith particle, where the subscripts ‘0’ and ‘1’ stand for the current time t₀ and the following time t₁, respectively. The idea of adaptive time steps is as follows: Depending on the order of the integrator, its local numerical error depends usually on some positive integer power of the step width. The local error should decrease monotonically with the step width. An adaptive time step criterion adjusts the time step in such a way, that the local error remains constant or, at least, has an upper bound during the integration. Time step criteria therefore involve an estimate of the local error.

Time step criteria: One may summarize a list of restrictions for the choice of an adaptive time step for the standard 4th-order Hermite integrator as follows:

1.

The time step must be a combination of dimensionless parameters and the quantities involved in (1.125) and (1.126),

2.

The time step must have the dimension of time and transform like a time under rescaling of units.

3.

The time step must be invariant under Galilei transformations of the N-body system.

4.

The time step must be well defined for special cases.

Aarseth criterion: It turns out that the Aarseth criterion

(1.127)

fulfills all requirements, where η is a dimensionless accuracy parameter. This criterion is well defined for the special case a_i = 0 at the Lagrangian points a tidal field or the “cold collapse” with _i = a_i⁽³⁾ = 0, where all particles are initially at rest (Aarseth, 1985).

Hierarchical block time steps: Usually one quantizes the time steps for N ≫ 1 in powers of 1∕2 such that blocks of particles with equal time steps can be integrated forwards simultaneously.

1.7 Runge-Kutta integrators (RK)

The Runge-Kutta method is in detail described in Press et al. (1986). As a standard numerical integrator, it “takes an initial condition and moves the objects in the direction specified by the differential equations”.

Suppose we have a system of N particles, which interact through gravity. Let m_i be the mass of the ith particle and Q_i,0 a 6-dimensional state vector, which contains position and velocity of the ith particle at time t₀. We have then

(1.128)

where the dot denotes the time derivative. The acceleration a_i,0 obeys the Newtonian law of gravitation,

(1.129)

where x_ij,0 = x_j,0 - x_i,0,r_ij,0 = |x_ij,0| and the summation is over all those pairs of particles, where one of the particles is the ith particle. Of course, the ith particle does not interact with itself.

Since a_i,0 depends on m₁,...,m_i-1,m_i+1,...,m_N and r_1,0,...,r_N,0, we have a system of ordinary differential equations

(1.130)

where F is the known vector mapping. Given the masses m_i and the vectors Q_i,0, we want to find the vector Q_i,1 at time t₁ = t₀ + Δt, where Δt is the time step.

For the pth-order Runge-Kutta mathod, we have to compute the following p 6-dimensional vectors:

(1.131)

with coefficients a_jm given in tables 4, 5 and 6. In the above definition, F is the known vector mapping from (1.130). The integration step is then

(1.132)

with coefficients b_j given in tables 4, 5 and 6. The coefficients for the second- and fourth-order scheme are given in Press et al. (1986). It is not known to the present author who calculated the coefficients for the eighth-order method.

1.8 Tasks

1.

How many summations have to be carried out for the calculation of a⁽²⁾ and a⁽³⁾ in Eqns. (1.45) and (1.46)?

2.

Scale the initial conditions of the Pythagorean problem to standard N-body units (Heggie & Mathieu, 1986).

3.

Integrate the initial conditions given in Section 1.2 with an integrator of your choice. Visualize the solution graphically, e. g. as shown in Appendix C, and/or animate it with the method described in Appendix D.

4.

Measure the growth behaviour of the relative energy error of an integrator of your choice and compare with that shown in Figure 8 for the standard Hermite method.

5.

Measure the time and distance of the closest encounter in the Pythagorean three-body problem before the binary forms and write down which particles are involved.

6.

Rotate the initial conditions for the Ducati orbit with an affine transformation by 45 degrees in the xy plane and integrate it.

7.

Think about the question, whether the 6th-order Hermite scheme described in Nitadori & Makino (2008) is getting time-symmetric or not after a few iterations of the corrector.

8.

Develop an adaptive time step criterion as in (1.127) for the 6th-order Hermite scheme described in Nitadori & Makino (2008).

2 Galactic dynamics

The aim of this chapter about orbital mechanics and galactic dynamics is to enable the reader to integrate orbits in different analytic gravitational potentials. It is left to the reader to write a program for the integration of orbits in these gravitational potentials. It may be useful to keep the program as general as possible to simplify the inclusion of new gravitational potentials without larger problems.

2.1 Spherically symmetric models

2.1.1 Mechanical foundations

Coordinates: A spherically symmetric self-consistent model of as stellar system, which is static in time and a solution of both the Poisson and the Boltzmann equation, is given by gravitational potential Φ(r), cumulative mass M(r), density ρ(r) and velocity dispersion σ(r). These scalar quantities depend only on radius r, which is defined as the distance from the center, not from time t or velocity v. We define

(2.133)

Poisson Gleichung: The Poisson equation for a spherically symmetric system is given by

(2.134)

where ∇² is the Laplacian operator and G is the gravitational constant.

Physical quantities: The force on a unit mass (the specific force or acceleration) is given by

(2.135)

where ∇ is the gradient operator.

Starting from the gravitational potential, we have

where we have used in the first step Newton’s second theorem (see Binney & Tremaine, 1987, p. 34f.) together with Newton’s law of gravitation.

Starting from the density, we have

where we assume lim_r→∞Φ(r) = 0 in the last step.

The circular velocity v_c(r) and escape velocity v_e(r) are given by

The total potential energy W is given by the equivalent expressions

The derivation of the above integrals is discussed in chapter 2 of the book by Binney & Tremaine (1987). If the system is in virial equilibrium, the value of the potential energy W is connected to the values of total kinetic energy K = -W∕2 and the total energy E = -K = W∕2.

From Eqn. (2.144) follows the root mean squared escape speed,

(2.145)

for a finite total mass M.

If the velocity dispersion tensor is isotropic, we can obtain the velocity dispersion σ(r) by integrating the Jeans equation of hydrostatic equilibrium,

(2.146)

The velocity dispersion is given by the integral

(2.147)

Virial theorem: The kinetic energy due to random motion is given by the integral

(2.148)

where a an integration by parts has been carried out in the first step and the condition (2.146) of hydrostatic equilibrium has been plugged in in the second step. If Π is finite it follows that C = 0 and the virial theorem is fulfilled (Ernst, 2009).

Integrals of motion: Two integrals of motion in static spherically symmetric stellar systems are

2.1.2 Kepler potential

The Kepler potential is solution of the Poisson equation for a point mass in the origin of coordinates and given by

(2.151)

Applications include, e. g., planetary dynamics or stellar orbits around a supermassive black hole. The acceleration is given by

(2.152)

The jerk, which is of relevance for the Hermite integration scheme, is given by

Integral of motion: An additional integral of motion in the Kepler potential is

(2.154)

2.1.3 Power-law models

The spherically symmetric power-law or scale-free models can be used to study the dynamics in the centers of galaxies. Note that power laws are connected to self-similarity, i.e. they are scale-invariant. In particular, the models discussed below are invariant under spatial magnification. They have a potential of the form

where r is the radius and r₀ is a length unit. The acceleration has the form

(2.156)

The jerk is given by

2.1.4 Plummer model

The gravitational potential of a Plummer model is given by

(2.158)

where a is the Plummer radius. The acceleration is then given by

(2.159)

The jerk, which is of relevance for the Hermite integration scheme, is given by

(2.160)

2.1.5 Dehnen models

The gravitational potential of a Dehnen model (Dehnen, 1993; Tremaine, 1994) has, in contrast to that of the Plummer model, a “cusp” in the center. It is given by

(2.161)

where a is a scale radius. The acceleration is given by

(2.162)

The jerk which is of relevance for the Hermite integration scheme, is given by

Special cases are the Hernquist model with γ = 1 (Hernquist, 1990) and the Jaffe model with γ = 2 (Jaffe, 1983).

2.1.6 Logarithmic potential

A logarithmic gravitational potential, which is often used for dark matter halos of galaxies, is given by

(2.164)

The acceleration is given by

(2.165)

The jerk is given by

(2.166)

2.2 Multi component models

2.2.1 The 3PLK Milky Way model

The parameters of the 3-component Plummer-Kuzmin model (Miyamoto & Nagai, 1975) of the Milky Way are given in table 7.¹ The total gravitational potential is given by a a superposition of three potentials for bulge, disk and halo. The circular velocity as a function of galactocentric radius is given in table 8. The gravitational potential has the form

(2.167)

The acceleration for the Plummer-Kuzmin model is given by the expression

(2.168)

The jerk is given by the lengthy expressions

(2.169)

(2.170)

which also depend on the velocity components. A few orbits are plotted in Figures 10 and 11.

2.2.2 JSH95 Milky Way model

The gravitational potential of the JSH95 model of the Milky Way (see Johnston, Spergel & Hernquist, 1995; Dinescu, Girard & van Altena, 1999) is given by

(2.172)

with M_b = 3.4 × 10¹⁰M_⊙, c = 0.7 kpc, M_d = 10¹¹M_⊙, a_d = 6.5 kpc, b_d = 0.26 kpc, v₀ = 128 km s^-1, d = 12.0 kpc. The bulge has a Hernquist profile, the disk is represented as a Plummer-Kuzmin model and the halo by a logarithmic potential.

The acceleration for the JSH95 model is given by the expression

2.2.3 P90 Milky Way model

The gravitational potential of the P90 model of the Milky Way (see Pacyński, 1990; Dinescu, Girard & van Altena, 1999) is given by

with M_b = 1.12 × 10¹⁰M_⊙, a_b = 0.0 kpc, b_b = 0.277 kpc, M_d = 8.07¹⁰M_⊙, a_d = 3.7 kpc, b_d = 0.20 kpc, M_h = 5 × 10¹⁰M_⊙, d = 6.0 kpc.

The acceleration for the P90 model is given by the expression

2.2.4 Zotos’ barred spiral galaxy model

The simple 2D effective potential of a barred spiral galaxy according to Zotos (2012), which is visualised in Figure 12, is given by the terms

with the parameters α = 8, β = 1.3, b = 2, v₀ = 15, M_d = 9500, M_b = 3000, M_n = 400, c_b = 1.5, c_n = 0.25 and c_h = 8.5. We have Ω_b = 1.25. The model consists of a logarithmic halo, a disc, a Plummer nucleus (bulge) and a bar. The potential rotates clockwise at constant angular velocity Ω_b. The effective potential in Eqn. (2.176) is the cut through the equatorial plane of a typical 3D potential of a barred spiral galaxy. We remark that for a axisymmetric halo the parameter β should be set to β = 1.

In Eqn. (2.176), the model units of length L₀, velocity V ₀, angular velocity Ω₀, time T₀ and mass M₀ of the parameters, in which the gravitational constant G = 1, can be scaled to physical units of a realistic barred spiral galaxy as NGC 1300:

The equations of motion in the rotating frame are given by

(2.180)

where r = (x,y,z) is the position vector, the dot denotes a derivative with respect to time and Φ is given by the first four terms in Eqn. (2.176). The Jacobi energy is an integral of motion and given by

(2.181)

A critical Jacobi energy is given by the value of the effective potential Φ_eff at the position of the two saddle points in Figure 12. Figures 13 and 14 illustrate fully-fledged spiral arms, which have been formed in the potential of Eqn. (2.176), as described in Ernst & Peters (2014). The spiral arms form out of escaping particles with overcritical Jacobi energies.

The density of points along one orbit in figures 13 and 14 is taken to be proportional to the velocity of a particle. This method can be called the “velocity-dependent plotting technique”. The spiral galaxies in this theory display spiral patterns which are persistent in time, in contrast to the transient spiral patterns in the old density wave theory (e.g. Pasha, 2004a,b). The winding problem in differentially rotating disks does not occur (“Weizsäcker’s cow”). Also, a large number of different galaxy morphologies can be reproduced in this theory.

2.3 Tasks

1.

Write an orbit integration program for the integration of orbits of a point mass in the gravitational potentials which have been presented in this chapter with one of the integrators which have been described in the previous chapter. The system of first-order differential equations to be integrated is given by

where x, v and a are position, velocity and acceleration of the point mass.

• According to the interest of the reader, a composition scheme, a Hermite scheme or a Runge-Kutta scheme is recommended for this purpose.
• Adaptive times steps may be implemented later after a corresponding modification of the program.

2.

Calculate analytical expressions for the physical quantities M(r),ρ(r),v_c(r),σ(r),... for some models..

3.

Calculate the rotation curve of the 3PLK model.

4.

Calculate the jerk for the JSH95 and the P90 model.

5.

Calculate spiral arms in the effective potential 2.176 according to the method described in Ernst & Peters (2014), i.e. by making use of the velocity-dependent plotting technique.

3 Linear stability analysis

3.1 Introduction to linear stability analysis

The linear stability analysis of equilibria is an important part of nonlinear dynamics. For this purpose one linearizes the dynamics around the equilibrium and considers the dynamics in this linearized theory, which is simpler than the full nonlinear theory. Methods of linear algebra are used. A German reference for linear algebra are the two books by Lorenz (1992a,b). In particular, the eigenvalue problem for the Jacobi matrix must be solved.

One-dimensional systems: Let

(3.182)

be a one-dimensional (1D) map. The Jacobi matrix consists of only one element,

(3.183)

If the modulus of the derivative, is smaller than one (larger than one) the corresponding fixed point is stable (unstable). If it is equal to one, we speak of a marginally stable fixed point.

Two-dimensional systems: Let, in a two-dimensional (2D) system,

The Jacobi matrix is given by

(3.187)

We can classify equilibrium points in two dimensions (2D) as follows (e. g. Izhikevich, 2007):

• Node: Both eigenvalues are real and have the same sign. The node is stable (unstable) when all eigenvalues are negative (positive).
• Saddle: Both eigenvalues are real and of opposite sign.
• Focus (or “spiral point”): The two eigenvalues are complex-conjugate. The focus is stable (unstable) when the eigenvalues have negative (positive) real parts.

Three-dimensional systems: Let, in a three-dimensional (3D) system,

The Jacobi matrix is given by

(3.191)

We can classify equilibrium points in three dimensions (3D) as follows (e. g. Izhikevich, 2007):

• Node: All eigenvalues are real and have the same sign. The node is stable (unstable) when all eigenvalues are negative (positive).
• Saddle: All eigenvalues are real and at least one differs in sign from the others.

• Focus-Node: The point has one real eigenvalue and one pair of complex-conjugate eigenvalues, and all eigenvalues have real parts of the same sign. The equilibrium is stable (unstable) when the sign is negative (positive).
• Saddle-Focus: The point has one real eigenvalue and one pair of complex-conjugate eigenvalues, and the real eigenvalue differs in sign from the real parts of the complex-conjugate ones.

3.2 Logistic Map (1D)

History: Historically, the non-linear logistic equation

(3.192)

here with a free parameter λ, was introduced by Pierre-François Verhulst for studying population dynamics (Verhulst, 1838). The solutions for n(τ) < 1 are characteristic S-shaped curves, so-called sigmoid functions.

Logistische Abbildung: The associated discrete logistic map is given by

(3.193)

with a free parameter λ. A plot of x_n vs. λ for a large n yields the so-called Feigenbaum diagram shown in Figure 15. By increasing λ, it can be seen that at the period doubling bifurcations the fixed point(s) get(s) unstable. This is the “route into chaos”.

Linear stability analysis: The Jacobi matrix element is given by

(3.194)

For a fixed point we have the condition

(3.195)

We obtain two solutions,

(3.196)

We find

• J₀₀(x_n,1) = 4λ. x_n,1 is stable for 0 < λ < 0.25.
• J₀₀(x_n,2) = 2 - 4λ. For 0.25 < λ < 0.75 x_n,1 is unstable and x_n,2 is stable.
• For λ > 0.75 both x_n,1 and x_n,2 are unstable and we obtain a series which oscillates between two (and later more) fixed points x_n,3 = f(x_n,4) and x_n,4 = f(x_n,3) as can be seen in Figure 15. This is called a 2-periodic or period-2 orbit. For the period-2 orbit we have the condition x_n+2 = x_n in Eqn. (3.193).

3.3 Hénon-Heiles system (2D)

History: Historically, the Hénon-Heiles system was the first system for which the existence of a so-called adelphic or non-classical integral of motion was demonstrated numerically. Figure 16 shows that the potential is not axisymmetric but has only a triangular symmetry according to the dihedral group D₃ (Cushman, 2005). Hence the z-component of angular momentum is only approximately conserved. The energy is conserved exactly since the Hénon-Heiles system is a Hamiltonian system. The adelphic integral of the Henon-Heiles system can be represented as a power series, whose lowest order is the z-component of angular momentum (Finkler, Jones & Sowell, 1990).

The Hénon-Heiles potential (Hénon & Heiles, 1964) is given by

(3.197)

Figure 16 shows the equipotential lines. The potential has a minimum H₁ and three saddle points H₂, H₃ and H₄.

Equations of motion: The equations of motion are

Figure 17 shows two Poincaré surfaces of section at two different energies. In the top panel of Figure 17 all orbits are quasiperiodic or periodic. In the bottom panel of Figure 17 there are large dotted chaotic regions in which smaller quasiperiodic islands are embedded.

Integrals of motion: The integrals of motion are

where I₃ is the adelphic integral to third order (Finkler, Jones & Sowell, 1990).

Linear stability analysis: Equilibrium points are determined via

(3.204)

The solutions are

(3.205)

The central potential is Φ(H₁) = 0. The potential at the saddle points is E_crit = Φ(H₂) = Φ(H₃) = Φ(H₄) = 1∕6. For energies E > E_crit particles can escape from the central region of the potential.

We proceed with a linear stability analysis. The Jacobi matrix is given by

We consider the 2 × 2 submatrix of J,

The characteristic polynomial is given by

(3.208)

The general solutions are given by

(3.209)

We find:

1.

The eigenvalue at H₁ is λ₁ = λ₂ = -1. It is two-fold degenerate, real and negative. Therefore H₁ is a stable equilibrium point.

2.

The eigenvalues at H₂ are λ₁ = -3 and λ₂ = 1. Therefore H₂ is a saddle point.

3.

The eigenvalues at H₃ are λ₁ = -3 and λ₂ = 1. Therefore H₃ is a saddle point.

4.

The eigenvalues at H₄ are λ₁ = -3 and λ₂ = 1. Therefore H₃ is a saddle point.

Computation of Poincaré surfaces of section:

The Poincaré surfaces of section, as shown in Figure 17, are 2D cuts through the 4D phase space. For example, for an y - v_y surface of section we take an initial condition y = y₀,v_y = v_y,0,x = x₀ = 0 and calculate v_x from v_x = , where E is the total energy, Φ is the potential and v_x is chosen to be always, say, the positive root. Then we integrate the initial condition forwards in time and plot a dot at each consequent (=“piercing point through the surface of section”) with x = 0 when v_x ≥ 0.

Interpretation of Poincaré surfaces of section:

Dotted regions mark chaotic regions in the phase space, the “chaotic sea”. Embedded are usually regular islands of invariant curves, which are formed by quasiperiodic orbits. Quasiperiodic orbits respect, besides the energy, another integral of motion, for example the z-component of angular momentum or, as is the case for the Hénon-Heiles system, an adelphic integral of motion. In the center of each island lies a stable fixed point, which corresponds to a stable periodic orbit of period n. Between the islands one can sometimes find hyperbolic fixed points, which belong to unstable periodic orbits. As is well-known (Contopoulos, 2002), these are a source of chaotic dynamics. One may also read the famous mathematical paper by Li & Yorke (1975), Period three implies chaos.

Algorithm for finding periodic orbits on a Poincaré surface of section:

In this connection, we remark that the periodic orbits, shown as colored dots on the title page of this book, can be found with a self-written program as follows: Let r₀ be the approximate phase space diameter of the Poincaré surface of section. We define an n-periodic orbit as a periodic orbit of period n.

1.

Calculate for an n-periodic orbit the nth consequent (“= piercing point through the surface of section”) of an orbit for an initial condition in the quasiperiodic island on the Poincaré surface of section close to the periodic orbit and measure the distance Δr₀∕r₀ of the nth consequent from the initial condition.

2.

Calculate the nth consequents of four neighbouring orbits in the directions north, south, east, west, which are located at a distance 0 < Δr′∕r₀ ≪ 1 away from the original orbit and measure for an n-periodic orbit the distances Δr_i∕r₀, i = 1,...,4, of the nth consequents with the corresponding initial condition.

3.

Go one step of size 0 < Δr′∕r₀ ≪ 1 in the direction, which belongs to the smallest of the four distances if min(Δr_i∕r₀) < Δr₀∕r₀. If the latter condition is not fulfilled the periodic orbit has be found to an accuracy of order (Δr′∕r₀).

4.

Decrease Δr′∕r₀ by a suitable factor.

5.

Repeat the procedure from 1. until the n-periodic orbit has been found to a reasonable accuracy.

3.4 Lotka-Volterra system (2D)

Population dynamics: Assume a system, in which two species interact, the predator and the prey populations. In the field of astrophysics at the university, the predator species might be the astrophysicists and the prey “species” the public funds (not: third-party funds). The theory of such a system is well-known and can near the fixed points be described by a system of coupled differential equations, the Lotka-Volterra equations,

where the variables are defined as follows:

• N₁: the number of astrophysicists,
• N₂: the amount of public money for astrophysics,
• α > 0: the constant decreasing rate of the number of astrophysicists per unit of public money,
• β > 0: the constant growth rate of the amount of public money per astrophysicist,
• γ > 0: the constant growth rate of the number of astrophysicists in the limit N₂ = 0,
• δ > 0: the reduction rate of the amount of public money in the limit that N₁ = 0,

Dimensionless formulation: We may introduce scaled variables (e. g. Kühn & Spurzem, 2005),

to arrive at the dimensionless formulation

with only one free parameter λ. In logarithmic notation, Eqns. (3.213) and (3.214) can be written more concisely as

Solutions: Figure (18) shows phase portraits for the Eqns. (3.213) - (3.214). These are similar to the Poincaré surfaces of section except for the fact that the solutions of the Eqns. (3.213) - (3.214) are all regular (non-chaotic), such that the phase can be identified easily. One solution is shown in Figure 19 for λ = 0.5.

Linear stability analysis: For stationary solutions we set the left-hand sides of Eqns. (3.213) and (3.214) to zero,

The solutions for arbitrary λ are

(3.219)

For λ = 0 other solutions exist.

We proceed with a linear stability analysis. For that purpose we consider the Jacobi matrix

(3.220)

with m = 1,2 and n = 1,2. The eigenvalues are given by the zeros of the characteristic polynomial, i.e.

(3.221)

The general solutions are given by

(3.222)

The stability of the fixed point E₁ can be seen as follows:

At E₁ the solutions are

(3.223)

We have Re = 0 for λ > 0 and and therefore E₁ is a marginally stable fixed point for this parameter range. The normalized eigenvector basis at E₁ is given by

(3.224)

The top panel of Figure 18 shows that for λ = 0.5 these eigenvectors point indeed in the right direction at E₁, i.e. that they are parallel to the n₁ and n₂ axes, respectively.

3.5 Lorenz attractor (3D)

History: Edward N. Lorenz was a meteorologist working on the theory of the weather. “He reduced the Navier-Stokes equations in the Boussinesq approximation for a plane, from below heated fluid model [...] to a system of three ordinary nonlinear differential equations, which governs the temporal evolution of three fundamental modes, a velocity and two temperature distributions” (Argyris, Faust & Haase, 1995).

Equations of motion: The system of equations of motion for the Lorenz attractor (Lorenz, 1963) are given by

Solution: They describe a so-called “strange attractor” (Ruelle & Takens, 1971), which is shown in the upper panel of Figure 20. A possible choice of the parameters is σ = 10, b = 8∕3, r = 28.

Linear stability analysis: The Jacobi matrix is given by

The fixed points are for b ⁄= 0 given by

(3.229)

The Jacobi matrix at L₁ is given by (cf. Argyris, Faust & Haase, 1995)

The eigenvalues of J^L₁ are given by

(3.231)

We find that

• For 0 < r < 1 all three eigenvalues are negative. L₁ is a stable node.
• For r = 1 we find λ₁ = 0,λ₂ = -(σ + 1) < 0,λ₃ = -b < 0. L₁ is a marginally stable node.
• For r > 1 we find λ₁ > 0,λ₂ < 0,λ₃ < 0. L₁ is a saddle point.

3.6 Rössler attractor (3D)

History: Otto E. Rössler searched for a simpler “strange Attractor” than the Lorenz attractor. He found one. Equations of motion: The system of equations of motion for the Rössler attractor (Rössler, 1976) are given by

Solution: The lower panel of Figure 20 shows the Rössler attractor. A possible choice of the parameters is a = 0.15, b = 0.20, c = 10.

Linear stability analysis: The Jacobi matrix is given by

The fixed points are for a ⁄= 0 given by

(3.236)

The Jacobi matrix at R₁ is given by

The characteristic polynomial at R₁ is given by

(3.238)

The eigenvalues of J^R₁ are given by

(3.239)

We find that

• For 0 < a < 2 λ₁ and λ₂ are complex-conjugate with Re(λ₁) = Re(λ₂) > 0 and λ₃ = -c < 0 is real. R₁ is an unstable spiral point in the e₁∕e₂ plane and stable in the direction of e₃, where e_i is the eigenvector corresponding to λ_i. Such a point is called a “saddle-focus”.
• For a = 2 we find λ₁ = λ₂ = a∕2 > 0,λ₃ = -c < 0. R₁ is unstable in the e₁∕e₂ direction and stable in the e₃ direction.
• For a > 2 we find λ₁ < 0,λ₂ < 0,λ₃ < 0. R₁ is a stable node.

3.7 Tasks

1.

Calculate the solutions for the condition x_n+2 = x_n in Eqn. (3.193) for the period-2 orbit.

2.

Calculate the time derivative of I₃ in Eqn. (3.203).

3.

Calculate a few Poincaré surfaces of section at different energies for the Hénon-Heiles system.

4.

Calculate a movie out of 400-1500 Poincaré surfaces of section at different energies, where the energy is continuously changed.

5.

Solve the equations of motion for the Lorenz attractor (3.225) - (3.227) and the Rössler attractor (3.232) - (3.234) numerically.

6.

Investigate their properties (Jacobi matrix, fixed points, their stability etc.).

7.

Is this possible? Find your own equations of motion which admit an attractor solution like the Lorenz or Rössler attractor.

4 Integrating the stellar initial mass function

More information can be found in the book “An introduction to the Theory of Stellar Structure and Evolution” by D. Prialnik (Prialnik, 2000). According to Salpeter (1955), the initial mass function ξ(m) (IMF) is defined as the number of stars per unit logarithmic mass interval, or, equivalently, as “the amount of mass in stars with masses in the interval [m,m + dm] formed at a given time in a given volume of space”,

(4.240)

where m = M∕M_⊙ is the mass in solar mass units and the IMF is assumed to be a simple power law. Salpeter found a single power law index of

(4.241)

for stars in the solar neighborhood. However, the lower and upper bounds of stellar masses dictated by the laws of stellar structure and evolution imply that the Salpeter IMF is not scale free as the power law (4.240) with the exponent (4.241) suggests.

One may as well define a function ϕ(m), which is sometimes called the “stellar birth function”, as the number of stars per unit mass interval, or, equivalently, the number of stars with masses in the range [m,m + dm] formed at a given time in a given volume of space,

(4.242)

Kroupa, Tout & Gilmore (1993) found the following double-logarithmic slopes of the IMF (4.240):

(4.243)

where the most drastic change in the power law index at 0.5 M_⊙ may indicate the existence of a characteristic mass scale in the star formation process, as Kroupa et al. notice. However, the IMF given by (4.243) is a relation which has been found empirically for stars in the solar neighborhood. The slope for the higher-mass stars (above one solar mass including the whole upper main sequence), is slightly higher than the Salpeter slope (4.241).

In Kroupa (2002), the IMF is extended to the very low mass regime:

(4.244)

The 4 constants in (4.244) are obtained by solving the system of 4 linear equations,

4.1 Integral solving

To obtain the result in Eqn. (4.245) we need to evaluate the corresponding integral. This can be done analytically for the case of the Kroupa (2002) IMF in Eqn. (4.244) or numerically. For the numerical computation, there are three main approaches:

4.1.1 Solving a differential equation

We may solve the corresponding differential equation (see e. g. Press et al., 1986),

(4.247)

for N = y(120.0) with the boundary condition y(0.01) = 0. For this purpose, we may use an appropriate integrator, e. g. a Runge-Kutta method.

4.1.2 Using geometrical rules

We may use geometrical rules such as the trapezoidal rule (see e. g. Press et al., 1986):

(4.248)

with the step size 0 < h ≪ 1. Higher-order geometrical rules are the Simpson’s rule, Simpsons 3∕8 rule or Bode’s rule (see e. g. Press et al., 1986).

4.1.3 Monte-Carlo integration

We need a suitable majorizing function Φ(m) of the function ϕ(m) with ϕ(m) < Φ(m) in the whole integration interval, whose integral is known. Now we “sprinkle” the area in the integration range below the majorizing function with random points. The knowledge of the ratio of the number of points below the function ϕ(m) and the number of points below the majorizing function Φ(m) allows the determination of the integral. We note that majorizing integrability is connected to the mathematical possibility to exchange integration and limiting procedure.

4.2 Tasks

1.

Find out whether the coefficients (0.035, 0.019 and 0.019) in Equation (13) of Kroupa, Tout & Gilmore (1993) are correct.

2.

Find the mean stellar mass for the Kroupa (2002) IMF numerically and compare with the analytical result.

3.

Derive the local error for the trapezoidal rule in Eqn. (4.248).

4.

Derive the 1σ,2σ and 3σ confidence levels by numerical integration over a normal (Gaussian) distribution using the trapezoidal rule.

5 Astrophysical Poisson equations

This chapter is based on lectures and exercises by Prof. W. M. Tscharnuter at the University of Heidelberg. The theory of the Lane-Emden equation can be found in the book “Galactic dynamics” of J. Binney and S. Tremaine (Binney & Tremaine, 1987, 2008).

5.1 Lane-Emden equation

The polytropic equation of state is

(5.249)

where K is the polytropic constant, n is the polytropic index and P and ρ are pressure and density, respectively. An integration of the equation of hydrostatic equilibrium

(5.250)

where Φ is the gravitational potential, leads to

(5.251)

The Poisson equation reads then

(5.252)

Introducing dimensionless quantities s and ψ(s),

we obtain the Lane-Emden equation

(5.256)

which has analytical solutions for n = 0,n = 1 and n = 5. For other values of n, it must be solved numerically using the boundary conditions ψ(z = 0) = 1, which means that the central potential is finite and dψ∕ds|_s=0 = 0, which corresponds to a flat core of the stellar polytrope.

To integrate it, we replace it into a system of two first-order differential equations,

The solutions are shown in Figure 22.

5.1.1 Chandrasekhar mass

From (5.251) follows a relation between central potential Φ_c and central density ρ_c,

(5.259)

and we have

(5.260)

On the other hand, from (5.251) and (5.259) follows

(5.261)

Solving (5.260) for ρ_c yields

(5.262)

We now have for the cumulative mass

|In addition to the dimensionless constant M_n = -_s1, we set another dimensionless constant R_n = s₁. These constants can be found by numerically integrating the Lane-Emden equation (5.256). s₁ is just the first zero of the Lane-Emden equation, i.e. it corresponds to the stellar radius.

Plugging in (5.262) in (5.266), we obtain after some tedious calculation the important relation

(5.267)

between the polytropic constant K, the radius R and the cumulative mass M. Solving for the polytropic constant K and plugging it in the polytropic equation of state (5.249), we find

(5.268)

For the case n = 3 (relativistic degenerate electron gas), we find for the central pressure

(5.269)

where G is the gravitational constant and the numerical constant M₃ ≈ 2.01824. On the other hand, quantum statistics of a relativistic degenerate electron gas tells us that

(5.270)

where h,c,m_p are the Planck constant, the speed of light and the proton mass, respectively, and μ_e = 2. Equating these two and solving for the mass gives the famous Chandrasekhar mass

(5.271)

5.1.2 Isothermal case

The corresponding equation for the isothermal sphere is given by

(5.272)

Their corresponding system of first-oder differential equations is given by

Illustrations of the solutions are shown in Figures 23, 24, 25 and 26.

5.2 Tasks

1.

Solve the Lane-Emden equation analytically for the cases n = 0,n = 1 und n = 5.

2.

Solve the Lane-Emden equation numerically for the case n = 7∕2 and plot density, cumulative mass and potential as functions of radius.

3.

Derive Eqn. (5.267).

4.

Calculate the constant M₃.

Afterword

The author would like to thank Dr. Thomas Peters for helping out with calculations and Dr. Patrick Glaschke for teaching him the “high art” of theoretical physics. Furthermore, the author would like to thank Prof. Dr. Rainer Spurzem and Prof. Dr. Andreas Just. For the method of writing bilingual scripts the author is grateful to Prof. Dr. em. Franz Wegner at the University of Heidelberg. Last, the author would like to thank Prof. Dr. em. Werner M. Tscharnuter for interesting lectures on celestial mechanics in the winter term 2006/2007 at the University of Heidelberg from which some of the material presented here stems.

The present book would not be without the existence of the stellar dynamics working group at the Astronomisches Rechen-Institut in Heidelberg, Germany, led by Prof. Dr. Rainer Spurzem and Prof. Dr. Andreas Just, their vivid exchange of knowledge with Dr. Sverre Aarseth at the Institute of Astronomy in Cambridge, UK (you may also read the book by Aarseth, 2003), and the possibility of being involved for years in this interesting field of research, which is beyond the mainstream topics of modern physics. The author hopes, that the tradition of doing N-body calculations at the Astronomisches Rechen-Institut, which once started with Dr. Sebastian von Hoerner (you may also read the article by von Hoerner, 2001) continues to be pursued in the future.

And your life, dear reader, and mine, how shall they be pursued? I dreamed: ... La vie, la route des routes ... On my balcony grows a tiny beech tree in a blue plastic bucket. Did a bird deposit a beechnut in the bucket for hard times? Or is it going to be a blackberry bush? How is it possible to comprehend this sprouting, growing, blossoming, maturing, reproducing, prospering and lastly dying life in its entirety? I think of Hanna’s words to the Homo Faber: »You do not treat life as a living shape, but as a mere addition, therefore no relation to time, because no relation to death. Life is living shape in time.«

    References

Wikipedia, the English-German dictionary of the LEO GmbH (http://dict.leo.org) as well as the Wörterbuch der philosophischen Begriffe of the Meiner Edition (Regenbogen & Meyer, 2013) have been used as further sources of information.