← Previous section Next section →

14.1. The tensor virial theorem

\label{sec-ellipequil-tensorvirial}

We introduced the virial theorem in Chapter 5.2 and specified there that the theorem that we derived was more specifically known as the scalar virial theorem. The scalar virial theorem is a simple, yet powerful tool for understanding the overall structure and mass distribution of a gravitating system. For example, using it to estimate the mass of a gravitating system allowed Fritz Zwicky to uncover the first evidence for dark matter from the large velocity dispersion of galaxies in the Coma cluster (see Chapter 6.1). Or later in Chapter 17.3.1, we will use the scalar virial theorem to determine the size of a just-formed, virialized dark-matter halo. But while the scalar virial theorem can give us much information about the overall mass of a system, it is a blunt tool to understand more complex systems such as elliptical galaxies and other flattened systems.

14.1.1. Collisional systems

We can derive a more sophisticated, tensorial form of the virial theorem by starting from the virial tensor \(\vec{G}_{\alpha\beta}\) rather than the scalar virial that we used in Chapter 5.2, defined as \begin{equation} \vec{G}_{\alpha\beta} = \frac{1}{2}\sum_{i=1}^N w_i\, \left[x^\alpha_i\,v^\beta_i+x_i^\beta\,v_i^\alpha\right]\,, \end{equation} where \(x^\alpha_i\) is the \(\alpha\)-th component of the position vector \(\vec{x}\) of particle \(i\), \(v^\beta_i\) is the \(\beta\)-th component of the velocity vector \(\vec{v}_i\) of particle \(i\), and the \(w_i\) are an arbitrary set of weights. Note that we have explicitly defined the tensor so as to be symmetric. By comparing to the definition of the scalar virial \(G\) in Equation (5.11), we see that \(G\) is simply the trace of \(\vec{G}_{\alpha\beta}\). Again assuming that the system of \(N\) particles is in equilibrium, we can set \(\mathrm{d} \vec{G}_{\alpha\beta} / \mathrm{d}t = 0\), which gives \begin{align} \frac{\mathrm{d} \vec{G}_{\alpha\beta}}{\mathrm{d} t} & = \frac{1}{2}\sum_{i=1}^N w_i\,\left[ \dot{x}^\alpha_i\,v^\beta_i + x^\alpha_i\,\dot{v}^\beta_i +\dot{x}_i^\beta\,v_i^\alpha + x_i^\beta\,\dot{v}_i^\alpha \right]\\ & = \sum_{i=1}^N w_i\,\left[ v^\alpha_i\,v^\beta_i + \frac{1}{2}x^\alpha_i\,g^\beta(\vec{x}_i) + \frac{1}{2}x_i^\beta\,g^\alpha(\vec{x}_i) \right]\,,\label{eq-tensorvirial-first} \end{align} where \(g^\beta(\vec{x})\) is the \(\beta\)-th component of the gravitational field at \(\vec{x}\).

The gravitational field \(\vec{g}(\vec{x})\) can either be an external field, in which case it can simply be plugged into the expression above, or it can be the field resulting from the mutual gravity between all of the \(N\) objects in the system. In the latter case, we can re-write the last two terms of Equation (14.3) with a similar set of manipulations as in Equations (5.16); we define this tensor to be the potential-energy tensor \(\vec{W}_{\alpha\beta}\) of a system of point particles \begin{align} \vec{W}_{\alpha\beta} & = -\frac{1}{2}\,\sum_{i=1}^N \sum_{j< i}\left[\left(w_i\,Gm_j\,x_i^\alpha-w_j\,Gm_i\,x_j^\alpha\right)\,\frac{x_i^\beta-x_j^\beta}{|\vec{x}_i-\vec{x}_j|^3} + (\alpha \leftrightarrow \beta)\right]\,, \end{align} where \((\alpha \leftrightarrow \beta)\) indicates that we repeat the first term with the \(\alpha\) and \(\beta\) indices swapped. When we set \(w_i \equiv m_i\) for all particles and use the symmetry under \(i\leftrightarrow j\), this simplifies to \begin{equation}\label{eq-potenergytensor-collisional} \vec{W}_{\alpha\beta} = -\frac{1}{2}\,\sum_{i=1}^N \sum_{j\neq i}Gm_i\,m_j\,\left[\frac{(x_i^\alpha-x_j^\alpha)\,(x_i^\beta-x_j^\beta)}{|\vec{x}_i-\vec{x}_j|^3}\right]\,. \end{equation} Similarly, we define the kinetic energy tensor \(\vec{K}_{\alpha\beta}\) as \begin{equation} \vec{K}_{\alpha\beta} = \frac{1}{2}\sum_{i=1}^N m_i\,v^\alpha_i\,v^\beta_i\,, \end{equation} and then Equation (14.3) can be written as \begin{align}\label{eq-tensorvirial-second} \frac{\mathrm{d} \vec{G}_{\alpha\beta}}{\mathrm{d} t} & = 2\,\vec{K}_{\alpha\beta} +\vec{W}_{\alpha\beta}\,. \end{align} This is the tensor virial theorem. In a steady state system, \(\mathrm{d} \vec{G}_{\alpha\beta} / \mathrm{d}t = 0\), such that this becomes \begin{align}\label{eq-tensorvirial-steady} 2\,\vec{K}_{\alpha\beta} +\vec{W}_{\alpha\beta} & = 0\,. \end{align} Taking the trace of this equation, we obtain the scalar virial theorem of Equation (5.20) again.

