Dynamical equilibrium in galaxies

Galaxies are dynamical systems. Stars in galaxies orbit in the gravitational potential contributed by other stars, even those which are further away. The potential and hence the orbits of stars are determined by the gross structure of the galaxy, rather than the action of some nearby stars. The force acting on stars thus does not fluctuate rapidly. We have also seen previously that collisions and close encounter are also quite rare in galaxies (although they can be important for systems such as globular clusters).

Galaxies like the Milky way have trillions of stars, each orbiting in a common potential. The behaviour of such a system can be described by the phase space distribution function \(f(\vec{x}, \vec{v}, t)\). Let us call the phase space coordinates \(\vec{w}\). The function \(f(\vec{w}, t)d^6\vec{w}\) gives the number of particles, \(dN\), within a given phase space volume. We have

\begin{equation} N = \int d^6\vec{w} f(\vec{w}, t)\,, \end{equation}

and that we can obtain the density by integrating the phase space function over velocities,

\begin{equation} n(\vec{x}, t) = \int d^3\vec{v} f(\vec{w}, t)\,. \end{equation}

The rate of change of particles occupying a small phase space volume \(dV_{\rm ps} = d^6\vec{w}\), will then be given by

\begin{equation} \frac{dN}{dt} = \int \frac{\partial}{\partial t}f(\vec{w}, t)d^6\vec{w}\,. \end{equation}

These particles must be leaving the phase space volume through its boundaries rather than jumping off into different locations in phase space (there are no collisions which will cause direct jumps between unconnected with each other), and hence must equal the divergence of the phase space flow \(f \dot{\vec{w}}\) through the same volume,

\begin{eqnarray} \frac{dN}{dt} &=& -\int \nabla_{\vec{w}} . \left[ f(\vec{w}, t)\,\vec{\dot{w}}\right] d^6 \vec{w} \end{eqnarray}

Equating the two we get

\begin{equation} \int d^6\vec{w} \left( \frac{\partial}{\partial t}f(\vec{w}, t) + \nabla_{\vec{w}} . \left[ f(\vec{w}, t)\,\dot{\vec{w}}\right] \right) = 0 \end{equation}

Since this holds for any phase space volume, the integrand itself must equal zero. Thus we can obtain the collisionless Boltzmann equation,

\begin{equation} \frac{\partial}{\partial t}f(\vec{w}, t) + \nabla_{\vec{w}} . \left[ f(\vec{w}, t)\,\dot{\vec{w}}\right] = 0 \end{equation}

The divergence can be separated out into the \(\{\vec{x},\,\vec{v}\}\) components. In our case, \(\partial{\vec{v}}/\partial{\vec{x}}=0\), and \(\partial{\dot{\vec{v}}}/\partial{\vec{v}}=0\), we obtain

\begin{eqnarray} 0&=&\frac{\partial}{\partial t}f(\vec{w}, t) + \dot{\vec{x}} . \nabla_{\vec{x}} f(\vec{w}, t) + \dot{\vec{v}} . \nabla_{\vec{v}} f(\vec{w}, t) \\ &=&\frac{\partial}{\partial t}f(\vec{w}, t) + \vec{v} . \nabla_{\vec{x}} f(\vec{w}, t) - \frac{\partial\Phi}{\partial\vec{x}} . \nabla_{\vec{v}} f(\vec{w}, t) \end{eqnarray}

It is very difficult to find a direct general solution to the Boltzmann equation in this form. But progress can be made by considering different moments of this equation. We define

\begin{eqnarray} n(\vec{x}) &=& \int f d^3\vec{v}\\ \langle{v_i}\rangle(\vec{x}) &=& \frac{1}{n(\vec{x})}\int f v_i d^3\vec{v} \\ \sigma^2_{ij}(\vec{x}) &=& \langle v_iv_j\rangle(\vec{x}) - \langle{v_i}\rangle\langle{v_j}\rangle = \frac{1}{n(\vec{x})}\int f \left[v_i v_j-\langle{v_i}\rangle\langle{v_j}\rangle\right] d^3\vec{v} \end{eqnarray}

The zeroth moment of the Boltzmann equation is given by

