5.4. The Jeans equations¶
5.4.1. Generalities¶
As a partial differential equation, Equation (5.33) admits a wide range of solutions depending on the initial and boundary conditions of the problem being considered. We will discuss two applications of solving the collisionless Boltzmann equation here: (i) to find an equilibrium distribution \(f(\vec{x},\vec{v})\) for a given potential \(\Phi(\vec{x})\) and a given density \(\nu(\vec{x})\) (which may or may not source the potential through the Poisson equation; we write the density as \(\nu(\vec{x})\) in case this density is not the full mass density \(\rho(\vec{x})\) that sources the potential) and (ii) to measure the mass distribution \(\Phi(\vec{x})\) from observations of \((\vec{x},\vec{v})\) assuming equilibrium.
In the next section, we will discuss direct solutions of \(f(\vec{x},\vec{v})\) of the collisionless Boltzmann equation, but in this section we first discuss a different and somewhat simpler approach. Rather than working with the collisionless Boltzmann equation (5.29) itself, we can multiply this equation with \(\vec{x}^\alpha\) or \(\vec{v}^\beta\) for some \((\alpha,\beta)\) and integrate over a part of phase-space to obtain moment equations which are known as Jeans equations.
The Jeans equations involve moments of the distribution function. Because astronomical observations are typically performed at a fixed \(\vec{x}\) (or at least at a fixed two-dimensional location on the sky), multiplying the collisionless Boltzmann equation by \(\vec{v}^\beta\) and integrating over all components of the velocity most directly connects to observable quantities. Typical moments that we will consider are the spatial number density \(\nu(\vec{x})\) \begin{equation} \nu(\vec{x}) = \int \mathrm{d}\vec{v}\,f(\vec{x},\vec{v})\,, \end{equation} and the mean velocity \(\bar{\vec{v}}(\vec{x})\) \begin{equation} \bar{\vec{v}}(\vec{x}) = \frac{1}{\nu(\vec{x})}\,\int \mathrm{d}\vec{v}\,\vec{v}\,f(\vec{x},\vec{v})\,, \end{equation} which has components \(\bar{v}_i\). We also consider the velocity dispersion tensor \(\vec{\sigma}(\vec{x})\) with components \(\sigma_{ij}(\vec{x})\) \begin{equation} \sigma_{ij}(\vec{x}) = \frac{1}{\nu(\vec{x})}\,\int \mathrm{d}\vec{v}\,(v_i-\bar{v}_i)\,(v_j-\bar{v}_j)\,f(\vec{x},\vec{v})\,. \end{equation} Higher-order moments are defined in a similar manner and quickly become much more complex.
By multiplying Equation (5.30) (or its other forms) by \(\vec{v}^\beta\) and integrating over \(\vec{v}\), we can derive a set of Jeans equations that relate these moments. For example, integrating Equation (5.30) over \(\vec{v}\) we obtain \begin{equation} \int \mathrm{d}\vec{v}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial t} +\int \mathrm{d}\vec{v}\,\dot{\vec{x}}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial \vec{x}} -\frac{\partial \Phi}{\partial \vec{x}}\,\int \mathrm{d}\vec{v}\,\frac{\partial f(\vec{x},\vec{v},t)}{\partial \vec{v}}=0\,. \end{equation} The final term vanishes after applying the divergence theorem from Equation (B.5) (assuming that \(f\) goes to zero as \(\vec{v} \rightarrow \infty\)) and the partial derivatives in the first two terms may be taken outside of the integral, such that this equation becomes \begin{equation} {\partial \nu(\vec{x}) \over \partial t} + \nabla \cdot \left[\nu(\vec{x})\,\bar{\vec{v}}(\vec{x})\right] = 0\,, \end{equation} which is the continuity equation for the density. Thus, we see that this Jeans equation involves both the density \(\nu(\vec{x})\) and the three-dimensional mean velocity \(\bar{\vec{v}}(\vec{x})\). This single equation therefore does not provide enough information to solve for the mean velocity given the density. Similarly multiplying Equation (5.30) by a component of the velocity \(v_j\), integrating over all \(\vec{v}\), and subtracting the result from \(\bar{v}_j\) times the continuity equation leads to the following equation \begin{equation} \nu\,{\partial \bar{v}_j \over \partial t} + \nu\,\bar{v}_i\,{\partial \bar{v}_j \over \partial x_i} +\nu\,{\partial \Phi \over \partial x_j} + {\partial \left[ \nu\sigma^2_{ij}\right] \over \partial x_i} = 0\,, \end{equation} for all \(i,j\) (we have dropped the functional dependence on \(\vec{x}\) of all of the quantities involved for notational simplicity). This Jeans equation therefore relates the density, mean velocity, and velocity dispersion tensor (as well as the gravitational field). However, while we now have an additional three equations, we have introduced the six components of the velocity dispersion tensor and there are therefore still not enough equations to solve for all of the kinematic quantities. Continuing this procedure by multiplying by higher powers of velocity (e.g., \(v_i\,v_j\)) and integrating over velocity leads to Jeans equations that involve higher and higher moments of the distribution function and at no point does the system of equations close, i.e., provide enough equations to uniquely solve for all kinematic properties. It is therefore not possible to solve for all of the moments of the distribution function through the Jeans equation starting from a known \(\Phi(\vec{x})\) and \(\nu(\vec{x})\) without additional assumptions to close the set of equations. This is a physical impossibility, not just a mathematical one: there are infinitely many equilibrium distribution functions in any gravitational potential and this continues to be the case when we demand that the distribution function sources the gravitational potential sources through Equation (5.32) (except for pathological cases such as the point-mass potential).
The Jeans equations are nevertheless useful, because they involve quantities that are readily observable, such as the spatial density (which can be traced through the surface brightness or through star counts), the mean velocity, and the velocity dispersion. By combining these observables with assumptions about the distribution function to close the system of equations, we can therefore constrain \(\Phi(\vec{x})\). In the Milky Way in particular, it is (in principle) possible to measure all moments of the velocity distribution (because we can obtain full six-dimensional phase-space information for many stars) and in this case the Jeans equation leads to direct measurements of \(\Phi(\vec{x})\) without any closure assumptions.
5.4.2. The spherical Jeans equation¶
To demonstrate the usefulness of the Jeans equations, we will start by considering them for spherical systems. To derive the spherical Jeans equation, we start from Equation (5.33) and we use the Hamiltonian appropriate for spherical coordinates \(\vec{q} = (r,\phi,\theta)\) (we also set \(m=1\)). This Hamiltonian has momenta \(\vec{p}\) that are easily derived by writing down the Lagrangian in spherical coordinates: \begin{align}\label{eq-spherequil-generalized-momenta-spherical} p_r & = \dot{r} = v_r\,;\quad p_\phi = r^2\,\sin^2 \theta\,\dot{\phi} = r\,\sin \theta \,v_\phi\,;\quad p_\theta = r^2\,\dot{\theta} = r\,v_\theta\,. \end{align} The Hamiltonian itself is \begin{equation} H = \frac{1}{2}\,\left(p_r^2 + \frac{p_\theta^2}{r^2}+\frac{p_\phi^2}{r^2\,\sin^2\theta}\right) + \Phi(r)\,. \end{equation} Using \(\dot{\vec{q}} = \partial H / \partial \vec{p}\) and \(\dot{\vec{p}} = -\partial H / \partial \vec{q}\) we find that Equation (5.29) becomes \begin{equation} p_r\,\frac{\partial f}{\partial r} + \frac{p_\theta}{r^2}\,\frac{\partial f}{\partial \theta} + \frac{p_\phi}{r^2\,\sin^2\theta}\,\frac{\partial f}{\partial \phi} -\left(\frac{\mathrm{d} \Phi}{\mathrm{d} r}-\frac{p_\theta^2}{r^3}-\frac{p_\phi^2}{r^3\,\sin^2\theta}\right)\,\frac{\partial f}{\partial p_r} +\frac{p_\phi^2\,\cos\theta}{r^2\,\sin^3\theta}\,\frac{\partial f}{\partial p_\theta} = 0\,.\\ \end{equation} We now multiply this by \(p_r\) and integrate over all \((p_r,p_\phi,p_\theta)\) using that \(\mathrm{d}p_r\,\mathrm{d}p_\phi\,\mathrm{d}p_\theta = r^2\,\sin\theta\,\mathrm{d}v_r\, \mathrm{d}v_\phi\,\mathrm{d}v_\theta\) and using partial integration to deal with the derivatives of \(f\) with respect to the momenta. In bringing derivatives outside of integrals, it is important to remember the quantities held constant for each partial derivative: the partial derivatives are for the canonical set of coordinates \((r,\phi,\theta,p_r,p_{\phi},p_{\theta})\) and a derivative like \(\partial / \partial r\) therefore in particular keeps the momenta constant and can therefore be brought outside of the momentum integral. In the end, we obtain \begin{align}\label{eq-spher-jeans-penult} \frac{\partial (r^2\,\sin\theta\,\nu\,\overline{v^2_r})}{\partial r} + \frac{\partial (\sin\theta\,\nu\,\overline{v_r\,v_\theta})}{\partial \theta} + \frac{\partial (\nu\,\overline{v_r\,v_\phi}/\sin\theta)}{\partial \phi} & +r^2\,\sin\theta\,\nu\,\left(\frac{\mathrm{d} \Phi}{\mathrm{d} r}-\frac{\overline{v_\theta^2}}{r}-\frac{\overline{v_\phi^2}}{r}\right) = 0\,. \end{align}
As we will see below, for a spherical system the equilibrium distribution function can only be a function of energy and angular momentum, of which only the energy depends on \(p_r\) and it does so in a quadratic manner. This implies that all odd moments of \(v_r\) must vanish, in particular \(\overline{v_r\,v_\theta} = \overline{v_r\,v_\phi} = 0\) and we can therefore simplify Equation (5.44) to \begin{equation}\label{eq-spher-jeans} \frac{\mathrm{d} (\nu\,\overline{v^2_r})}{\mathrm{d} r} +\,\nu\,\left(\frac{\mathrm{d} \Phi}{\mathrm{d} r}+\frac{2\overline{v_r^2}-\overline{v_\theta^2}-\overline{v_\phi^2}}{r}\right) = 0\,. \end{equation} We can define a parameter \(\beta\) that quantifies the amount of orbital anisotropy (the difference in the velocity distribution in different directions) \begin{align}\label{eq-anisotropy} \beta & \equiv 1 - \frac{\sigma_\theta^2 + \sigma_\phi^2}{2\,\sigma_r^2} = 1 - \frac{\overline{v_\theta^2} + \overline{v_\phi^2}}{2\,\overline{v_r^2}}\,, \end{align} where the first equation is the actual definition and the second equation follows if we assume that the system has no net rotation. This parameter measures a system’s radial anisotropy in the following sense: if the spread in the motions in the radial direction \(r\) is the same as that in the tangential directions \((\phi,\theta)\) then \(\sigma_r = \sigma_\phi = \sigma_\theta\), \(\beta = 0\), and the system is isotropic; if the spread in radial motions is much larger than that in the perpendicular dimensions, \(\sigma_r \gg \sigma_\theta, \sigma_\phi\) and \(\beta \rightarrow 1\) and we say that the system is radially biased; if all orbits are circular, \(\sigma_r = 0\), \(\beta = -\infty\), and we say that the system is tangentially biased (note the unfortunate asymmetry inherent in this definition: the radially biased range for \(\beta\) is \(0 < \beta \leq 1\), while the tangentially biased range is \(-\infty \leq \beta < 0\)). The anisotropy parameter in a realistic system will in general depend on radius \(\beta = \beta(r)\).
In terms of this anisotropy parameter, we can write the spherical Jeans equation for a non-rotating system as \begin{equation}\label{eq-spher-jeans-anisotropy} \frac{\mathrm{d} (\nu\,\overline{v^2_r})}{\mathrm{d} r} + 2\frac{\beta}{r}\,\nu\,\overline{v^2_r} = -\nu\,\frac{\mathrm{d} \Phi}{\mathrm{d} r}\,. \end{equation}
We can also directly express the enclosed-mass profile in terms of the density \(\nu\), the velocity dispersion \(\sigma_r = \sqrt{\overline{v_r^2}}\), and the anisotropy \(\beta\), because \(G\,M(<r)/r^2 = \mathrm{d}\Phi(r) / \mathrm{d} r\) \begin{equation}\label{eq-spher-jeans-anisotropy-mass} M(< r) = -\frac{r\,\sigma_r^2}{G}\,\left(\frac{\mathrm{d} \ln (\nu\,\sigma_r^2)}{\mathrm{d} \ln r} + 2\,\beta\right)\,. \end{equation} For spherical systems, we can therefore directly measure the enclosed-mass profile from measurements of the radial velocity dispersion and the logarithmic decline of the density times the dispersion squared as long as the anisotropy is known or can be constrained. Unfortunately, for most systems of interest it is difficult to measure the anisotropy, because it requires measurements of at least two of the velocity dimensions, while we typically can only access one through measurements of the line-of-sight velocity (proper motions can only be measured within a few to tens of kpc from the Sun and are therefore inaccessible at the present time for extra-galactic systems or for the far reaches of the Milky Way). This gives rise to the famous mass–anisotropy degeneracy in galactic mass modeling. One way around this is to derive higher-order spherical Jeans equations, which make use of higher-order moments of the velocity distribution (but as we discussed above, these equations do not close and assumptions are always necessary to make the equations work with only the line-of-sight velocity).
When the anisotropy \(\beta\) is assumed to be constant and the system is not rotating, the spherical Jeans equation can be integrated to give \begin{equation}\label{eq-spher-jeans-anisotropy-integrated} \nu\,\sigma^2_r = r^{-2\beta}\,\int_r^\infty\,\mathrm{d}r'\,r'^{2\beta}\,\nu(r')\,\frac{\mathrm{d}\Phi}{\mathrm{d}r}= G\,r^{-2\beta}\,\int_r^\infty\,\mathrm{d}r'\,r'^{2\beta-2}\,\nu(r')\,M(< r')\,. \end{equation} The integration range is such that the solution satisfies the boundary condition that \(\lim_{r\rightarrow\infty}\nu\,\sigma_r^2 = 0\).
If \(\beta = \beta(r)\) and the system is not rotating, we can consider Equation (5.47) to be a differential equation for \(\nu\,\overline{v^2_r} = \nu\,\sigma_r^2\). Differential equations of the type in Equation (5.47) can be solved in terms of an integrating factor (see Equation B.114 in Appendix B.4.2) \begin{equation} \gamma(r) = 2\,\int_{r_0}^r\mathrm{d}r'\,\frac{\beta(r')}{r'}\,, \end{equation} where \(r_0\) is some arbitrary lower limit; the solution is then \begin{equation}\label{eq-spher-jeans-anisotropy-integrated-general} \nu\,\sigma^2_r = \frac{\int_r^\infty\mathrm{d} r\,e^{2\,\int_{r_0}^r\mathrm{d}r'\,\beta(r')/r'}\,\nu\,\frac{\mathrm{d} \Phi}{\mathrm{d} r}}{e^{2\,\int_{r_0}^r\mathrm{d}r'\,\beta(r')/r'}}\,. \end{equation} These solutions for constant or radially-varying \(\beta\) can be used to determine the velocity dispersion for a mass model \(M(<r')\) and an assumed \(\beta(r)\) when the three-dimensional density is known. In external systems we can typically only measure the two-dimensional surface-brightness profile \(I(R)\) or the projected two-dimensional stellar density \(\Sigma(R)\); we discuss how to deal with this in Jeans analyses in Chapter 6.5.
As an example of the spherical Jeans equation, we can compute the velocity dispersion for a dark-matter halo modeled as an NFW profile for different values of the anisotropy. We use a halo with virial mass \(M_\mathrm{vir} = 10^{12}\,M_\odot\) and a concentration of 8.37 (following the Dutton & Macciò 2014 relation from Equation 2.66 as in Chapter 2.4.6). The radial velocity dispersion out to the virial radius is shown in Figure 5.4.
[15]:
from galpy.potential import NFWPotential
from galpy.df import jeans
from matplotlib.cm import winter
from matplotlib.ticker import FuncFormatter
np= NFWPotential(mvir=1,conc=8.37,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
betas= [-5.,-3.,-1.,-0.5,0.,0.5,0.75,0.99]
rs= numpy.geomspace(0.1,np.rvir(H=67.1,overdens=200.,wrtcrit=True).to_value(u.kpc),41)*u.kpc
figure(figsize=(10,5))
for ii,beta in enumerate(betas):
sigs= numpy.array([jeans.sigmar(np,r,beta=beta).to_value(u.km/u.s) for r in rs])
subplot(1,2,1)
galpy_plot.plot(rs.to_value(u.kpc),
sigs,
semilogx=False,
xrange=[0.,270.],
yrange=[0.,200.],
gcf=True,
overplot=beta!=betas[0],
color=winter(ii/(len(betas)-1.)),
xlabel=r'$r\,(\mathrm{kpc})$',
ylabel=r'$\sigma_r\,(\mathrm{km\,s}^{-1})$')
# Bit of fiddling with the label placement
text(255.,sigs[-1]-(beta < 0)*(3-ii)*3.5+(beta <= -3)*5.,
rf'${beta:.2f}$',fontsize=15.,
ha='right')
subplot(1,2,2)
galpy_plot.plot(rs.to_value(u.kpc),
sigs,
loglog=True,
gcf=True,
overplot=beta!=betas[0],
xrange=[0.03,400.],
yrange=[1.1,1450.],
color=winter(ii/(len(betas)-1.)),
xlabel=r'$r\,(\mathrm{kpc})$')
subplot(1,2,1)
text(240.,135.,r'$\beta$',fontsize=18.,
ha='center')
subplot(1,2,2)
axvline(np.a*np._ro,ymin=0.4,ymax=0.85,color='k',ls='--',lw=2.,zorder=0)
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
gca().yaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
tight_layout();
Figure 5.4: The radial velocity dispersion of a Milky-Way-like NFW halo.
We show the radial velocity dispersion both on a linear scale as on a logarithmic scale so we can best study its behavior. We see that the velocity dispersion increases as we go from the virial radius inwards, but it only keeps increasing for \(\beta > 0.5\) and for such radially-anisotropic halos, the radial velocity dispersion increases without bounds. For \(\beta < 0.5\), the radial velocity dispersion decreases as one approaches the center. On the left, we indicate the scale radius of the profile by the dashed line, which shows that the transition between increasing and decreasing velocity dispersion occurs around the scale radius. The behavior of the velocity dispersion can be understood from solving the Jeans equation analytically in the regions \(r \ll a\) and \(r \gg a\). The overall scale of the radial velocity dispersion is higher for higher values of \(\beta\).
Using Equation (5.46), we can also compute the tangential dispersion \(\sigma_t = \sqrt{\sigma_\theta^2+\sigma_\phi^2}\) and this is displayed in Figure 5.5.
[16]:
figure(figsize=(8,5))
for ii,beta in enumerate(betas):
sigs= numpy.sqrt(2.*(1.-beta))\
*numpy.array([jeans.sigmar(np,r,beta=beta).to_value(u.km/u.s) for r in rs])
galpy_plot.plot(rs.to_value(u.kpc),
sigs,
semilogx=False,
xrange=[0.,290.],
yrange=[0.,180.],
gcf=True,
overplot=beta!=betas[0],
color=winter(ii/(len(betas)-1.)),
xlabel=r'$r\,(\mathrm{kpc})$',
ylabel=r'$\sigma_t\,(\mathrm{km\,s}^{-1})$')
plot(rs.to_value(u.kpc),np.vcirc(rs).to_value(u.km/u.s),'k--');
Figure 5.5: The tangential velocity dispersion of a Milky-Way-like NFW halo.
The magnitude of the tangential dispersion behaves opposite that of the radial velocity dispersion, increasing as \(\beta\) decreases, reflecting the increasing importance of tangential motions as \(\beta\) decreases. The dashed curve shows the rotation curve of the NFW profile for the chosen parameters and we see that the tangential dispersion approaches the rotation curve as \(\beta \rightarrow -\infty\). Note that while we talk about “tangential dispersion” here and use the symbol \(\sigma_t\), what we actually calculate is \(\sqrt{\langle v_\theta^2\rangle + \langle v_\phi^2\rangle}\) and the tangential support can be provided either through dispersion or rotation. This makes it clear why the tangential dispersion approaches the circular velocity as \(\beta \rightarrow -\infty\): a purely-anisotropic, spherical system must have all bodies on circular orbits.
For more details on the solution of the spherical Jeans equation for NFW profiles, see Lokas & Mamon (2001).
.
.
.