14.1.2. Collisionless systems

\label{sec-tensorvirial-collisionless}

For collisionless systems, such as dark-matter halos or smooth models of galaxies, we can derive a continuous version of the tensor (and scalar) virial theorem. To do this, we define the virial tensor with weights \(w(\vec{x}) = \rho(\vec{x},t)\,\mathrm{d} \vec{x}\), using the mean velocity \(\bar{\vec{v}}(\vec{x},t)\) for the velocity, and we integrate over all of position space: \begin{equation} \vec{G}_{\alpha\beta} = \frac{1}{2}\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\left[x^\alpha\,\bar{v}^\beta+x^\beta\,\bar{v}^\alpha\right]\frac{1}{2}\int\mathrm{d}\vec{x}\,\mathrm{d}\vec{v}\,f(\vec{x},\vec{v},t)\,\left[x^\alpha\,v^\beta+x^\beta\,v^\alpha\right]\,, \end{equation} where in the second equality we have used that \(\bar{\vec{v}}(\vec{x},t) = \int\mathrm{d}\vec{v}\,f(\vec{x},\vec{v},t)/\rho(\vec{x},t)\), where \(f(\vec{x},\vec{v},t)\) is the distribution function. This is a natural definition of the virial tensor for a collisionless system (for which the phase-space distribution function is the fundamental quantity). Again taking the time derivative, we obtain now \begin{equation} \frac{\mathrm{d}\vec{G}_{\alpha\beta}}{\mathrm{d}t} = \frac{1}{2}\int\mathrm{d}\vec{x}\,\mathrm{d}\vec{v}\,f(\vec{x},\vec{v},t)\,\left[2\,v^\alpha\,v^\beta+x^\alpha\,\dot{v}^\beta+x^\beta\,\dot{v}^\alpha\right]\,,\label{eq-tensorvirial-cont-first} \end{equation} where we have used the fact that the total time derivative of the distribution function vanishes (Liouville’s theorem; Equation 5.34). We again split the right-hand side of this equation into a kinetic- and potential-energy tensor. First, we define the kinetic energy tensor for continuous, collisionless systems as \begin{equation} \vec{K}_{\alpha\beta} = \frac{1}{2}\,\int\mathrm{d}\vec{x}\,\mathrm{d}\vec{v}\,f(\vec{x},\vec{v},t)\,v^\alpha\,v^\beta\,. \end{equation} We can also write this as \begin{align} \vec{K}_{\alpha\beta} & = \frac{1}{2}\,\int\mathrm{d}\vec{x}\,\,\rho(\vec{x},t)\,\overline{v^\alpha\,v^\beta} = \vec{T}_{\alpha\beta} +\frac{1}{2}\,\vec{\Pi}_{\alpha\beta}\,, \end{align} where \(\vec{T}_{\alpha\beta}\) and \(\vec{\Pi}_{\alpha\beta}\) are the contributions from ordered and random motion to the kinetic energy \begin{align}\label{eq-ellipequil-disptensor-ordered} \vec{T}_{\alpha\beta} & = \frac{1}{2}\,\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\bar{v}^\alpha\,\bar{v}^\beta\,,\\ \vec{\Pi}_{\alpha\beta} & = \phantom{\frac{1}{2}\,}\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\overline{(v-\bar{v})^\alpha\,(v-\bar{v})^\beta} = \,\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\vec{\sigma}_{\alpha\beta}(\vec{x})\,,\label{eq-ellipequil-disptensor-random} \end{align} where \(\vec{\sigma}_{\alpha\beta}\) is the velocity dispersion tensor.

