Consider the motion of a star in a spherically-symmetric potential, Phi = Phi(|r|). The orbit of the star remains in a plane perpendicular to the angular momentum vector, and it's natural to adopt a polar coordinate system; call the coordinates R = |r| and phi. The system has n = 2 degrees of freedom, so the phase space has 4 dimensions.
The equations of motion can be derived by starting with the lagrangian,
1 2 2 2
(1) L(R,phi,Rdot,phidot) = - (Rdot + R phidot ) - Phi(R) ,
2
where Rdot = dR/dt and phidot = d(phi)/dt.
Differentiating with respect to Rdot and phidot
yields the momenta conjugate to R and phi,
dL
(2) ------- = Rdot = v_R ,
d(Rdot)
dL 2
(3) --------- = R phidot = R v_phi = J ;
d(phidot)
here v_R and v_phi are velocities in the radial and
azimuthal directions. The hamiltonian may now be expressed as a
function of the coordinates and conjugate momenta:
1 2 2 2
(4) H(R,phi,v_R,J) = - (v_R + J / R ) + Phi(R) .
2
Then the equations of motion are
dR dH
(5) -- = ------ = v_r ,
dt d(v_r)
d(phi) dH J
(6) ------ = -- = --- ,
dt dJ R^2
d(v_r) dH d(Phi) J^2
(7) ------ = - -- = - ------ + --- ,
dt dR dR R^3
dJ dH
(8) -- = - ------ = 0 .
dt d(phi)
Here dJ/dt = 0 because the conjugate coordinate phi
does not appear in H; phi is called a cyclic
coordinate.
The system has two integrals of motion. One, of course, is the total energy E, numerically equal to the value of H. The other is the angular momentum J. These quantities are given by
1 2 2
(9) E = - ( v_R + v_phi ) + Phi(R) ,
2
(10) J = R v_phi .
Each of these integrals defines a hypersurface in phase space, and the
orbit must remain in the intersection of these hypersurfaces. This
can be visualized by ignoring the phi coordinate and drawing
surfaces of constant E and J in the
three-dimensional space (R,v_R,v_phi). Surfaces of constant
E are figures of revolution about the R axis, while
surfaces of constant J are hyperbolas in the
(r,v_phi) plane. The intersection of these surfaces is a
closed curve, and an orbiting star travels around this curve.
For an orbit of a given J, the system may be reduced to one degree of freedom by defining the effective potential,
J^2
(11) Psi(R) = Phi(R) + ----- ;
2 R^2
the corresponding equations of motion are then just
dR
(12) -- = v_R ,
dt
d(v_R) d(Psi)
(13) ------ = - ------ .
dt dR
Because Psi(R) diverges as R -> 0, the star is
energetically prohibited from coming too close to the origin, and
shuttles back and forth between turning points R_min and
R_max.
In addition to its periodic radial motion described by Eqs. 4 & 5, a star also executes a periodic azimuthal motion as it orbits the center of the potential. If the radial and azimuthal periods are incommensurate, as is usually the case, the resulting orbit never returns to its starting point in phase space; in coordinate space such an orbit is a rosette. The Keplerian potential is a very special case in which the radial and azimuthal periods of all bound orbits are equal. The only other potential in which all orbits are closed is the harmonic potential generated by a uniform sphere; here the radial period is half the azimuthal one and all bound orbits are ellipses centered on the bottom of the potential well. Thus in the Keplerian case all stars advance in azimuth by Delta phi = 2 pi between successive pericenters, while in the harmonic case they advance by Delta phi = pi. Galaxies typically have mass distributions intermediate between these extreme cases, so most orbits in spherical galaxies are rosettes advancing by pi < Delta phi < 2 pi between pericenters.
Any system in which stellar collisions are rare may be idealized as a collection of N point-sized bodies, each with mass m_i, position r_i, and velocity v_i. The hamiltonian for such a system is
-- --
1 \ 2 G \ m_i m_j
(14) H(r,v) = - | m_i v_i - - | ----------- ,
~ ~ 2 / 2 / |r_i - r_j|
-- --
i i, j#i
where H depends on all body positions and velocities, the
first sum runs over all N bodies, the second runs over all
pairs of bodies (twice, hence the factor of 1/2), and
G is the gravitational constant. Then the equations of
motion are
--
d(r_i) d(v_i) \ m_j (r_i - r_j)
(15) ------ = v_i , ------ = -G | --------------- ,
dt dt / |r_i - r_j|^3
--
j#i
where the sum runs over all bodies except body i.
The famous Aarseth code in Fortran (Aarseth.f)
or C (Aarseth.c) calculates the general motion
of the N-body system in computer simulations.
N-body systems obey several basic conservation laws. A symmetry of the hamiltonian is a transformation which leaves the physical system unchanged. For example, translation in time, t -> t + Delta t, is a symmetry of Eq. 14 because H is not an explicit function of time; consequently the total system energy E = K + U = H is conserved. Likewise, symmetry with respect to translation in space, r -> r + Delta r, implies conservation of total linear momentum, and symmetry with respect to rotation gives rise to conservation of total angular momentum.
Another general result shown by manipulating Eq. 15 is the scalar virial theorem which states that for a system in equilibrium,
(16) 2 < K > + < U > = 0 ,where K and U are the total kinetic and potential energy, respectively, and the angle-brackets indicate time-averages. Since E = K+U, the time-averaged kinetic and potential energies are related to the conserved total energy by
(17) < K > = - E , < U > = 2 E .
The total mass M and total energy E of an N-body system thus define characteristic velocity and length scales
2 2
2 < K > |E| M M
(18) V = 2 --- = 2 --- , R = - G ------= G ---- .
M M < U > 2|E|
These are sometimes known as the virial velocity and radius,
respectively.The quantity t_c = R/V is an estimate of the time a typical star takes to cross the system. This timescale may be expressed in several different ways; for example, in terms of the total mass M and energy E, it is
1/2
(19) t_c = G (M^5 / 8 |E|^3) .
Note that M and E are conserved, so t_c is
a constant even for systems which are far from dynamical equilibrium.
In such cases t_c approximates the time-scale over which the
system evolves toward equilibrium.
Another expression for t_c follows from the substitution V^2 = G M/R valid for systems near equilibrium:
-1/2
(20) t_c = (G M / R^3) .
Here the quantity M/R^3, which has units of density, appears.
In systems with galaxy-like density profiles, the virial radius is
approximately proportional to the half-mass radius: R = 2.5
R_h. Using this relationship, it follows that
-1/2
(21) t_c = ~1.36 (G rho_h) ,
where rho_h is the mean density within R_h. Since
the crossing time is just supposed to indicate a typical time-scale
for orbital motion, it is usual to drop the numerical constant, and
define
-1/2
(22) t_c = (G rho_h) .