\begin{eqnarray} 0 &=& \int d^3\vec{v} \frac{\partial}{\partial t}f(\vec{w}, t) + \int d^3\vec{v} \,\vec{v} . \nabla_{\vec{x}} f(\vec{w}, t) - \int d^3\vec{v} \frac{\partial\Phi}{\partial\vec{x}} . \nabla_{\vec{v}} f(\vec{w}, t) \\ &=& \frac{\partial}{\partial t}\int d^3\vec{v} f(\vec{w}, t) + \nabla_{\vec{x}} . \left[\int d^3\vec{v} \,\vec{v} f(\vec{w}, t)\right] - \int d^3\vec{v} \nabla_{\vec{v}}. \left(\frac{\partial\Phi}{\partial\vec{x}} f(\vec{w}, t)\right) \\ &=& \frac{\partial n}{\partial t} + \nabla_{\vec{x}} . \left[n \langle \vec{v}\rangle\right] - \int d^3\vec{v} \nabla_{\vec{v}}. \left(\frac{\partial\Phi}{\partial\vec{x}} f(\vec{w}, t)\right) \frac{\partial n}{\partial t} + \nabla_{\vec{x}} . \left[n \langle \vec{v}\rangle\right] \\ &=& \frac{\partial n}{\partial t} + \nabla_{\vec{x}} . \left[n \langle \vec{v}\rangle\right] \end{eqnarray}

The last term can be shown to be zero by using the divergence theorem and noting that all realistic distribution functions which should die down to zero as \(|v|\rightarrow\infty\). The zeroth moment thus recovers the continuity equation.

The next order moment of the Boltzmann equation is

\begin{eqnarray} 0 &=& \int \vec{v} \frac{\partial}{\partial t}f(\vec{w}, t) d^3\vec{v} + \int \vec{v} \,\vec{v} . \nabla_{\vec{x}} f(\vec{w}, t) d^3\vec{v} - \int \vec{v} \frac{\partial\Phi}{\partial\vec{x}} . \nabla_{\vec{v}} f(\vec{w}, t) d^3\vec{v} \\ &=& \frac{\partial n \langle \vec{v}\rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2+n\langle \vec{v}\rangle \langle \vec{v}\rangle\right] - \int \vec{v} \nabla_{\vec{v}}. \left(\frac{\partial\Phi}{\partial\vec{x}} f(\vec{w}, t)\right) d^3\vec{v}\\ &=& \frac{\partial n \langle \vec{v}\rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2+n\langle \vec{v}\rangle \langle \vec{v}\rangle\right] - \int \nabla_{\vec{v}}. \left(\vec{v}\frac{\partial\Phi}{\partial\vec{x}} f(\vec{w}, t)\right) d^3\vec{v} + \int \frac{\partial\Phi}{\partial\vec{x}} f(\vec{w}, t) \nabla_{\vec{v}}. \vec{v} d^3\vec{v}\\ &=& \frac{\partial n \langle \vec{v}\rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2+n\langle \vec{v}\rangle \langle \vec{v}\rangle\right] - \int \nabla_{\vec{v}}. \left(\vec{v}\frac{\partial\Phi}{\partial\vec{x}} f(\vec{w}, t)\right) d^3\vec{v} + n \frac{\partial\Phi}{\partial\vec{x}}\\ &=& \frac{\partial n \langle \vec{v}\rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2+n\langle \vec{v}\rangle \langle \vec{v}\rangle\right] + n \frac{\partial\Phi}{\partial\vec{x}}\\ &=& \langle\vec{v}\rangle \frac{\partial n}{\partial t} + n\frac{\partial \langle \vec{v} \rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2\right] + \langle \vec{v}\rangle \nabla_{\vec{x}} . \left[n\langle \vec{v}\rangle\right] + n \langle \vec{v}\rangle\nabla_{\vec{x}} . \left[\langle \vec{v}\rangle\right] + n \frac{\partial\Phi}{\partial\vec{x}}\\ &=& \langle\vec{v}\rangle \left(\frac{\partial n}{\partial t} + \nabla_{\vec{x}} . \left[n \langle \vec{v}\rangle\right]\right) + n\frac{\partial \langle \vec{v} \rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2\right] + n \langle \vec{v}\rangle\nabla_{\vec{x}} . \left[\langle \vec{v}\rangle\right] + n \frac{\partial\Phi}{\partial\vec{x}}\\ &=& n\frac{\partial \langle \vec{v} \rangle}{\partial t} + \nabla_{\vec{x}} . \left[n \sigma^2\right] + n \left[\nabla_{\vec{x}} . \langle \vec{v}\rangle\right] \langle \vec{v}\rangle + n \frac{\partial\Phi}{\partial\vec{x}} \end{eqnarray}