The last two terms of Equation (14.10) are the potential energy tensor for a continuous system and can be worked out as \begin{align}\label{eq-pot-energy-tensor} \vec{W}_{\alpha\beta} & = \frac{1}{2}\int\mathrm{d}\vec{x}\,\mathrm{d}\vec{v}\,f(\vec{x},\vec{v},t)\,\left[x^\alpha\,\frac{\partial \Phi}{\partial x^\beta}+x^\beta\,\frac{\partial \Phi}{\partial x^\alpha}\right]= -\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,x^\alpha\,\frac{\partial \Phi}{\partial x^\beta}\,, \end{align} where we used that \(\dot{v}^\alpha = \vec{g}^\alpha(\vec{x}) = -\partial \Phi/\partial x^\alpha\) and, in the second equality, that the middle expression is symmetric under \(\alpha \leftrightarrow\beta\). With these expressions for the components of the kinetic energy tensor and the potential energy tensor, we can write the tensor virial theorem for a collisionless system as \begin{equation} \frac{\mathrm{d}\vec{G}_{\alpha\beta}}{\mathrm{d}t} = 2\,\vec{K}_{\alpha\beta} +\vec{W}_{\alpha\beta} = 2\,\vec{T}_{\alpha\beta} +\vec{\Pi}_{\alpha\beta}+\vec{W}_{\alpha\beta}\,.\label{eq-tensorvirial-cont-full} \end{equation}

For an equilibrium, collisionless system, we obtain \begin{equation} 2\,\vec{K}_{\alpha\beta} +\vec{W}_{\alpha\beta} = 2\,\vec{T}_{\alpha\beta} +\vec{\Pi}_{\alpha\beta}+\vec{W}_{\alpha\beta} = 0\,.\label{eq-tensorvirial-cont-equil} \end{equation} Taking the trace, we obtain the scalar virial theorem for a continuous, collisionless system \begin{equation} 2\,K + W = 2\,T +\Pi+W = 0\,,\label{eq-tensorvirial-cont-equil-scalar} \end{equation} where the kinetic terms are \begin{align}\label{eq-tensorvirial-T} T & = \frac{1}{2}\,\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\bar{\vec{v}}\cdot \bar{\vec{v}}\,;\quad \Pi = \phantom{\frac{1}{2}\,}\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\overline{(\vec{v}-\bar{\vec{v}})\cdot\,(\vec{v}-\bar{\vec{v}})}\,, \end{align} and they combine to the scalar kinetic energy \begin{align} K = T + \frac{\Pi}{2} & = \frac{1}{2}\,\int\mathrm{d}\vec{x}\,\,\rho(\vec{x},t)\,\overline{v^2}= \frac{M\,\langle v^2\rangle}{2}\,, \end{align} where \(v = |\vec{v}|\) and the \(\langle \cdot \rangle\) averaging is over all dimensions of phase space (rather than just the velocity averaging of \(\bar{\cdot}\)); \(M\) is the total mass of the system.

The scalar potential energy of the system is \begin{equation}\label{eq-tensorvirial-scalar-pot-energy} W = -\int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\vec{x}\cdot \nabla \Phi\,. \end{equation}

