Editing Stellar dynamics (section)

== Applications and examples ==

Stellar dynamics is primarily used to study the mass distributions within stellar systems and galaxies. Early examples of applying stellar dynamics to clusters include [[Albert Einstein]]'s 1921 paper applying the [[virial theorem]] to spherical star clusters and [[Fritz Zwicky]]'s 1933 paper applying the virial theorem specifically to the [[Coma Cluster]], which was one of the original harbingers of the idea of [[dark matter]] in the universe.<ref>{{Cite journal|last=Einstein|first=Albert |date=2002 |title=A Simple Application of the Newtonian Law of Gravitation to Star Clusters|url=http://alberteinstein.info/vufind1/images/einstein/5-166.tr.pdf|journal=The Collected Papers of Albert Einstein|volume=7|pages=230–233|via=Princeton University Press}}</ref><ref>{{Cite journal|last=Zwicky|first=Fritz|date=2009|title=Republication of: The redshift of extragalactic nebulae |journal=General Relativity and Gravitation|volume=41|issue=1|pages=207–224|doi=10.1007/s10714-008-0707-4 |bibcode=2009GReGr..41..207Z |s2cid=119979381}}</ref> The Jeans equations have been used to understand different observational data of stellar motions in the Milky Way galaxy. For example, [[Jan Oort]] utilized the Jeans equations to determine the average matter density in the vicinity of the solar neighborhood, whereas the concept of asymmetric drift came from studying the Jeans equations in cylindrical coordinates.<ref>{{Cite book|title=Astrophysics for Physicists|last=Choudhuri|first=Arnab Rai |publisher=Cambridge University Press|year=2010|isbn=978-0-521-81553-6|location=New York|pages=213–214}}</ref>

Stellar dynamics also provides insight into the structure of galaxy formation and evolution. Dynamical models and observations are used to study the triaxial structure of elliptical galaxies and suggest that prominent spiral galaxies are created from [[Galaxy merger|galaxy mergers.]]<ref name=":0" /> Stellar dynamical models are also used to study the evolution of active galactic nuclei and their black holes, as well as to estimate the mass distribution of dark matter in galaxies.

[[File:Potential of a thick disk.pdf|thumb|Note the somewhat pointed end of the equal potential in the (R,z) meridional plane of this R0=5z0=1 model]]

=== A unified thick disk potential ===

Consider an oblate potential in cylindrical coordinates 
<math display="block">\begin{align} \Phi(R,z) &  ={G M_0 \over 2z_0} 
 \left[2\sinh^{-1}\!\! Q - \sinh^{-1} \!\!Q_{+} - \sinh^{-1} \!\! Q_{-}\right] \\ 
&={G M_0 \over 2 z_0} \log { (\sqrt{1+  Q^2} + Q )^2 \over 
\left[\sqrt{1+ Q_{+}^2}+ Q_{+}\right] \left[\sqrt{1+Q_{-}^2} + Q_{-} \right]},\\
Q_{\pm} & \equiv {R_0 + \left|~ |z| \pm z_0~  \right| \over R}, \\
Q & \equiv {R_0 + [0, |z| - z_0 ]_\max \over R}, \\
\end{align} </math>
where <math> z_0, R_0</math> are (positive) vertical and radial length scales.  
Despite its complexity, we can easily see some limiting properties of the model.

First we can see the total mass of the system is <math>M_0 </math> because  
<math display="block"> \Phi(R,z) 
\rightarrow {G M_0 \over 2z_0} (2 Q_{-} - Q_{-} -Q_{+}) = -{G M_0 \over R} , 
</math>
when we take the large radii limit 
<math> R \rightarrow \infty, ~|z| \ge z_0, </math>, so that <math> Q = Q_{-}=Q_{+}-{2z_0 \over R} = {|z| + (R_0 - z_0)  \over R} \rightarrow 0.</math>

We can also show that some special cases of this unified potential become the potential of the Kuzmin razor-thin disk, that of the Point mass <math> M_0 </math>, and that of a uniform-Needle mass distribution: 
<math display="block"> \Phi_{KM}(R,z) = -{G M_0 \over \sqrt{ R^2 + (|z|+R_0)^2}}, ~~ z_0=0, </math>
<math display="block"> \Phi_{PT}(R,z) = -{G M_0 \over \sqrt{R^2+z^2}} , ~~ z_0=R_0=0, </math>
<math display="block"> \Phi_{UN}^{R_0=0}(R, z) = {G M_0 \over 2z_0} \left[2\sinh^{-1}\!\! {(0, |z| - z_0 )_\max \over R} - \sinh^{-1} \!\!{z_0 + |z| \over R} - \sinh^{-1} \!\!{\left|~z_0 - |z|~\right| \over R}\right]. </math>

=== A worked example of gravity vector field in a thick disk ===

First consider the vertical gravity at the boundary, 
<math display="block"> g_z(R,z) = - \partial_z \Phi(R,z) = -{G M_0 z \over 2z_0^2} \left[ {1 \over \sqrt{R_0^2+ R^2}} - { 1 \over \sqrt{(R_0+2z_0)^2 + R^2} } \right] , ~~ z= \pm z_0, </math>

Note that both the potential and the vertical gravity are continuous across the boundaries, hence no razor disk at the boundaries.
Thanks to the fact that at the boundary, <math>\partial_{|z|} (2 Q) - \partial_{|z|} Q_{-} = \partial_{|z|} \left(Q_{+} - \frac{2z_0}{R}\right) = {1 \over R} </math> is continuous.  Apply Gauss's theorem by integrating the vertical force over the entire disk upper and lower boundaries, we have 
<math display="block"> 2 \int_0^\infty (2 \pi R dR) |g_z(R,z_0)| = 4 \pi G M_0, </math>
confirming that <math>M_0</math> takes the meaning of the total disk mass.

The vertical gravity drops with <math display="block"> -g_z \rightarrow G M_0 z (1+R_0/z_0)/R^3 </math> at large radii, which is enhanced over the vertical gravity of a point mass <math> G M_0 z/R^3 </math> due to the self-gravity of the thick disk.

=== Density of a thick disk from Poisson Equation ===

Insert in the cylindrical Poisson eq. 
<math display="block"> \rho(R,z) ={\partial_z \partial_z \Phi \over 4 \pi G} + {\partial_R (R\partial_R \Phi) \over 4 \pi G R} = { M_0 R_0/z_0 \over 4\pi (R^2+R_0^2)^{3/2}} H(z_0-|z|), </math>
which drops with radius, and is zero beyond <math> |z| > z_0</math> and uniform along the z-direction within the boundary.

[[File:Density of a thick disk model.pdf|thumb|Note the vertically uniform thick disk density contour in this R0=5z0=1 model]]

=== Surface density and mass of a thick disk ===

Integrating over the entire thick disc of uniform thickness <math> 2 z_0 </math>, we find the surface density and the total mass as 
<math display="block"> \Sigma(R) = (2 z_0)\rho(R,0), ~~ M_0 = \int_0^\infty (2\pi R dR) \Sigma(R). </math>

This confirms that the absence of extra razor thin discs at the boundaries.  In the limit, <math> z_0 \rightarrow 0 </math>, this thick disc potential reduces to that of a razor-thin Kuzmin disk, for which we can verify <math>{|g_z (R,0+)| \over 2\pi G} \rightarrow \Sigma(R) \rightarrow {M_0 R_0 \over 2\pi (R^2+R_0^2)^{3/2}} </math>.

=== Oscillation frequencies in a thick disk ===

To find the vertical and radial oscillation frequencies, we do a Taylor expansion of potential around the midplane.
<math display="block"> \Phi (R_1, z) \approx \Phi(R,0) + {\omega^2 R} (R_1-R) + {\kappa^2 \over 2} (R_1-R)^2 + {\nu^2 \over 2} z^2 </math>
and we find the circular speed <math> V_\text{cir}</math> and the vertical and radial epicycle frequencies to be given by
<math display="block"> (R \omega)^2 \equiv V^2_\text{cir}  = \left[{(1+R_0/z_0) G M_0\over \sqrt{R^2+(R_0+z_0)^2} } - {(R_0/z_0) G M_0 \over \sqrt{R^2+R_0^2}} \right], </math>
<math display="block"> \nu^2 = {G M_0 (R_0/z_0 + 1) \over (R^2+(R_0+z_0)^2)^{3/2}}, </math>
<math display="block"> \kappa^2 + \nu^2 - 2 \omega^2 = 4 \pi G \rho(R,0) = {G M_0 R_0/z_0 \over (R^2+R_0^2)^{3/2}}. </math>
Interestingly the rotation curve <math> V_\text{cir} </math> is solid-body-like near the centre <math> R \ll R_0 </math>, and is Keplerian far away.

At large radii three frequencies satisfy 
<math display="inline"> \left.\left[\omega, \nu, \kappa, \sqrt{4 \pi G \rho}\right]\right|_{R\to \infty} \to [1,1+R_0/z_0,1, R_0/z_0]^{1\over 2} \sqrt{G M_0\over R^3} </math>.
E.g., in the case that <math> R \to \infty </math> and <math> R_0 / z_0 = 3</math>, the oscillations <math> \omega: \nu: \kappa = 1: 2 : 1 </math> forms a resonance.

In the case that <math> R_0 =0 </math>, the density is zero everywhere except uniform needle between <math>|z| \le z_0 </math> along the z-axis.

If we further require <math> z_0=0</math>, then we recover a well-known property for closed ellipse orbits in point mass potential,
<math display="block"> \omega: \nu: \kappa = 1: 1 : 1 .</math>

=== A worked example for neutrinos in galaxies ===

For example, the phase space distribution function of non-relativistic neutrinos of mass m anywhere will not exceed the maximum value set by
<math display="block"> f(\mathbf{x},\mathbf{v},t) = {dN \over dx^3 dv^3} \le {6 \over (2\pi \hbar/m)^3},  ~~~  </math>
where the Fermi-Dirac statistics says there are at most 6 flavours of neutrinos within a volume <math> dx^3 </math> and a velocity volume <math display="block"> dv^3  = (dp/m)^3 = [(2\pi\hbar/dx)/m]^3,</math>.

Let's approximate the distribution is at maximum, i.e., 
<math display="block"> f(x, y, z, V_x, V_y, V_z) = {6 \over (2\pi \hbar/m)^3} q^{\alpha \over 2}, ~~ 0 \le q(E)={\Phi_{\max}- E \over V_0^2/2} \le 1, </math>
where <math> 0 \ge \Phi_{\max} \ge E= \Phi(x,y,z) + {V_x^2 + V_y^2 + V_z^2 \over 2} \ge \Phi_{\min} \equiv \Phi_{\max}- {V_0^2 \over 2} </math> such that <math> E_{\min}, E_{\max} </math>, respectively, is the potential energy of at the centre or the edge of the gravitational bound system.  The corresponding neutrino mass density, assume spherical, would be
<math display="block">\rho(r) = n(x,y,z) m = \int dV_x \int dV_y \int dV_z ~m~ f(x,y,z,V_x,V_y,V_z), </math>
which reduces to
<math display="block">\rho(r) = { C  (\Phi_{\max}-\Phi(r))^{3+\alpha \over 2} \over (\Phi_{\max}-\Phi_{\min} )^{\alpha \over 2} }, ~~~ C={6 m \pi 2^{5/2} B\left(1+{\alpha \over 2}, {3 \over 2}\right) \over (2\pi \hbar/m)^3} </math>

Take the simple case <math> \alpha \to 0 </math>, and estimate the density at the centre <math> r=0 </math> with an escape speed <math> V_0 </math>, we have
<math display="block"> \rho(r) \le \rho(0) \rightarrow { m^4 V_0^3 \over \pi^2 \hbar^3} \approx m_\mathrm{eV}^4 V_{200}^3 \times \text{[Cosmic Critical Density]}.</math>

Clearly eV-scale neutrinos with <math> m_{eV} \sim 0.1-1 </math> is too light to make up the 100–10000 over-density in galaxies with escape velocity <math> V_{200} \equiv V/(\mathrm{200km/s}) \sim 0.1-3.4 </math>, while 
neutrinos in clusters with <math> V \sim \mathrm{2000 km/s} </math> could make up <math> 100-1000 </math> times cosmic background density.

By the way the freeze-out cosmic neutrinos in your room have a non-thermal random momentum <math display="inline"> \sim {(\mathrm{2.7 K}) k \over c} \sim (1~\mathrm{eV}/c^2) (\mathrm{70 km/s}) </math>, and do not follow a Maxwell distribution, and are not in thermal equilibrium with the air molecules because of the extremely low cross-section of neutrino-baryon interactions.

=== A Recap on Harmonic Motions in Uniform Sphere Potential ===

Consider building a steady state model of the fore-mentioned uniform sphere of density <math> \rho_0 </math> and potential <math> \Phi(r)</math>
<math display="block"> \begin{align}
\rho(|\mathbf{r}|) &=\rho_0 \equiv M_\odot n_0, ~~ |\mathbf{r}|^2=x^2+y^2+z^2 \le r_0^2, ~~ \Omega\equiv \sqrt{4 \pi G \rho_0 \over 3} \equiv {V_0 \over r_0} \\
\Phi(|\mathbf{r}|) &= {\Omega^2 (x^2+y^2+z^2) -3 V_0^2 \over 2}= {V_e(r)^2 \over 2} - \Phi(r_0), 
\end{align}</math>
where <math> V_e(r) = V_0 \sqrt{1-{r^2 \over r_0^2}} =\sqrt{2\Phi(r_0)-2\Phi(r)}</math> is the speed to escape to the edge <math> r_0</math>.

First a recap on motion "inside" the uniform sphere potential.
Inside this constant density core region, individual stars go on resonant harmonic oscillations of angular frequency <math> \Omega </math> with <math display="block">  \begin{align}
\ddot{x} = & - \Omega^2 x =-\partial_x \Phi,  \\
\ddot{y} = & - \Omega^2 y, ~~~  {\dot{y}(t)^2 \over 2}+{\Omega^2 y(t)^2 \over 2} \equiv I_y(y,\dot{y}) ={\dot{y}(0)^2 \over 2}+ {\Omega^2 y(0)^2 \over 2} \le {(\Omega r_0)^2 \over 2} \\
\ddot{z} = & - \Omega^2 z, \rightarrow \dot{z}(t)=  \dot{z}(0) \cos (\Omega t) + \Omega z(0) \sin (\Omega t).
\end{align}</math>
Loosely speaking our goal is to put stars on a weighted distribution of orbits with various energies <math> f\left(I_x(x,\dot{x}), I_y(y,\dot{y}), I_z(z,\dot{z}\right) = DF(\mathbf{r},\mathbf{V})</math>, i.e., the phase space density or distribution function, such that their overall stellar number density reproduces the constant core, hence their collective "steady-state" potential.  Once this is reached, we call the system is a self-consistent equilibrium.

=== Example on Jeans theorem and CBE on Uniform Sphere Potential ===

Generally for a time-independent system, Jeans theorem predicts that <math> f(\mathbf{x},\mathbf{v}) </math> is an implicit function of the position and velocity through a functional dependence on "constants of motion".

For the uniform sphere, a solution for the Boltzmann Equation, written in spherical coordinates <math> (r,\theta,\phi) </math> and its velocity components <math> (V_r,V_\theta,V_\phi) </math> is
<math display="block"> f(r,\theta,\varphi,V_r,V_\theta,V_\varphi) = {C_0 \over V_0^3} \sqrt{V_0^2 \over 2Q}, </math>
where <math>C_0 = 2\pi^{-2} \rho_0 </math> is a normalisation constant, which has the dimension of (mass) density. And we define a (positive enthalpy-like dimension <math> \text{km}^2/\text{s}^2 </math>) Quantity 
<math display="block"> Q[\mathbf{x},\mathbf{v}] \equiv \left[0, \left(-V_0^2 - E \right) + {J^2 \over 2 r_0^2} \right]_\max \left[{J_z \over |J_z|}, 0\right]_\max .</math>
Clearly anti-clockwise rotating stars with <math> J_z \le 0,~~ Q=0</math> are excluded.

It is easy to see in spherical coordinates that

<math> J^2 = r^2 V_t^2 = r^2 (V_\theta^2+V_\varphi^2),  </math>

<math> J_z = V_\varphi r \sin\theta, </math>

<math> E = {V_r^2+V_t^2 \over 2} + \Phi(r), ~ V_t \equiv \sqrt{V_\theta^2+V_\varphi^2}</math>

Insert the potential and these definitions of the orbital energy E and angular momentum J and its z-component Jz along every stellar orbit, we have
<math display="block"> 2Q= \text{Heaviside}\left({V_\varphi \over |V_\varphi|}\right) \times  \left[ V_0^2 \left(1-{r^2 \over r_0^2}\right) - V_r^2  - \left(1 - {r^2 \over r_0^2}\right) {\left(V_\theta^2+V_\varphi^2\right)}, 0 \right]_\max, </math> which implies <math> |V_r| \le V_e(r)</math>, and <math> |V_\theta|, V_\varphi </math> between zero and <math> V_0 </math>.

To verify the above <math> E, ~J_z</math> being constants of motion in our spherical potential, we note
<math display="block"> dE/dt = {\partial E\over \partial t} + \mathbf{v} {\partial E \over \partial \mathbf{x}} + (\mathbf{-\nabla \Phi}) {\partial E \over \partial \mathbf{v}} </math>

<math display="block"> dE/dt 
= {\partial \Phi\over \partial t} + \mathbf{v} {\partial \Phi \over \partial \mathbf{x}} + (\mathbf{-\nabla \Phi}) \mathbf{v} = {\partial \Phi\over \partial t} =0 </math> for any "steady state" potential.

<math display="block"> dJ_z/dt
= {\partial J_z\over \partial t} + {\partial J_z \over \partial \mathbf{x}} \cdot \mathbf{v} - (\mathbf{\nabla \Phi}) \cdot {\partial J_z \over \partial \mathbf{v}}, </math> which reduces to <math> dJ_z/dt = 0 + [(V_y)V_x + (-V_x)V_y] - \left[(-y) {x \over R}{\partial \Phi(R,z) \over \partial R} + (x) {y\over R}{\partial \Phi(R,z) \over \partial R}\right]
= 0 </math> around the z-axis of any axisymmetric potential, where <math display="inline">R=\sqrt{x^2+y^2} </math>.

Likewise the x and y components of the angular momentum are also conserved for a spherical potential.  Hence <math> dJ/dt =0 </math>.

So for any time-independent spherical potential (including our uniform sphere model), 
the orbital energy E and angular momentum J and its z-component Jz along every stellar orbit satisfy
<math display="block"> dE[\mathbf{x},\mathbf{v}]/dt = dJ[\mathbf{x},\mathbf{v}]/dt= dJ_z[\mathbf{x},\mathbf{v}]/dt =0 .</math>

Hence using the chain rule, we have
<math display="block"> {d \over dt} Q(E[\mathbf{x},\mathbf{v}],J[\mathbf{x},\mathbf{v}],J_z[\mathbf{x},\mathbf{v}]) = {\partial Q \over \partial E} {dE \over dt} + {\partial Q \over \partial J_z} {dJ_z \over dt} + {\partial Q \over \partial J} {dJ \over dt} = 0 ,</math>
i.e., <math display="inline"> {d \over dt} f= f'(Q) {d Q[\mathbf{x},\mathbf{v}]\over dt}  =0 </math>, so that CBE is satisfied, i.e., our 
<math display="block"> f(\mathbf{x},\mathbf{v}) = f(E[\mathbf{x},\mathbf{v}],J[\mathbf{x},\mathbf{v}],J_z[\mathbf{x},\mathbf{v}]) </math>
is a solution to the Collisionless Boltzmann Equation for our static spherical potential.

=== A worked example on moments of distribution functions in a uniform spherical cluster ===

We can find out various moments of the above distribution function, reformatted as with the help of three Heaviside functions, 
<math display="block"> f(|\mathbf{r}|,V_r,V_\theta,V_\varphi) = 
{C_0 \over V_0^3} \left.{\text{H}(1-x) 
\over \left(1-x^2\right)^{1 \over 2}}\right|_{x \equiv {|\mathbf{r}| \over r_0}} { \text{H}(V_\varphi) \text{H}(1-q) \over (1-q)^{1 \over 2} } , 
~~ q(\mathbf{r},\mathbf{V}) \equiv {V_r^2 \over V_e(|\mathbf{r}|)^2} + {V_\theta^2  \over V_0^2} +  {V_\varphi^2 \over V_0^2}, </math> 
once we input the expression for the earlier potential <math>\Phi(r) </math> inside <math> r \le r_0 </math>, or even better the speed to "escape from r to the edge" <math> r_0 </math> of a uniform sphere <math> V_e(r)=V_0 \sqrt{1-{r^2 \over r_0^2}}. </math>  
Clearly the factor <math> {V_e(|\mathbf{r}|) \over \sqrt{2Q} } = \sqrt{\max[0,{1 \over 1-q}]}</math> in the DF (distribution function) is well-defined only if <math>Q \ge 0 \rightarrow q\le 1  </math>, which implies a narrow range on radius <math> 0 \le |\mathbf{r}|<r_0 </math> and excludes high velocity particles, e.g., <math> V_t > V_0 > V_e(r) </math>, from the distribution function (DF, i.e., phase space density).

In fact, the positivity carves the (<math> V_\varphi \ge 0 </math>) left-half of an ellipsoid in the <math>[V_r, V_\theta, V_\varphi]</math> velocity space ("velocity ellipsoid"), 
<math display="block">  
q(\mathbf{r},\mathbf{V}) \equiv {V_r^2 \over V_0^2 (1-r^2/r_0^2)}  +  \left({V_\theta^2  \over V_0^2} + {V_\varphi^2 \over V_0^2} \right) \equiv u_r^2 + u_\theta^2 + u_\varphi^2 \le 1,  
</math>
where <math> (u_r,u_\theta,u_\varphi)</math> is <math> (V_r, V_\theta,V_\varphi)</math> rescaled by the function <math> V_e(r)=V_0 \sqrt{1-r^2/r_0^2} </math> or <math> V_0 </math> respectively.

The velocity ellipsoid (in this case) has rotational symmetry around the r axis or <math> V_r </math> axis.  It is more squashed (in this case) away from the radial direction, hence more tangentially anisotropic because everywhere <math> V_e(r) < V_0 </math>, except at the origin, where the ellipsoid looks isotropic. Now we compute the moments of the phase space.
 
E.g., the resulting density (moment) is
<math display="block">\begin{align}
 \rho(r,\theta,\varphi) & = 
\int_{-V_e(r)}^{V_e(r)} dV_r \int_{-V_0}^{V_0} dV_\theta \int_{0}^{V_0} dV_\varphi
{C_0 \over V_0^3} \left({2Q \over V_0^2}\right)^{-1/2} \\
& = \int_{-1}^{1} \int_{-1}^{1} \int_0^1 { (V_e du_r) (V_0 du_\theta) (V_0 du_\varphi) C_0 
\over V_0^3 (1- r^2/r_0^2)^{1/2} (1 - q)^{1/2} }\left.\right|_{q=u_r^2 + u_\theta^2 + u_\varphi^2}\\
& = C_0  { \int_0^1 (1-u^2)^{-1/2} (2\pi u^2 du)} = \rho_0
\end{align}</math>
is indeed a spherical (angle-independent) and uniform (radius-independent) density inside the edge, where the normalisation constant <math> C_0 =2 \pi^{-2} \rho_0 </math>.

The streaming velocity is computed as the weighted mean of the velocity vector
<math display="block">\begin{align}
 \langle\mathbf{V}\rangle (\mathbf{x}) & \equiv 
{\int f d\mathbf{V}^3 \mathbf{V} \over \int f d\mathbf{V}^3 } \\
& = {1 \over \rho} \int f d\mathbf{V}^3  [V_r, V_\theta, V_\varphi] {C_0 V_0^2 (2Q)^{-1/2}} \\
& = \left[{ 
  \int_{-1}^{1} \!\!u_r...du_r, 
~~\int_{-1}^{1} \!\!u_\theta...du_\theta , 
~~\int_0^1 (2du_r) \int_{0}^{\sqrt{1-u_r^2}} \!\!(2du_\theta) \int_{0}^{\sqrt{1-u_r^2-u_\theta^2}} \!\!\!\!\!\!\!\!\!\!{du_\varphi  u_\varphi V_0 \over (1 - u_r^2 - u_\theta^2 - u_\varphi^2)^{1/2} } 
\over \int_{0}^1 (2\pi U dU) \int_{0}^{\sqrt{1-U^2}}  du_\varphi (1 - U^2 - u_\varphi^2)^{-1/2}  }\right] \\
& = \left[0,0,{4 V_0 \over 3\pi}\right] = \overline{\mathbf{V}(\mathbf{x})}, 
\end{align}</math>
where the global average (indicated by the overline bar) of flow implies uniform pattern of flat azimuthal rotation, but zero net streaming everywhere in the meridional <math> (r,\theta)</math> plane.

Incidentally, the angular momentum global average of this flat-rotation sphere is 
<math display="block"> \overline{\mathbf{r} \times \langle\mathbf{V}\rangle} = \int_0^{r_0}  {(\rho 4\pi r^2 dr) \over M_0} [0,0,r \langle V_\varphi\rangle] = [0, 0, {3 r_0\over 4} \overline{V_\varphi}]. </math>
Note global average of centre of mass does not change, so <math> \overline{\mathbf{V}_i(\mathbf{x})} =0 </math> due to global momentum conservation in each rectangular direction <math> i=x,y,z</math>, and this does not contradict the global non-zero rotation.

Likewise thanks to the symmetry of <math> f(r,\theta,\varphi,V_r,V_\theta,V_\varphi) = f(r,\theta,\pm \varphi, \pm V_r, \pm V_\theta,V_\varphi) </math>, we have
<math>  \langle\mathbf{(\pm V_r) V_\varphi}\rangle  =0 </math>, <math>~ \langle \mathbf{(\pm V_\theta) V_\varphi}\rangle  =0 </math>, <math>~ 
\langle\mathbf{(\pm V_r) V_\theta}\rangle  =0 </math> everywhere}.

Likewise the rms velocity in the rotation direction is computed by a weighted mean as follows, E.g.,
<math display="block">\begin{align}
\langle\mathbf{V}_\varphi^2\rangle(|\mathbf{x}|)
&\equiv {\int f d\mathbf{V}^3 V_\varphi^2 \over \rho(|\mathbf{r}|)} \\
& = {\int_0^1 (2du_r) \int_{0}^{\sqrt{1-u_r^2}} (2du_\theta) \int_{0}^{\sqrt{1-u_r^2-u_\theta^2}} du_\varphi { (u_\varphi V_0)^2 \over (1 - q)^{1/2} } \over \int_0^1 { (2\pi u^2 du) (1 - u^2 )^{-1/2} } } \\
& = 0.25V_0^2 = 0.5 \langle V_t^2 \rangle \\
& = {\!\!\int_0^1 (2du_r)\!\! \int_{0}^{\sqrt{1-u_r^2}} (2du_\varphi) \!\!\int_{0}^{\sqrt{1-u_r^2-u_\varphi^2}} du_\theta { (u_\theta V_0)^2 \over (1 - q)^{1/2} }
\over \int_0^1 { (2\pi u^2 du) (1 - u^2 )^{-1/2} } } \\
& =\langle\mathbf{V}_\theta^2\rangle(|\mathbf{x}|),  \\
\end{align}</math>

Here
<math> \langle V_t^2 \rangle = \langle V_\theta^2 +  V_\varphi^2\rangle =0.5V_0^2. </math>

Likewise
<math display="block"> \langle\mathbf{V}_r^2\rangle(\mathbf{x}) = {\!\!\int_0^1 (du_\varphi) \int_{0}^{\sqrt{1-u_\varphi^2}} \!\!(2du_\theta) \!\!\int_{0}^{\sqrt{1-u_\varphi^2-u_\theta^2}} \!\!\!{ (2du_r)(u_r V_e(r))^2 \over (1 - q)^{1/2} } \over \int_0^1 { (2\pi u^2 du) (1 - u^2 )^{-1/2} } } = \left({V_0 \over 2} \sqrt{1-{r^2 \over r_0^2}} \right)^2. </math>

So the pressure tensor or dispersion tensor is
<math display="block"> \begin{align}
\sigma^2_{ij}(\mathbf{r})= &  {P_{ij}(\mathbf{r}) \over \rho(\mathbf{r})} \\
  =&\langle\mathbf{V}_i\mathbf{V}_j\rangle- \langle\mathbf{V}_i\rangle\langle\mathbf{V}_j\rangle \\
=  & \begin{bmatrix}
    \left[1-({r \over r_0})^2\right]\left({V_0\over 2}\right)^2 & 0 & 0 \\
    0 & \left({V_0\over 2}\right)^2 & 0 \\
    0 & 0 & \left[1- ({8 \over 3 \pi})^2\right]\left({V_0\over 2}\right)^2 
    \end{bmatrix} 
\end{align}
</math>
with zero off-diagonal terms because of the symmetric velocity distribution.  
Note while there is no Dark Matter in producing the previous flat rotation curve, 
the price is shown by the reduction factor <math> {8 \over 3 \pi} = 0.8488</math> in the random velocity spread in the azimuthal direction. Among the diagonal dispersion tensor moments, <math> \sigma_\theta \equiv \sqrt{\sigma^2_{\theta\theta}} = 0.5V_0</math> is the biggest among the three at all radii, and <math> \sigma_\varphi \equiv \sqrt{\sigma^2_{\varphi\varphi}} \ge \sigma_r \equiv \sqrt{\sigma^2_{rr}} </math> only near the edge between <math> 0.8488 r_0 \le r \le r_0 </math>.

The larger tangential kinetic energy than that of radial motion seen in the diagonal dispersions is often phrased by an anisotropy parameter
<math display="block"> \beta(r) \equiv 
1 -  { 0.5\langle {\mathbf V_t}^2(|\mathbf{r}|) \rangle \over \langle{\mathbf V_r}^2\rangle(|\mathbf{r}|) } = 1 -  {\langle {\mathbf V_\theta}^2(|\mathbf{r}|) \rangle \over \langle{\mathbf V_r}^2\rangle(|\mathbf{r}|) } = - {r^2 \over r_0^2 - r^2} \le 0; </math>
a positive anisotropy would have meant that radial motion dominated, and a negative anisotropy means that tangential motion dominates (as in this uniform sphere).

=== A worked example of Virial Theorem ===

Twice kinetic energy per unit mass of the above uniform sphere is

<math display="block">\begin{align}
{2K \over M_0} & =  \overline{\langle V^2\rangle} \equiv  \langle \overline{V^2} \rangle \\
& = M_0^{-1} \int_0^{M_0} \langle V_\theta^2+V_\varphi^2+V_r^2 \rangle dM \\
& = M_0^{-1} \int_0^1 \left({V_0^2 \over 4}+{V_0^2 \over 4}+{ (1-x^2)V_0^2 \over 4}\right) d(x^3 M_0) = 0.6 V_0^2 , ~~ x \equiv {r \over r_0} = \left({M \over M_0}\right)^{1 \over 3},
\end{align}</math>
which balances the potential energy per unit mass of the uniform sphere, inside which <math>M \propto r^3 \propto x^3 </math>.

The average Virial per unit mass can be computed from averaging its local value <math> \mathbf{r} \cdot (-\mathbf{\nabla} \Phi)</math>, which yields
<math display="block">\begin{align}
 {W \over M_0} & = \overline{\mathbf{r} \cdot (-\mathbf{\nabla} \Phi)}\\
&=M_0^{-1} \int_0^{r_0}  \mathbf{r} \cdot {- G M \mathbf{r} \over |\mathbf{r}|^3} (\rho d\mathbf{r}^3) = -M_0^{-1} \int_0^{M_0} {G M \over |\mathbf{r}|} dM  \\
& = -M_0^{-1} \int_{0}^{M_0}  { G M ~ dM \over r_0~(M/M_0)^{1 \over 3} } = -{3 G M_0 \over 5 r_0} = -0.6 V_0^2,
\end{align} </math>
as required by the Virial Theorem.  For this self-gravitating sphere, we can also verify that the virial per unit mass equals the averages of half of the potential 
<math display="block"> \begin{align}
{E_\text{pot} \over M_0} & = \overline{\langle {\Phi \over 2}\rangle} \\
& = M_0^{-1} \int_{x>0}^{x<1}  {\Phi(r_0 x) \over 2} d (M_0 x^3) \\
& = {W \over M_0} = {-2K \over M_0}.\end{align}
</math>
Hence we have verified the validity of Virial Theorem for a uniform sphere under self-gravity, i.e., the gravity due to the mass density of the stars is also the gravity that stars move in self-consistently; no additional dark matter halo contributes to its potential, for example.

=== A worked example of Jeans Equation in a uniform sphere ===

Jeans Equation is a relation on how the pressure gradient of a system should be balancing the potential gradient for an equilibrium galaxy.  In our uniform sphere, the potential gradient or gravity is
<math display="block"> \nabla \Phi = {d \Phi \over dr} = {\Omega^2 r} \ge 0, ~~\Omega = {V_0 \over r_0}.  </math>

The radial pressure gradient
<math display="block"> -{d (\rho \sigma_r^2) \over \rho dr} = -{d \sigma_r^2 \over dr} - {\sigma_r^2 \over r} {d \log \rho \over d\log r} = {\Omega^2 r \over 2} + 0 \ge 0.  </math>

The reason for the discrepancy is partly due to centrifugal force
<math display="block"> {\bar{V}_\varphi^2 \over r} =  {(0.8488V_0)^2 \over r} > 0, </math>
and partly due to anisotropic pressure
<math display="block"> \begin{align}
{(\sigma_\theta^2 -\sigma_r^2) \over r} &=0.25 \Omega^2 r \ge 0\\
{(\sigma_\varphi^2 -\sigma_r^2) \over r} & = 0.25\Omega^2 r - {0.1801 V_0^2 \over r} = \pm , \end{align}</math>
so <math> 0.2643V_0 = \sigma_\varphi < \sigma_r =0.5V_0</math> at the very centre, 
but the two balance at radius <math> r=0.8488r_0</math>, and then 
reverse to <math> 0.2643 V_0 = \sigma_\varphi > \sigma_r =0 </math> at the very edge.

Now we can verify that
<math display="block">\begin{align}
{\partial \langle V_r \rangle \over \partial t} & = (-\sum_{i=x,y,z} V_i \partial_i \langle V_r\rangle ) - \cancel{ \langle V_r\rangle \over t_\text{fric} } - \nabla_r \Phi + \sum_{i=x,y,z}{-\partial_i (n \sigma^2_{ir}) \over n}
\\
&={\bar{V}_\theta^2 +\bar{V}_\varphi^2 -2\bar{V}_r^2 \over r} - 0 
-{\partial \Phi \over \partial r} 
+ \left[-{d (\rho \sigma_r^2) \over \rho dr} + {\sigma_\theta^2 + \sigma_\varphi^2 -2\sigma_r^2 \over r} \right]  \\
& =  {0 + (0.4244V_0)^2 - 2 \times 0 \over r} -(\Omega^2 r) + \\
& \left[ {\Omega^2 r \over 2} + 
{(0.5V_0)^2 + (0.2643V_0)^2- 2 \times 0.25 \Omega^2 (r_0^2-r^2) \over r} \right] \\
& = 0.
\end{align} </math>
Here the 1st line above is essentially the Jeans equation in the r-direction, which reduces to the 2nd line, the Jeans equation in an anisotropic (aka <math> \beta \ne 0 </math>) rotational (aka <math>\langle V_\varphi\rangle  \ne 0</math>) axisymmetric (<math> \partial_\varphi \Phi(\mathbf{x},t) =0</math> ) sphere (aka <math> \partial_\theta n(\mathbf{x},t) =0</math>) after much coordinate manipulations of the dispersion tensor; similar equation of motion can be obtained for the two tangential direction, e.g., <math> {\partial \langle V_\varphi\rangle \over \partial t} </math>, which are useful in modelling   ocean currents on the rotating earth surface or angular momentum transfer in accretion disks, where the frictional term <math> -{ \langle V_\varphi\rangle \over t_\text{fric} } </math> is important. 
The fact that the l.h.s. <math> {\partial V_r \over \partial t} =0 </math> means that 
the force is balanced on the r.h.s. for this uniform (aka <math> \nabla_\mathbf{x} m n(\mathbf{x},t) =0</math>) spherical model of a galaxy (cluster) to stay in a steady state (aka time-independent equilibrium <math> {\partial n(\mathbf{x},t) \over \partial t} = 0 </math> everywhere) statically (aka with zero flow <math> \langle \mathbf{V}(\mathbf{x},t) \rangle =0</math> everywhere).  Note systems like accretion disk can have a steady net radial inflow <math> \langle \mathbf{V}(\mathbf{x}) \rangle < 0</math> everywhere at all time.

=== A worked example of Jeans equation in a thick disk ===

Consider again the thick disk potential in the above example. 
If the density is that of a gas fluid, then the pressure would be zero at the boundary <math> z=\pm z_0 </math>.  To find the peak of the pressure, we note that 
<math display="block"> P(R,z) = \int^{z_0}_z \partial_z\Phi \rho(R) dz = \rho(R) [\Phi(R,z_0) - \Phi(R,z)].</math>

So the fluid temperature per unit mass, i.e., the 1-dimensional velocity dispersion squared would be 
<math display="block"> \sigma^2(R,z) = {P(R,z) \over \rho(R)},  ~~ |z| \le z_0  </math>

<math display="block"> \sigma^2= {G M_0 \over 2z_0} \log{Q(z) Q(-z) \over Q (z_0)Q(-z_0) },  ~~ Q(z) \equiv R_0 + z_0 + z +\sqrt{R^2+(R_0+z_0+z)^2}. </math>

Along the rotational z-axis, 
<math display="block"> \sigma^2(0,z) = {G M_0 \over 2z_0} \log {4(R_0+z_0+z) (R_0+z_0-z) \over 4R_0 (R_0+2z_0) } </math>
<math display="block"> \sigma(0,z) = \sqrt{G M_0 \over 2z_0} \sqrt{\log { (R_0+z_0)^2-z^2 \over (R_0+z_0)^2-z_0^2 } } , </math>
which is clearly the highest at the centre and zero at the boundaries <math> z=\pm z_0 </math>.  Both the pressure and the dispersion peak at the midplane <math> z=0 </math>.  In fact the hottest and densest point is the centre, where <math> P(0,0) = { M_0 \over 4 \pi R_0^2 z_0 } {-G M_0 \log [1- (1+R_0/z_0)^{-2}] \over 2 z_0} . </math>

=== A recap on worked examples on Jeans Eq., Virial and Phase space density ===

Having looking at the a few applications of Poisson Eq. and Phase space density and especially the Jeans equation, we can extract a general theme, again using the Spherical cow approach.

Jeans equation links gravity with pressure gradient, it is a generalisation of the Eq. of Motion for single particles.  While Jeans equation can be solved in disk systems, the most user-friendly version of the Jeans eq. is the spherical anisotropic version for a static <math> \langle{v_j}\rangle =0 </math> frictionless system <math> t_\text{fric} \rightarrow \infty</math>, hence the local velocity speed 
<math>  \sigma_j^2 (r) = \langle{v_j^2}\rangle(r)  - \underbrace{\langle{v_j}\rangle^2(r)}_{=0}  = { \int\limits_\infty\!\! dv_r dv_\theta dv_\varphi  ({v}_j-\overbrace{\langle{v}\rangle^p_j}^{=0})^2 f_p \over 
\int\limits_\infty\!\! dv_r dv_\theta dv_\varphi f_p},
</math> everywhere for each of the three directions <math> ~_j =~ _r, ~_\theta, ~_\varphi</math>.
One can project the phase space into these moments, which is easily if in a highly spherical system, which admits conservations of energy <math> E = </math> and angular momentum J.  The boundary of the system sets the integration range of the velocity bound in the system.

In summary, in the spherical Jeans eq., 
<math display="block"> \begin{align} {d \Phi \over d r} = & {GM(r) \over r^2}  \\
= & -{d(n \langle{v_r^2}\rangle ) \over n(r) dr}+ {\langle{v_\theta^2}\rangle  + \langle{ v_\phi^2}\rangle -2 \langle{v_r^2}\rangle \over r} ,  \\
= & -{d(n \langle{v_r^2}\rangle ) \over n(r) dr}, ~~\text{hydrostatic equilibrium if isotropic velocity } \\
= & { \langle v_t^2 \rangle \over r}, ~~\text{if purely centrifugal balancing of gravity with no radial motion}, \langle v_t^2 \rangle \equiv \langle{v_\theta^2}\rangle  + \langle{ v_\phi^2}\rangle  \end{align} </math>
which matches the expectation from the Virial theorem <math>\overline{r \partial_r \Phi} = \overline{v_\text{cir}^2} = \overline{G M \over r}= \overline{\langle v_t^2 \rangle} </math>, 
or in other words, the <math>\overline{\text{global average}}</math> kinetic energy of an equilibrium equals the average kinetic energy on circular orbits with purely transverse motion.