The component forms of the last vector equation are called the Jeans equations. These correspond to the Euler equations in fluid dynamics with the tensor \(n\sigma^2\) acting as the pressure (per unit mass). Unlike fluids which have a uniform pressure in all directions, in our case at each location there is a dispersion ellipsoid with principal axes.

It is difficult to also solve the Jeans equations as such. Given a density distribution \(n(x)\), they have nine unknowns, 3 streaming velocities and 6 components of the symmetric \(\sigma^2\) tensor. Thus we have an incomplete set of equations. Adding higher order moments does not help, and we have to resort to symmetries to make further progress.

Spherically symmetric case

For the spherically symmetric case, the Boltzmann equation will be

\begin{equation} \frac{df}{dt} = \frac{\partial f}{\partial t} + \dot{r} \frac{\partial f}{\partial r} + \dot{\theta} \frac{\partial f}{\partial \theta} + \dot{\phi} \frac{\partial f}{\partial \phi} \end{equation}

One can follow the steps above in the spherically symmetric coordinate system to derive the three Jeans equations for spherically symmetric case. If we assume that the system is static, then the partial derivatives with respect to time vanish. Spherical symmetry will also mean that there are no streaming motions in any directions, and the \(\sigma^2\) tensor is diagonal with \(\sigma_{\theta}^2=\sigma_{\phi}^2\). This finally leads us to

\begin{equation} \frac{\partial \left(n\sigma^2_r\right)}{\partial r} + \frac{2n}{r}\left[\sigma^2_r-\sigma^2_\theta \right] + n\frac{d\Phi}{dr} = 0 \end{equation}

This equation still has two unknowns and is a single equation. If we define

\begin{equation} \beta(r) = 1-\frac{\sigma^2_\theta}{\sigma^2_r}\, \end{equation}

then we obtain

\begin{equation} \frac{\partial \left(n\sigma^2_r\right)}{\partial r} + \frac{2n\sigma_r^2\beta}{r} + n\frac{d\Phi}{dr} = 0\,. \end{equation}

This simplifies to

\begin{eqnarray} \frac{r}{n\sigma^2_r}\frac{\partial \left(n\sigma^2_r\right)}{\partial r} + 2\beta &=& - \frac{r}{\sigma^2_r}\frac{d\Phi}{dr} \\ &=& \frac{GM(\lt r)}{\sigma_r^2r}\,, \end{eqnarray}

which implies that the mass profile can be obtained as

\begin{eqnarray} M(\lt r) = \frac{r\sigma^2_r}{G} \left( \frac{\partial \ln n}{\partial \ln r} + \frac{\partial \ln\sigma^2_r}{\partial \ln r} + 2\beta \right) \end{eqnarray}

If we assume a constant value for the velocity anisotropy \(\beta=0\), we can use the integrating factor method to solve the above equation. We will obtain