By using the Poisson equation, we can obtain an alternative form as (dropping the [\(\vec{x},t\)] arguments to reduce clutter) \begin{align} W & = -\frac{1}{4\pi G} \int\mathrm{d}\vec{x}\,\nabla^2 \Phi\,\vec{x}\cdot \nabla \Phi\,,\\ & = -\frac{1}{4\pi G} \sum_{i} \left[\int\mathrm{d}\vec{x}\,\frac{\partial^2 \Phi}{\partial x_i^2}\,x_i\,\frac{\partial \Phi}{\partial x_i} + \sum_{j\neq i}\int\mathrm{d}\vec{x}\,\frac{\partial^2 \Phi}{\partial x_j^2}\,x_i\,\frac{\partial \Phi}{\partial x_i}\right]\,,\\ & = -\frac{1}{4\pi G} \sum_{i} \left[\int\mathrm{d}\vec{x}\,\frac{\partial^2 \Phi}{\partial x_i^2}\,x_i\,\frac{\partial \Phi}{\partial x_i} - \sum_{j\neq i}\int\mathrm{d}\vec{x}\,\frac{\partial \Phi}{\partial x_j}\,x_i\,\frac{\partial^2 \Phi}{\partial x_j \partial x_i}\right]\,,\\ & = -\frac{1}{8\pi G} \sum_{i} \left[ \int\mathrm{d}\vec{x}\,\frac{\partial }{\partial x_i}\left(\left|\frac{\partial \Phi}{\partial x_i}\right|^2\right)\,x_i - \sum_{j\neq i} \int\mathrm{d}\vec{x}\,\frac{\partial }{\partial x_i}\left(\left|\frac{\partial \Phi}{\partial x_j}\right|^2\right)\,x_i\right]\,, \end{align} where we used integration-by-parts to re-write the second term. After another round of integration by parts we get \begin{align} 8\pi G\,W & = \sum_{i} \left[ \int\mathrm{d}\vec{x}\,\left|\frac{\partial \Phi}{\partial x_i}\right|^2 - \sum_{j\neq i} \int\mathrm{d}\vec{x}\,\left|\frac{\partial \Phi}{\partial x_j}\right|^2\right]=- \sum_{i} \int\mathrm{d}\vec{x}\,\left|\frac{\partial \Phi}{\partial x_i}\right|^2 = \sum_{i} \int\mathrm{d}\vec{x}\,\frac{\partial^2 \Phi}{\partial x_i^2}\,\Phi \end{align} because there are twice as many \(j \neq i\) as there are \(j = i\) in the sums in the second expression. After applying the Poisson equation again, we find the following simple expression: \begin{align}\label{eq-tensorvirial-scalar-pot-energy-dens-pot} W & = \frac{1}{2} \int\mathrm{d}\vec{x}\,\rho(\vec{x},t)\,\Phi(\vec{x},t)\,. \end{align}

By considering the amount of work necessary to move a small amount of mass from infinity to position \(\vec{x}\), it can be shown that this scalar potential energy is the total work necessary to assemble the mass distribution \(\rho(\vec{x})\).

14.1.3. Properties and examples of the potential energy tensor

\label{sec-tensorvirial-potenergytensor}

Because we will be using the potential energy tensor and its scalar form later in this chapter and also in later chapters, it is worth discussing it in a little more detail. In this section, we discuss some general properties of the potential energy tensor and we calculate it for a few example cases that we will use later.

For a spherical mass distribution, the potential energy tensor is particularly simple. The general expression from Equation (14.15) for a spherical mass distribution becomes \begin{align} \vec{W}_{\alpha\beta} & = -\int\mathrm{d}\vec{x}\,\rho(r,t)\,x^\alpha\,\frac{\partial \Phi}{\partial r}\,\frac{\partial r}{\partial x^\beta} = -\int\mathrm{d}r\,r\,\rho(r,t)\,\frac{\partial \Phi}{\partial r}\,\int \mathrm{d}\theta\,\mathrm{d}\phi\,\sin\theta\,x^\alpha\,x^\beta\,. \end{align} From the definition of spherical coordinates in Equations (A.1, it is straightforward to compute that the inner integral \(\int \mathrm{d}\theta\,\mathrm{d}\phi\,\sin\theta\,x^\alpha\,x^\beta = -4\pi\,r^2\delta^{\alpha\beta}/3\), where \(\delta^{\alpha\beta}\) is the Kronecker delta, and therefore \begin{align}\label{eq-potenergytensor-sphere} \vec{W}_{\alpha\beta} & = \frac{4\pi}{3}\,\delta^{\alpha\beta}\int\mathrm{d}r\,r^3\,\rho(r,t)\,\frac{\partial \Phi}{\partial r}= -\frac{4\pi\,G}{3}\,\delta^{\alpha\beta}\int\mathrm{d}r\,r\,\rho(r,t)\,M(< r,t)\,, \end{align} where we used the relation between the radial component of the gravitational field and the enclosed mass (Equation 2.22). Thus, the potential energy tensor for a spherical mass distribution is diagonal and furthermore, because the expression does not depend on the \(\alpha\) and \(\beta\) indices except for the Kronecker delta, it is isotropic with all elements on the diagonal equal. The scalar potential energy is therefore \begin{align} W & = -4\pi\,G\int\mathrm{d}r\,r\,\rho(r,t)\,M(< r)\,. \end{align}

An example that we will use in Chapter 17 is the potential energy of a homogeneous sphere (see Chapter 2.4.2). Because the enclosed mass for the homogeneous sphere is \(M(<r) = 4\pi\,\rho_0\,r^3/3\) and the density is zero at \(r \geq R\), we have \begin{align} W^{\mathrm{homog.\ sphere}} & = -4\pi\,G\int_0^R\mathrm{d}r\,r\,\rho_0\,\frac{4\pi\,\rho_0\,r^3}{3} = -\frac{3}{5}\,\frac{GM^2}{R}\,.\label{eq-potenergytensor-homogsphere} \end{align} where \(M = 4\pi\,\rho_0\,R^3/3\) is the total mass.

Another example is that of a thin spherical shell, which is useful to compare to the thin oblate shell that we will discuss below. For a thin shell of mass \(M\) at radius \(R\), which we first encountered in Chapter 2.2, we have that \(\rho(r) = M\delta(r-R)/[4\pi\,R^2]\) and \(\partial \Phi/\partial r = GM\,\Theta(r-R)/r^2\) where \(\Theta(r)\) is the Heaviside function. Thus, Equation (14.29) gives \begin{align} \vec{W}^{\mathrm{oblate\ shell}}_{\alpha\beta} & = \frac{GM^2}{3\,R^2}\,\delta^{\alpha\beta}\int\mathrm{d}r\,r\,\delta(r-R)\,\Theta(r-R) = \frac{GM^2}{6\,R}\,\delta^{\alpha\beta}\,,\label{eq-potenergytensor-sphericalshell} \end{align} because \(\int\mathrm{d}r\,r\,\delta(r-R)\,\Theta(r-R) = R/2\) (see Equation B.33).

For an axisymmetric mass distribution symmetric with respect to rotation around \(z\), the potential energy tensor similarly simplifies. The general expression from Equation (14.15) now becomes \begin{align} \vec{W}_{\alpha\beta} & = -\int\mathrm{d}\vec{x}\,\rho(R,z,t)\,x^\alpha\,\begin{cases} \frac{\partial \Phi}{\partial R}\,\frac{x^\beta}{R}\,, & \beta \in \{x,y\}\\ \frac{\partial \Phi}{\partial z}\,, & \beta = z\end{cases}\,,\\ & = -\delta^{\alpha\beta}\int\mathrm{d}R\,\mathrm{d}z\,\rho(R,z,t)\,\begin{cases} \pi R^2\,\frac{\partial \Phi}{\partial R}\,, & \beta \in \{x,y\}\\ 2\pi R z\,\frac{\partial \Phi}{\partial z}\,, & \beta = z\end{cases}\,.\label{eq-potenergy-axi} \end{align} Thus, for an axisymmetric mass distribution, the potential energy tensor is diagonal as well and \(W_{xx} = W_{yy}\).

To better understand the potential energy tensor for non-spherical systems, we can write \(\vec{W}_{\alpha\beta}\) in a form similar to that in Equation (14.5) (for simplicity, we drop the time dependence in the following equations) \begin{equation}\label{eq-pot-energy-tensor-likecollisional} \vec{W}_{\alpha\beta} = -\frac{1}{2}\,G\int\mathrm{d}\vec{x}\,\mathrm{d}\vec{x}'\,\rho(\vec{x})\,\rho(\vec{x}')\,\frac{(x'^\alpha-x^\alpha)(x'^\beta-x^\beta)}{|\vec{x}'-\vec{x}|^3}\,, \end{equation} which is the continuous limit of Equation (14.5) and can be derived from Equation (14.15) by using Equation (2.12). Because the numerator of the final factor is the only place that the indices \(\alpha\) and \(\beta\) appear, this expression makes it clear that the relative magnitude of the diagonal elements is set by the distribution of \(|x'^\alpha-x^\alpha|\) of the mass distribution. For a system that is flattened along the \(z\) axis, the typical magnitude of \(|z'-z|\) is less than that of \(|x'-x'|\) and \(|y'-y|\) and therefore \(W_{zz}/W_{xx} < 1\) and \(W_{zz}/W_{yy} < 1\). For an oblate, axisymmetric system, \(W_{zz}/W_{xx} = W_{zz}/W_{yy} < 1\).