\begin{equation} \sigma_r^2(r) = \frac{1}{n(r)r^{2\beta}} \int n(r') r'^{2\beta} \frac{GM(\lt r')}{r'^2} dr' \end{equation}

The observables that we can have access to are the projected surface density of stars,

\begin{equation} \Sigma(R) = \int_{0}^{\infty} n(r') \frac{2r'dr'}{\sqrt{r^2-R^2}}\, \end{equation}

and the line of sight velocity dispersion,

\begin{equation} \sigma^2_{\rm los}(R) = \frac{1}{\Sigma(R)} \int_{0}^{\infty} \left(1-\beta\frac{R^2}{r'^2} \right) n(r') \sigma^2_r(r')\frac{2r'dr'}{\sqrt{r^2-R^2}}\, \end{equation}

Depending upon the system one might have only access to aperture averaged velocity dispersions only. In that case one has to further integrate the above equation within an aperture to obtain the aperture averaged velocity dispersion,

\begin{equation} \sigma^2_{\rm los}(\lt R) = \frac{\int_{0}^{R} dR' R' \sigma^2_{\rm los}(R) \Sigma(R')}{\int_0^R dR' R' \Sigma(R')} \end{equation}

These equations then allow to relate parametric mass profiles to the observables \(\Sigma\) and \(\sigma_{\rm los}(\lt R)\).

Axially symmetric case

Asymmetric drift

Tensor virial theorem

We integrated the collisionless Boltzmann equation by taking the first velocity moment and integrating over velocity space. This gave (now in index notation)

\begin{equation} \frac{\partial n\langle v_j \rangle}{\partial t} + \frac{\partial n\langle v_i\rangle\langle v_j\rangle}{\partial x_i} + \frac{\partial n\sigma^2_{ij}}{\partial x_i} + n\frac{\partial \Phi}{\partial x_j} = 0\,. \end{equation}

One can then take the first position moment of this equation and integrate over all positions,

\begin{equation} \int d\vec{x} \,\,x_k\frac{\partial n\langle v_j \rangle}{\partial t} + \int d\vec{x} \,\,x_k\frac{\partial n\langle v_i\rangle\langle v_j\rangle}{\partial x_i} + \int d\vec{x} \,\,x_k\frac{\partial n\sigma^2_{ij}}{\partial x_i} + \int d\vec{x} \,\,x_k n\frac{\partial \Phi}{\partial x_j} = 0\,. \end{equation}

We can integrate the second and the third terms,

\begin{eqnarray} {\rm Term\,II\,}:\int d\vec{x} \,\,x_k\frac{\partial n\langle v_i\rangle\langle v_j\rangle}{\partial x_i} &=& \int d\vec{x} \,\,\frac{\partial n x_k\langle v_i\rangle\langle v_j\rangle}{\partial x_i} - \int d\vec{x} \,n \langle v_i\rangle\langle v_j\rangle \frac{\partial x_k}{\partial x_i} \\ &=& \int d\vec{x} n \langle v_k\rangle\langle v_j\rangle\\ {\rm Term\,III\,}:\int d\vec{x} \,\,x_k\frac{\partial n\sigma^2_{ij}}{\partial x_i} &=& \int d\vec{x} \,\,\frac{\partial n x_k\sigma^2_{ij}}{\partial x_i} - \int d\vec{x} \,n \sigma^2_{ij} \frac{\partial x_k}{\partial x_i} \\ &=& \int d\vec{x} n \sigma^2_{kj}\\ \end{eqnarray}

These two terms help define the kinetic energy tensor

\begin{equation} {\cal K}_{jk} = \frac{1}{2} \int d\vec{x}\,\, n\left[\langle v_k\rangle\langle v_j\rangle + \sigma^2_{kj}\right] \end{equation}

The last term can be used to define the potential energy tensor

\begin{equation} {\cal W}_{jk} = \int d\vec{x}\,\,x_k n\frac{\partial \Phi(x)}{\partial x_j} \end{equation}

We can substitute the equation for the potential as an explicit integral over another variable \(x'\).

\begin{eqnarray} {\cal W}_{jk} = -\int d\vec{x}\,\,x_k n\frac{\partial}{\partial x_j}\left[ \int d{\vec x'} G\frac{\rho(\vec{x'})}{|\vec{x'}-\vec{x}|} \right] \end{eqnarray}

\begin{equation} {\cal W}_{jk} = -m\int\int d\vec{x} d\vec{x'} \,\,x_k n(\vec{x}) \left( G\frac{n(\vec{x'})(x'_j-x_j)}{|\vec{x'}-\vec{x}|^3} \right) \end{equation}

But note that the same equation can be also written as (by starting the definition with an integral over \(\vec{x}'\) and then expressing potential as an integral over \(\vec{x}\))

\begin{eqnarray} {\cal W}_{jk} &=& -m\int\int d\vec{x'} d{\vec x} \,\,x'_k n(\vec{x'}) \left( G\frac{n(\vec{x})(x_j-x'_j)}{|\vec{x}-\vec{x'}|^3} \right) \\ \end{eqnarray}

Thus \(W_{jk}\) can be expressed as half the sum of the above two expressions

\begin{eqnarray} {\cal W}_{jk} &=& -\frac{m}{2}\int\int d\vec{x'} d{\vec x} \,\,(x_k-x'_k) (x_j-x'_j) n(\vec{x})n(\vec{x'}) \left( G\frac{1}{|\vec{x}-\vec{x'}|^3} \right) \\ \end{eqnarray}

Finally one can show that the first term in our starting equations can be shown to be equal to the second derivative of the moment of inertia tensor.

\begin{equation} {\cal I}_{ij} = \int n x_i x_j d{\vec{x}} \end{equation}

We get the tensor form of the virial theorem,

\begin{equation} \frac{d^2I_{jk}}{dt^2} = 2{\cal K}_{jk} + {\cal W}_{jk} \end{equation}

The trace of this equation gives the well known scalar virial theorem.

Applications of the virial theorem

Masses of galaxy clusters

The flattening of rotating oblate elliptical galaxies

Jeans theorem

We had considered integrals of motion along an orbit in our lectures about orbit theory. These integrals of motion are just functions of \(\vec{x}\) and \(\vec{v}\). Thus,

\begin{equation} I(\vec{x}, \vec{v}) = {\rm constant} \end{equation}

The total derivative of the integral of motion (derivative along the orbit) is zero. Hence

\begin{equation} \frac{d}{dt}I(\vec{x}, \vec{v}) = 0 = \vec{v}.\nabla I - \nabla \Phi.\frac{\partial I}{\partial \vec{v}} \end{equation}

Comparing with the collisionless Boltzmann equation, we see that \(I\) is a solution to this equation. Given this, one can see that the following theorem due to Jeans should hold true.

Jeans theorem: Any function of the integrals of motion \(f(\{I_n\})\) is a steady state solution to the collisionless Boltzmann equation, and any steady state solution of the equation depends on \([\vec{v}, \vec{x}]\) only through integral of motion.

Now remember that we had two types of integrals of motion, the isolating and the non-isolating type. The non-isolating type do not reduce the dimensions of phase space available to the particle, while the isolating ones do. These non-isolating integrals of motion are not important. Recall regular orbits are those which have N=3 isolating integrals of motion. There is a stronger version of Jeans theorem which corresponds to systems which have regular orbits.

Strong Jeans theorem: The distribution function in a system which has regular orbits can be presumed to be a function only of the three independent isolating integrals of motion.

Proving this theorem requires use of action angle variables from classical mechanics. Refer to BT Appendix 4.A for the proof. But the consequence of this theorem is that we can write down the distribution function as a function of the isolating integrals of motion.

Spherically symmetric case

Let us consider spherical symmetry for example. We know of 4 integrals of motion \((E, \vec{L})\). Assuming that all properties of the system are spherically symmetric (this is not necessarily a fully justified assumption, see e.g., Maxwell’s demon from Lynden-Bell 1960), then the distribution function \(f\) will only depend upon the magnitude of \(L\) and \(E\). Let us first consider the simple case where \(f\) depends only upon \(E\).

The velocity dispersion in the spherically symmetric case will be obtained as

\begin{equation} \langle v_{\rm r}^2\rangle = \int dv_r dv_\theta dv_\phi v_r^2 f\left(\frac{1}{2}\left[v_r^2 + v_\theta^2 + v_\phi^2 \right] + \Phi \right) \end{equation}

However, the distribution function depends similarly upon \((v_r, v_\theta, v_\phi)\). Therefore, in such a system we should have \(\sigma_r^2=\sigma_\theta^2=\sigma_\phi^2\), i.e., we have an isotropic velocity dispersion. The average streaming motion is equal to zero.

We can consider a relative potential \(\Psi = -\Phi+\Phi_0\), where \(\Phi_0\) is some constant. The relative energy is defined to be \({\cal E} = -E+\Phi_0 = \Psi-1/2v^2\). For isolated systems one can choose \(\Phi_0=0\), in general it is chosen such that the distribution function \(f({\cal E})\gt 0\) for \({\cal E}\gt 0\), and is equal to zero otherwise.

Let us say we want to find density-distribution function pairs (just like we found density-potential pairs). The density distribution is the zeroth moment of the distribution function. Thus,

\begin{equation} n(\vec{x}) = \int fd^3{\vec v} = 4\pi \int_0^{\Psi} f({\cal E})\sqrt{2(\Psi-{\cal E})} d{\cal E} \end{equation}

One can differentiate this equation with respect to \(\Psi\) to obtain

\begin{equation} \frac{d n}{d {\Psi}}=2\sqrt{2}\pi \int_0^{\Psi} f({\cal E})\frac{1}{\sqrt{\Psi-{\cal E}}} d{\cal E} \end{equation}

This is an Abel integral which can then be inverted to obtain \(f({\cal E})\) once we specify a density.

\begin{equation} f({\cal E}) = \frac{1}{2\sqrt{2}\pi^2} \frac{d}{d{\cal E}} \left(\int_0^{\cal E} d\Psi \frac{1}{\sqrt{{\cal E}-\Psi}} \frac{dn}{d\Psi} \right) \end{equation}

As we need \(f({\cal E})\gt 0\) for \({\cal E}\gt 0\), the quantity inside the brackets needs to be an increasing function of \({\cal E}\). Thus one could always find a distribution function \(f({\cal E})\) for a given \(\rho\), but it is not necessarily guaranteed to be physical. In this case, one can try to posit some form of \(f(E,L)\) instead to find a solution for the distribution function.

Another option to obtain density-distribution function pairs, is by starting from a parametric form for the distribution function. We will consider some special cases, say the distribution function \(f\) is such that

\begin{equation} f({\cal E}) = \frac{\rho_1}{{(2\pi\sigma^2)}^{3/2}}\exp\left[ \frac{{\cal E}}{\sigma^2} \right] = \frac{\rho_1}{{(2\pi\sigma^2)}^{3/2}}\exp\left[ \frac{\Psi-v^2/2}{\sigma^2} \right] \end{equation}

The density distribution is given by integrating over velocities

\begin{equation} \rho = \rho_1 \exp\left[ \frac{\Psi}{\sigma^2} \right] \end{equation}

From the definition of \(\Psi\), it can be shown that it obeys the Poisson equation with a negative sign

\begin{equation} \nabla^2 \Psi = -4\pi G\rho \end{equation}

This will help us to obtain the dependence of \(\rho(r)\)

\begin{equation} \frac{1}{r^2}\frac{d}{d r} \left( r^2 \sigma^2 \frac{d \ln \rho}{d r}\right) = -4\pi G\rho \end{equation}

and the equation corresponds to the density profile

\begin{equation} \rho(r) = \frac{\sigma^2}{2\pi G r^2} \end{equation}

Poisson equation

The Poisson equation for a system which itself generates the potential will be given by

\begin{equation} \nabla^2 \Phi = 4\pi G\rho = 4\pi G \int f(E) d^3{\vec v} \end{equation}

Given the spherical symmetry, we can write \(d^3v = 4\pi v^2dv = 4\pi \sqrt{2(E-\Phi)} dE\), because \(dE = vdv\), and hence

\begin{equation} \frac{1}{r^2} \frac{\partial }{\partial r} \left( r^2 \frac{\partial \Phi}{\partial r}\right) = 16\pi^2 G \int f(E) \sqrt{2(E-\Phi)} dE \end{equation}

Appendix

Notes on partial derivatives

While deriving the above equations, one has to remember that in phase space, \(\vec{x}\) and \(\vec{v}\) are independent coordinates. When we take partial derivatives of \(\vec{x}\) with respect to any components of \(\vec{v}\), it should be zero as we will keep the rest of the phase space coordinates fixed, i.e.,

\begin{equation} \left.\frac{\partial x_i}{\partial v_j}\right|_{\rm fixed\,\,x_i}= 0 \end{equation}

For the partial time derivatives as well, if we are sitting at a fixed phase space point, the partial derivatives of the phase space coordinates \({\vec{x}, \vec{v}}\) with respect to time will be equal to zero.

Note however, that if we integrate over velocities and obtain some velocity averaged quantity as a function of position, then we can take its partial time derivative with respect to time. For example, the average streaming velocities \(\langle\vec{v}\rangle\) could depend upon position, and so can have partial derivatives with respect to both time and space.