To understand the dynamical structure of elliptical galaxies, it is useful to approximate them as oblate spheroids (see Chapter 12.2), where the density is only a function of the coordinate \(m\) defined through Equation (12.8): \(m^2 = R^2 + z^2/q^2\) with \(q \leq 1\). As discussed in Chapter 12.2, this means that the density in \((R,z)\) is constant on the ellipses with eccentricity \(e = \sqrt{1-q^2}\). Because this model is axisymmetric, the potential-energy tensor is diagonal and its diagonal components can be calculated using Equation (14.34) above. This gives \begin{align}\label{eq-potenergytensor-oblateshell} W_{xx} & = -\frac{GM^2\,A_x(q)}{4\,m_0\,q}\,;\quad W_{zz} = -\frac{GM^2\,q\,A_z(q)}{4\,m_0}\,, \end{align} where \(M\) is the mass of the shell, \(m_0\) is its spheroidal radius, and \(A_i(q)\) are functions that are typically expressed as a function of the eccentricity \(e = \sqrt{1-q^2}\) of the shell and which are given by \begin{align} A_x(e) & = \phantom{2}\frac{\sqrt{1-e^2}}{e^2}\,\left(\frac{\sin^{-1} e}{e}-\sqrt{1-e^2}\right)\,;\quad A_z(e) = 2\frac{\sqrt{1-e^2}}{e^2}\,\left(\frac{1}{\sqrt{1-e^2}}-\frac{\sin^{-1} e}{e}\right)\,. \end{align} Because \(A_x(0) = A_z(0) = 2/3\), it is clear that in the limit of a spherical shell, Equations (14.36) reduce to the corresponding expression of Equation (14.32) for a spherical shell.

The ratio of the \(zz\) and \(xx\) components of the potential-energy tensor for a thin oblate spheroidal shell is then given by (we express all occurrences of the axis ratio \(q\) in terms of the eccentricity) \begin{equation}\label{eq-Wzz-over-Wxx-oblateshell} \frac{W_{zz}}{W_{xx}} = 2\,(1-e^2)\,\frac{\left(\frac{1}{\sqrt{1-e^2}}-\frac{\sin^{-1} e}{e}\right)}{\left(\frac{\sin^{-1}e}{e}-\sqrt{1-e^2}\right)}\,. \end{equation} The functional behavior of this ratio for eccentricities between 0 (spherical) and 1 (razor-thin) is as shown in Figure 14.1.

def Wzz_over_Wxx(e):
    return 2.*(1.-e**2.)*(1./numpy.sqrt(1.-e**2.)-numpy.arcsin(e)/e)/(numpy.arcsin(e)/e-numpy.sqrt(1.-e**2.))
# To better resolve the behavior near e=1, use a denser grid there
e_WzzWxx= numpy.hstack((numpy.linspace(1e-4,0.5,51),
                     sorted(1.-numpy.geomspace(1e-4,0.5,51))))
figure(figsize=(6,4))
plot(e_WzzWxx,Wzz_over_Wxx(e_WzzWxx),label=r'$\mathrm{Equation\  (14.38)}$')
# Direct calculation discussed in text below for Hernquist and PerfectEllipsoid
from scipy import integrate
from galpy.potential import (evaluateRforces, evaluatezforces,
                             evaluateDensities, TriaxialHernquistPotential,
                             PerfectEllipsoidPotential)
def potenergy_tensor_axi(pot,component='xx',Rmin=0.,Rmax=numpy.inf,zmin=lambda R: 0., zmax= lambda R: numpy.inf):
    if component[0] != component[1]:
        return 0.
    if component == 'xx' or component == 'yy':
        def integrand(z,R):
            return evaluateDensities(pot,R,z,phi=0.)*R**2.*evaluateRforces(pot,R,z,phi=0.)
    elif component == 'zz':
        def integrand(z,R):
            return 2.*evaluateDensities(pot,R,z,phi=0.)*R*z*evaluatezforces(pot,R,z,phi=0.)
    else:
        return RuntimeError(f'Component {component} not understood')
    return 2.*numpy.pi*integrate.dblquad(integrand,Rmin,Rmax,zmin,zmax)[0]
def q_from_e(e):
    return numpy.sqrt(1.-e**2.)
def e_from_q(q):
    return numpy.sqrt(1.-q**2.)
# Hernquists
# a=3
es_hernq_a3= numpy.arange(0.1,1.0,0.1)
a= 3.
Wxx_hernq_a3= numpy.zeros_like(es_hernq_a3)
Wzz_hernq_a3= numpy.zeros_like(es_hernq_a3)
for ii,e in enumerate(es_hernq_a3):
    thp= TriaxialHernquistPotential(a=a,c=q_from_e(e))
    Wxx_hernq_a3[ii]= potenergy_tensor_axi(thp,component='xx')
    Wzz_hernq_a3[ii]= potenergy_tensor_axi(thp,component='zz')
plot(es_hernq_a3,Wzz_hernq_a3/Wxx_hernq_a3,'o',label=r'$\mathrm{Hernquist},\ a=3$')
# a= 6
es_hernq_a6= numpy.arange(0.05,1.05,0.1)
a= 6.
Wxx_hernq_a6= numpy.zeros_like(es_hernq_a6)
Wzz_hernq_a6= numpy.zeros_like(es_hernq_a6)
for ii,e in enumerate(es_hernq_a6):
    thp= TriaxialHernquistPotential(a=a,c=q_from_e(e))
    Wxx_hernq_a6[ii]= potenergy_tensor_axi(thp,component='xx')
    Wzz_hernq_a6[ii]= potenergy_tensor_axi(thp,component='zz')
plot(es_hernq_a6,Wzz_hernq_a6/Wxx_hernq_a6,'d',label=r'$\mathrm{Hernquist},\ a=6$')
# PerfectEllipsoid
es_pellips_a4= numpy.arange(0.075,1.075,0.1)
a= 4.
Wxx_pellips_a4= numpy.zeros_like(es_pellips_a4)
Wzz_pellips_a4= numpy.zeros_like(es_pellips_a4)
for ii,e in enumerate(es_pellips_a4):
    thp= PerfectEllipsoidPotential(a=a,c=q_from_e(e))
    Wxx_pellips_a4[ii]= potenergy_tensor_axi(thp,component='xx')
    Wzz_pellips_a4[ii]= potenergy_tensor_axi(thp,component='zz')
plot(es_pellips_a4,Wzz_pellips_a4/Wxx_pellips_a4,'p',label=r'$\mathrm{Perfect\ ellipsoid}$')
xlabel(r'$e$')
ylabel(r'$W_{zz}/W_{xx}$')
ylim(0.0,1.05)
legend(frameon=False,fontsize=18.);

Figure 14.1: Ratio of the \(zz\) and \(xx\) components of the potential-energy tensor for a thin oblate spheroidal shell from Equation (14.38) compared to a direct calculation for different oblate potentials.

As expected, the ratio is equal to one for a spherical shell and smaller than one for flattened shells.

Equation (14.38) furthermore demonstrates that the ratio \(W_{zz}/W_{xx}\) is only a function of the axis ratio of the oblate shell, not of its radius or mass. Another way of stating this is to say that the ratio \(W_{zz}/W_{xx}\) depends only on the shape of the oblate shell, but not on its radial structure. Remarkably, this result continues to hold for an extended density distribution that is constant on similar spheroids, as was shown by Roberts (1962) (indeed, Roberts 1962 actually showed that this holds for any extended density distribution that is constant on similar ellipsoids). Thus, the ratio \(W_{zz}/W_{xx}\) of any oblate density distribution that is constant on similar spheroids is given by Equation (14.38).

We can test that this is the case by computing the potential energy tensor for different, non-trivial oblate density distributions. First, we implement the calculation of the potential-energy tensor for an axisymmetric distribution using Equation (14.34):

from scipy import integrate
from galpy.potential import evaluateRforces, evaluatezforces, evaluateDensities
def potenergy_tensor_axi(pot,component='xx',Rmin=0.,Rmax=numpy.inf,zmin=lambda R: 0., zmax= lambda R: numpy.inf):
    if component[0] != component[1]:
        return 0.
    if component == 'xx' or component == 'yy':
        def integrand(z,R):
            return evaluateDensities(pot,R,z,phi=0.)*R**2.*evaluateRforces(pot,R,z,phi=0.)
    elif component == 'zz':
        def integrand(z,R):
            return 2.*evaluateDensities(pot,R,z,phi=0.)*R*z*evaluatezforces(pot,R,z,phi=0.)
    else:
        return RuntimeError(f'Component {component} not understood')
    return 2.*numpy.pi*integrate.dblquad(integrand,Rmin,Rmax,zmin,zmax)[0]

While for density distributions without sharp edges, we will generally integrate over all of space, we have allowed the option to specify the integration range so we can test this function against the expression for the potential-energy tensor of the homogeneous sphere in Equation (14.31).

from galpy.potential import HomogeneousSpherePotential
Rmax= 3.
hp= HomogeneousSpherePotential(amp=2.,R=Rmax)
print(f'Wxx = {potenergy_tensor_axi(hp,component='xx',Rmax=Rmax,zmax=lambda R: numpy.sqrt(Rmax**2.-R**2.)):.3f}, '
      f'Wzz = {potenergy_tensor_axi(hp,component='zz',Rmax=Rmax,zmax=lambda R: numpy.sqrt(Rmax**2.-R**2.)):.3f}, '
      f'analytic Wxx = Wzz = {-3./5.*hp.mass(Rmax)**2./Rmax/3.:.3f}'
)

Wxx = -777.600, Wzz = -777.600, analytic Wxx = Wzz = -777.600

Thus, we see that we get the correct result for the homogeneous sphere.

Next, we compute the \(W_{xx}\) and the ratio \(W_{zz}/W_{xx}\) for flattened Hernquist profiles (for which the density \(\rho(m)\) is that of the Hernquist sphere of Chapter 2.4.6). First, we define some utility functions to go back and forth between \(q\) and \(e\).

def q_from_e(e):
    return numpy.sqrt(1.-e**2.)
def e_from_q(q):
    return numpy.sqrt(1.-q**2.)

and then we compute the tensor for \(e = 0.1\) to \(e=0.9\) in \(\delta e = 0.1\) increments. We can then compare these values to what Equation (14.38) predicts. This is shown as the orange markers in Figure 14.1. We see there that they agree very well! When we change the radial density structure by choosing a different scale parameter \(a\) (choosing a different grid in \(e\) so they do not fully overlap), we get the green diamonds in Figure 14.1. All of the points for the different models agree with the analytic prediction. Showing instead the scalar potential energy (the trace of the tensor, \(W = 2W_{xx}+W_{zz}\) in the axisymmetric case) as a function of \(e\) for the different models, we see that they are different in Figure 14.2.

figure(figsize=(6,4))
plot(es_hernq_a3,2.*Wxx_hernq_a3+Wzz_hernq_a3,'o',label=r'$\mathrm{Hernquist},\ a=3$')
plot(es_hernq_a6,2.*Wxx_hernq_a6+Wzz_hernq_a6,'d',label=r'$\mathrm{Hernquist},\ a=6$')
xlabel(r'$e$')
ylabel(r'$W$')
xlim(-0.03,1.03)
legend(frameon=False,fontsize=18.);

Figure 14.2: Scalar potential energy for flattened Hernquist profiles with different scale parameters.

To completely change the radial structure, we can also consider a completely different density profile, that of the perfect ellipsoid. The ratio \(W_{zz}/W_{xx}\) in that case is shown as the red pentagons in Figure 14.1. We again find good agreement, further demonstrating that the ratio \(W_{zz}/W_{xx}\) is independent of the radial structure of the oblate mass model and only depends on the distribution’s shape, that is, its axis ratio. In the next section we will exploit this property to better understand the dynamical structure of elliptical galaxies.

← Previous section Next section →