SUPPORT THE WORK

GetWiki

stellar dynamics

ARTICLE SUBJECTS
aesthetics  →
being  →
complexity  →
database  →
enterprise  →
ethics  →
fiction  →
history  →
internet  →
knowledge  →
language  →
licensing  →
linux  →
logic  →
method  →
news  →
perception  →
philosophy  →
policy  →
purpose  →
religion  →
science  →
sociology  →
software  →
truth  →
unix  →
wiki  →
ARTICLE TYPES
essay  →
feed  →
help  →
system  →
wiki  →
ARTICLE ORIGINS
critical  →
discussion  →
forked  →
imported  →
original  →
stellar dynamics
[ temporary import ]
please note:
- the content below is remote from Wikipedia
- it has been imported raw for GetWiki
Stellar dynamics is the branch of astrophysics which describes in a statistical way the collective motions of stars subject to their mutual gravity. The essential difference from celestial mechanics is that the number of body N gg 10. (File:Animation of Juno trajectory.gif|thumb|Slingshot of a test body in a two-body potential)(File:Hamiltonian flow classical.gif|thumb|N-particles in quasi-periodic motion in the phase space (x, mv) of an essentially static potential)Typical galaxies have upwards of millions of macroscopic gravitating bodies and countless number of neutrinos and perhaps other dark microscopic bodies. Also each star contributes more or less equally to the total gravitational field, whereas in celestial mechanics the pull of a massive body dominates any satellite orbits.BOOK, Encyclopedia of Astronomy and Astrophysics, Murdin, Paul, Nature Publishing Group, 2001, 978-0750304405, 1, Stellar Dynamics,

Connection with fluid dynamics

Stellar dynamics also has connections to the field of plasma physicscds.cern.ch/record/1053485/files/p37.pdf {{Bare URL PDF|date=March 2022}} The two fields underwent significant development during a similar time period in the early 20th century, and both borrow mathematical formalism originally developed in the field of fluid mechanics.In accretion disks and stellar surfaces, the dense plasma or gas particles collide very frequently, and collisions result in equipartition and perhaps viscosity under magnetic field. We see various sizes for accretion disks and stellar atmosphere, both made of enormous number of microscopic particle mass,
(L/V,M/N)
  • sim (10^{-8}text{pc}/500text{km/s},1M_odot/10^{55}=m_{p}) at stellar surfaces,
  • sim (10^{-4}{text{pc}}/10{text{km/s}},0.1M_{odot }/10^{54}sim m_{p}) around Sun-like stars or km-sized stellar black holes,
  • sim (10^{-1}{text{pc}}/100{text{km/s}},10M_{odot }/10^{56}sim m_{p}) around million solar mass black holes (about AU-sized) in centres of galaxies.
The system crossing time scale is long in stellar dynamics, where it is handy to note that
1000text{pc}/1text{km/s} = 1000 text{Myr} = text{HubbleTime}/14.
The long timescale means that, unlike gas particles in accretion disks, stars in galaxy disks very rarely see a collision in their stellar lifetime. However, galaxies collide occasionally in galaxy clusters, and stars have close encounters occasionally in star clusters.As a rule of thumb, the typical scales concerned (see the Upper Portion of P.C.Budassi’s Logarithmic Map of the Universe) are (L/V,M/N)
  • sim (mathrm{10pc/10km/s},1000M_odot/1000) for M13 Star Cluster,
  • sim (mathrm{100kpc/100km/s}, 10^{11}M_{odot }/10^{11}) for M31 Disk Galaxy,
  • sim (mathrm{10Mpc/1000km/s},10^{14}M_{odot }/10^{77}=m_nu) for neutrinos in the Bullet Clusters, which is a merging system of N = 1000 galaxies.

Connection with Kepler problem and 3-body problem

At a superficial level, all of stellar dynamics might be formulated as an N-body problemby Newton’s second law, where the equation of motion (EOM) for internal interactions of an isolated stellar system of N members can be written down as,m_ifrac{d^{2} mathbf{r_i}}{dt^{2}} = sum_{i=1 atop i ne j}^N frac{G m_i m_j left(mathbf{r}_j - mathbf{r}_iright)}{left| mathbf{r}_j - mathbf{r}_iright|^3}. Here in the N-body system, any individual member, m_i is influenced by the gravitational potentials of the remaining m_j members.In practice, except for in the highest performance computer simulations, it is not feasible to calculate rigorously the future of a large N system this way. Also this EOM gives very little intuition. Historically, the methods utilised in stellar dynamics originated from the fields of both classical mechanics and statistical mechanics. In essence, the fundamental problem of stellar dynamics is the N-body problem, where the N members refer to the members of a given stellar system. Given the large number of objects in a stellar system, stellar dynamics can address both the global, statistical properties of many orbits as well as the specific data on the positions and velocities of individual orbits.

Concept of a gravitational potential field

Stellar dynamics involves determining the gravitational potential of a substantial number of stars. The stars can be modeled as point masses whose orbits are determined by the combined interactions with each other. Typically, these point masses represent stars in a variety of clusters or galaxies, such as a Galaxy cluster, or a Globular cluster. Without getting a system’s gravitational potential by adding all of the point-mass potentials in the system at every second, stellar dynamicists develop potential models that can accurately model the system while remaining computationally inexpensive.BOOK, Galactic Dynamics, Binney, James, Tremaine, Scott, Princeton University Press, 2008, 978-0-691-13027-9, Princeton, 35, 63, 65, 698, The gravitational potential, Phi, of a system is related to the acceleration and the gravitational field, mathbf{g} by:frac {d^{2}mathbf {r_{i}} }{dt^{2}}}=mathbf {vec {g}} =-nabla _{mathbf {r_{i}} }Phi (mathbf {r_{i}} ),~~Phi (mathbf {r} _{i})=-sum _{k=1 atop kneq i}^{N}{frac {Gm_{k}}{left|mathbf {r} _{i}-mathbf {r} _{k}right|}, whereas the potential is related to a (smoothened) mass density, rho , via the Poisson’s equation in the integral form
Phi(mathbf {r}) = - int {G rho(mathbf{R}) d^3mathbf{R} over left|mathbf {r}-mathbf {R} right|}
or the more common differential form nabla^2Phi = 4pi G rho.

An example of the Poisson Equation and escape speed in a uniform sphere

Consider an analytically smooth spherical potential
begin{align}
Phi(r) & equiv left(-V_0^2right) + left[{r^2 -r_0^2 over 2r_0^2}, ~~ 1 -{r_0 over r} right]_{max} !!!! V_0^2 equiv Phi(r_0)-{V_e^2(r) over 2}, ~~Phi(r_0) = - V_0^2 , mathbf{g} &= -mathbf{nabla} Phi(r) = -Omega^2 r H(r_0 - r) - { G M_0 over r^2}H(r-r_0), ~~Omega={V_0 over r_0}, ~~M_0 = {V_0^2 r_0 over G},end{align}where V_e(r) takes the meaning of the speed to “escape to the edge” r_0, and sqrt{2}V_0 is the speed to “escape from the edge to infinity”. The gravity is like the restoring force of harmonic oscillator inside the sphere, and Keplerian outside as described by the Heaviside functions.We can fix the normalisation V_0 by computing the corresponding density using the spherical Poisson Equation
Grho = {d over 4 pi r^2 dr} {r^2 dPhi over dr} = { d (G M) over 4 pi r^2 dr} = {3 V_0^2 over 4 pi r_0^2}H(r_0-r),
where the enclosed mass
M(r) = {r^2 dPhi over G dr} = int_0^{r} dr int_0^{pi} (r dtheta) int_0^{2 pi} (r sintheta dvarphi) rho_0 H(r_0-r) = left. M_0 x^3right|_{x={r over r_0}}.
Hence the potential model corresponds to a uniform sphere of radius r_0 , total mass M_0 with
{V_0 over r_0} equiv sqrt{4pi G rho_0 over 3} = sqrt{G M_0 over r_0^3}.

Key concepts

While both the equations of motion and Poisson Equation can also take on non-spherical forms, depending on the coordinate system and the symmetry of the physical system, the essence is the same: The motions of stars in a galaxy or in a globular cluster are principally determined by the average distribution of the other, distant stars. The infrequent stellar encounters involve processes such as relaxation, mass segregation, tidal forces, and dynamical friction that influence the trajectories of the system’s members.JOURNAL, de Vita, Ruggero, Trenti, Michele, MacLeod, Morgan, 2019-06-01, Correlation between mass segregation and structural concentration in relaxed stellar clusters,doi.org/10.1093/mnras/stz815, Monthly Notices of the Royal Astronomical Society, 485, 4, 5752–5760, 10.1093/mnras/stz815, 1903.07619, 0035-8711,

Relativistic Approximations

There are three related approximations made in the Newtonian EOM and Poisson Equation above.

SR and GR

Firstly above equations neglect relativistic corrections, which are of order of
(v/c)^2 ll 10^{-4} as typical stellar 3-dimensional speed, v sim 3-3000 km/s, is much below the speed of light.

Eddington Limit

Secondly non-gravitational force is typically negligible in stellar systems. For example, in the vicinity of a typical star the ratio of radiation-to-gravity force on a hydrogen atom or ion,
Q^text{Eddington} = { {sigma_e over 4pi m_H c} {Lodot over r^2} over {G M_odot over r^2} } = {1 over 30,000},
hence radiation force is negligible in general, except perhaps around a luminous O-type star of mass 30M_odot, or around a black hole accreting gas at the Eddington limit so that its luminosity-to-mass ratio L_bullet / M_bullet is defined by Q^text{Eddington} =1 .

Loss cone

Thirdly a star can be swallowed if coming within a few Schwarzschild radii of the black hole. This radius of Loss is given by s le s_text{Loss} = frac{6 G M_bullet}{c^2} The loss cone can be visualised by considering infalling particles aiming to the black hole within a small solid angle (a cone in velocity). These particle with small theta ll 1 have small angular momentum per unit mass J equiv r v sintheta le J_text{loss} = frac{4G M_bullet}{c}. Their small angular momentum (due to ) does not make a high enough barrier near s_text{Loss} to force the particle to turn around.The effective potential
Phi_text{eff}(r) equiv E- {dot{r}^2 over 2} = {J^2 over 2r^2} + Phi(r) , is always positive infinity in Newtonian gravity. However, in GR, it
nosedives to minus infinity near frac{6 G M_bullet}{c^2} if J le frac{4G M_bullet}{c}. Sparing a rigorous GR treatment, one can verify this s_text{loss}, J_text{loss} by computing the last stable circular orbit, where the effective potential is at an inflection point Phi’’_text{eff}(s_text{loss})=Phi’_text{eff}(s_text{loss})=0 using an approximate classical potential of a Schwarzschild black hole
Phi(r) = - {(4G M_bullet/c)^2 over 2r^2} left[1+{3 (6 G M_bullet/c^2)^2 over 8 r^2 }right] - frac{G M_bullet}{r} left[1 - {(6G M_bullet/c^2)^2 over r^2}right].

Tidal disruption radius

A star can be tidally torn by a heavier black hole when coming within the so-called Hill’s radius of the black hole, inside which a star’s surface gravity yields to the tidal force from the black hole,WEB, Binney, James, Galaxy Dynamics,cds.cern.ch/record/1371794/files/9780691130279_TOC.pdf, Princeton University Press, 4 January 2022, i.e.,
(1-1.5) ge Q^text{tide}
equiv { G M_odot /R_odot^2 over [G M_bullet/s^2_text{Hill} - G M_bullet/(s_text{Hill}+R_odot)^2] }, ~~~s_text{Hill} rightarrow R_odot left({ (2-3) GM_bullet over GM_odot}right)^{1 over 3}, For typical black holes of M_bullet = (10^0-10^{8.5}) M_odot the destruction radius max[s_text{Hill}, s_text{Loss}] = 400R_odot maxleft[left({M_bullet over 3 times 10^7 M_odot}right)^{1/3}, {M_bullet over 3 times 10^7 M_odot}right] = (1-4000) R_odot ll 0.001 mathrm{pc}, where 0.001pc is the stellar spacing in the densest stellar systems (e.g., the nuclear star cluster in the Milky Way centre). Hence (main sequence) stars are generally too compact internally and too far apart spaced to be disrupted by even the strongest black hole tides in galaxy or cluster environment.

Radius of sphere of influence

A particle of mass m with a relative speed V will be deflected when entering the (much larger) cross section pi s^2_bullet of a black hole. This so-called sphere of influence is loosely defined by, up to a Q-like fudge factor sqrt{lnLambda} ,
1 sim sqrt{lnLambda} equiv frac{V^2/2}{G (M_bullet + m)/s_bullet},
hence for a Sun-like star we have,
s_bullet = {G (M_bullet +M_odot) sqrt{lnLambda} over V^2/2 } approx {M_bullet over M_odot} {V^2_odot over V^2} R_odot >[s_text{Hill}, s_text{Loss}]_{max} = (1-4000) R_odot,
i.e., stars will neither be tidally disrupted nor physically hit/swallowed in a typical encounter with the black hole thanks to the high surface escape speed V_odot =sqrt{2 G M_odot/R_odot} = 615mathrm{km/s} from any solar mass star, comparable to the internal speed between galaxies in the Bullet Cluster of galaxies, and greater than the typical internal speed V sim sqrt{2 G (N M_odot)/R} ll mathrm{300 km/s} inside all star clusters and in galaxies.

Connections between star loss cone and gravitational gas accretion physics

First consider a heavy black hole of mass M_bullet is moving through a dissipational gas of (rescaled) thermal sound speed text{Ï‚’} and density rho_text{gas} , then every gas particle of mass m will likely transfer its relative momentum m V_bullet to the BH when coming within a cross-section of radius s_bullet equiv {(G M_bullet+ G m) sqrt{lnLambda} over (V_bullet^2+text{Ï‚’}^2)/2}, In a time scale t_text{fric} that the black hole loses half of its streaming velocity, its mass may double by Bondi accretion, a process of capturing most of gas particles that enter its sphere of influence s_bullet , dissipate kinetic energy by gas collisions and fall in the black hole. The gas capture rate is
{M_bullet over t_text{Bondi}^{gas} }

sqrt{text{Ï‚’}^2 + V_bullet^2}(pi s_bullet^2) rho_text{gas}

4pi rho_text{gas} left[ {(G M_bullet)^2 over (text{Ï‚’}^2 + V_bullet^2)^{3 over 2} } right] lnLambda, ~~ text{Ï‚’} equiv sigma sqrt{1+ gamma^3 over 2 (9/8)^{2/3}} approx [text{Ï‚} , gamma sigma]_text{max},

where the polytropic index gamma is the sound speed in units of velocity dispersion squared, and the rescaled sound speed text{Ï‚’} allows us to match the Bondi spherical accretion rate, dot{M}_bullet approx pi rho_text{gas} text{Ï‚} left[ {(G M_bullet) over text{Ï‚}^2} right]^2 for the adiabatic gas gamma=5/3, compared to dot{M}_bullet approx 4pi rho_text{gas} text{Ï‚} left[ {(G M_bullet) over text{Ï‚}^2} right]^2 of the isothermal case gamma=1 .Coming back to star tidal disruption and star capture by a (moving) black hole, setting ln Lambda =1 , we could summarise the BH’s growth rate from gas and stars,
{M_bullet over t_text{Bondi}^{gas} } + {M_bullet over t_text{loss}^{*} } with,
dot{M}_bullet

sqrt{text{Ï‚’}^2 + V_bullet^2} m n (pi s_bullet^2, pi s_text{Hill}^2 , pi s_text{Loss}^2)_text{max}, ~~s_bullet approx {(G M_bullet+ G m) over (V_bullet^2+text{Ï‚’}^2)/2},

because the black hole consumes a fractional/most part of star/gas particles passing its sphere of influence.

Gravitational dynamical friction

Consider the case that a heavy black hole of mass M_bullet moves relative to a background of stars in random motion in a cluster of total mass (N M_odot) with a mean number density n sim (N-1)/(4pi R^3/3) within a typical size R .Intuition says that gravity causes the light bodies to accelerate and gain momentum and kinetic energy (see slingshot effect). By conservation of energy and momentum, we may conclude that the heavier body will be slowed by an amount to compensate. Since there is a loss of momentum and kinetic energy for the body under consideration, the effect is called dynamical friction.After certain time of relaxations the heavy black hole’s kinetic energy should be in equal partition with the less-massive background objects. The slow-down of the black hole can be described as
-{M_bullet dot{V}_bullet } = {M_bullet V_bullet over t_text{fric}^text{star} } ,
where t_text{fric}^text{star} is called a dynamical friction time.

Dynamical friction time vs Crossing time in a virialised system

Consider a Mach-1 BH, which travels initially at the sound speed text{Ï‚} = V_0 , hence its Bondi radius s_bullet satisfies
{GM_bullet sqrt{lnLambda} over s_bullet} = V_0^2 = text{Ï‚}^2 = { 0.4053 G M_odot (N-1) over R},
where the sound speed is text{Ï‚} = sqrt{ 4 G M_odot (N-1) over pi^2 R} with the prefactor {4 over pi^2} approx {4 over 10}=0.4 fixed by the fact that for a uniform spherical cluster of the mass density rho = n M_odot approx {M_odot (N-1) over 4.19 R^3} , half of a circular period is the time for “sound” to make a oneway crossing in its longest dimension, i.e., 2t_{text{Ï‚}} equiv 2t_{text{cross}} equiv {2R over text{Ï‚}} = pi sqrt{R^3 over G M_odot (N-1)} approx (0.4244 G rho)^{-1/2}. It is customary to call the “half-diameter” crossing time t_{text{cross}} the dynamical time scale.Assume the BH stops after traveling a length of l_text{fric} equiv text{Ï‚} t_text{fric} with its momentum M_bullet V_0=M_bullet text{Ï‚} deposited to {M_bullet over M_odot} stars in its path over l_text{fric}/(2R) crossings, then the number of stars deflected by the BH’s Bondi cross section per “diameter” crossing time is
N^text{defl} = { ({M_bullet over M_odot}) } {2R over l_text{fric}} = N {pi s_bullet^2 over pi R^2} = N left({M_bullet over 0.4053 M_odot N}right)^2 lnLambda. More generally, the Equation of Motion of the BH at a general velocity mathbf{V}_bullet in the potential Phi of a sea of stars can be written as

-{dover dt} (M_bullet V_bullet) - M_bullet nabla Phi equiv {(M_bullet V_bullet) over t_text{fric}} =
overbrace{ N pi s_bullet^2 over pi R^2}^{N^text{defl}} {(M_odot V_bullet) over 2t_text{Ï‚}} = { 8 lnLambda’ over N t_text{Ï‚}} M_bullet V_bullet, {pi^2 over 8} approx 1 and the Coulomb logarithm modifying factor {lnLambda’ over lnLambda} equiv left[{pi^2 over 8}right]^2 left[(1+ {V_bullet^2 over text{Ï‚’}^2})right]^{-2} (1+{M_odot over M_bullet}) le left[{text{Ï‚’} over V_bullet}right]^4 le 1 discounts friction on a supersonic moving BH with mass M_bullet ge M_odot . As a rule of thumb, it takes about a sound crossing t_text{Ï‚’} time to “sink” subsonic BHs, from the edge to the centre without overshooting, if they weigh more than 1/8th of the total cluster mass. Lighter and faster holes can stay afloat much longer.

More rigorous formulation of dynamical friction

The full Chandrasekhar dynamical friction formula for the change in velocity of the object involves integrating over the phase space density of the field of matter and is far from transparent.It reads as{M_bullet d (mathbf{V}_bullet) over dt} = -{M_bullet mathbf{V}_bullet over t_text{fric}^text{star} } = - {m mathbf{V}_bullet ~ n(mathbf{x}) dmathbf{x}^3 over dt} lnLambda_text{lag}, where
~~ n(mathbf{x}) dx^3 = dt V_{bullet} (pi s_bullet^2) n(mathbf{x}) = dt n(mathbf{x}) |V_{bullet}| pi left[{G(m+M_bullet) over |V_{bullet}|^2/2}right]^2
is the number of particles in an infinitesimal cylindrical volume of length |V_{bullet} dt| and a cross-section pi s_bullet^2 within the black hole’s sphere of influence.Like the “Couloumb logarithm” lnLambda factors in the contribution of distant background particles, here the factor ln(Lambda_text{lag}) also factors in the probability of finding a background slower-than-BH particle to contribute to the drag. The more particles are overtaken by the BH, the more particles drag the BH, and the greater is ln(Lambda_text{beaten}) . Also the bigger the system, the greater is lnLambda .A background of elementary (gas or dark) particles can also induce dynamical friction, which scales with the mass density of the surrounding medium, m~ n; the lower particle mass m is compensated by the higher number density n. The more massive the object, the more matter will be pulled into the wake.Summing up the gravitational drag of both collisional gas and collisionless stars, we have
M_bullet {d ( mathbf{V}_{bullet}) over M_bullet dt} =
- 4pi left[{GM_bullet over |V_{bullet}|}right]^2 mathbf{hat{V}}_{bullet} (rho_text{gas} lnLambda_text{lag}^{gas} + m n_text{*} lnLambda_text{lag}^{*}).~~Here the “lagging-behind” fraction for gas JOURNAL, Ostriker, Eva, Dynamical Friction in a Gaseous Medium, The Astrophysical Journal, 1999, 513, 1, 252, 10.1086/306858, astro-ph/9810324, 1999ApJ...513..252O, 16138105,ui.adsabs.harvard.edu/abs/1999ApJ...513..252O/abstract, and for stars are given by
begin{align}
lnLambda_text{lag}^{gas}(u) & = ln~ { left[{1+uover lambda}right]^{1 over 2} left[{|1-u|over lambda}right]^{H[u-lambda-1]-H[1-lambda-u] over 2} over exp{ [u+lambda,1]_min^2 - [u-lambda,1]_min^2 over 4 lambda} }, & approx ln left[ {sqrt{ (u^3 - 1)^2 + lambda^3 } + u^3 -1 over sqrt{1+lambda^3}-1 } right]^{1 over 3}, ~~ u equiv {|V_bullet| t over text{Ï‚’} t}, ~~ lambda equiv({s_bullet over text{Ï‚’}t}) {lnLambda_text{lag}^{*} over lnLambda} & equiv int_{0}^{|m V_{bullet}|} !!!! { (4pi p^2 dp) e^{-{p^2 over 2 (m sigma)^2}}over (sqrt{2pi} m sigma)^3 } left.right|_{p=m |v|} approx { |mathbf{V}_{bullet}|^3 over |mathbf{V}_{bullet}|^3 + 3.45 sigma^3 }, lnLambda &= int{dmathbf{x_1}^3 ~2 Heaviside[{n(mathbf{x_1}) over n(mathbf{x})} - 1 - {M_bullet over N M_odot} ] over (s_bullet^2 + |mathbf{x_1}-mathbf{x}|^2)^{3 over 2} } approx lnsqrt{1+left({0.123 N M_odot over M_bullet}right)^2 }, end{align} where we have further assumed that the BH starts to move from time t=0; the gas is isothermal with sound speed text{Ï‚} ; the background stars have of (mass) density m n(mathbf{x}) in a Maxwell distribution of momentum p=m v with a Gaussian distribution velocity spread sigma (called velocity dispersion, typically sigma le text{Ï‚} ).Interestingly, the G^2 (m+M_bullet) (m n(mathbf{x})) dependence suggests that dynamical friction is from the gravitational pull of by the wake, which is induced by the gravitational focusing of the massive body in its two-body encounters with background objects.We see the force is also proportional to the inverse square of the velocity at the high end, hence the fractional rate of energy loss drops rapidly at high velocities.Dynamical friction is, therefore, unimportant for objects that move relativistically, such as photons. This can be rationalized by realizing that the faster the object moves through the media, the less time there is for a wake to build up behind it. Friction tends to be the highest at the sound barrier, where lnLambda_text{lag}^{gas}left.right|_{u=1} =ln {text{Ï‚’}t over s_bullet } .

Gravitational encounters and relaxation

Stars in a stellar system will influence each other’s trajectories due to strong and weak gravitational encounters. An encounter between two stars is defined to be strong/weak if their mutual potential energy at the closest passage is comparable/minuscule to their initial kinetic energy. Strong encounters are rare, and they are typically only considered important in dense stellar systems, e.g., a passing star can be sling-shot out by binary stars in the core of a globular cluster.BOOK, Galaxies in the Universe, Sparke, Linda, Linda Sparke, Gallagher, John, Cambridge, 2007, 978-0521855938, New York, 131, This means that two stars need to come within a separation,
s_* = {G M_odot + G M_odot over V^2/2} = { 2 over 1.5}{G M_odot over text{Ï‚}^2} = {3.29 R over N-1},
where we used the Virial Theorem, “mutual potential energy balances twice kinetic energy on average”, i.e., “the pairwise potential energy per star balances with twice kinetic energy associated with the sound speed in three directions”,
1 sim Q^text{virial} equiv {overbrace{2K}^{(N M_odot) V^2} over |W|} = {N M_odottext{Ï‚}^2 + N M_odottext{Ï‚}^2 + N M_odottext{Ï‚}^2
over {N (N-1) over 2} {G M_odot^2 over R_{pair} } },where the factor N (N-1)/2 is the number of handshakes between a pair of stars without double-counting, the mean pair separation R_text{pair} ={pi^2 over 24} R approx 0.411234 R is only about 40% of the radius of the uniform sphere.Note also the similarity of the Q^text{virial} leftarrow rightarrow sqrt{lnLambda}.

Mean free path

The mean free path of strong encounters in a typically (N-1) = 4.19 n R^3 gg 100 stellar system is then
l_text{strong} = {1 over (pi s_*^2)n } approx {(N-1) over 8.117} R gg R ,
i.e., it takes about 0.123 N radius crossings for a typical star to come within a cross-section pi s_*^2 to be deflected from its path completely. Hence the mean free time of a strong encounter is much longer than the crossing time R/V .

Weak encounters

Weak encounters have a more profound effect on the evolution of a stellar system over the course of many passages. The effects of gravitational encounters can be studied with the concept of relaxation time. A simple example illustrating relaxation is two-body relaxation, where a star’s orbit is altered due to the gravitational interaction with another star.Initially, the subject star travels along an orbit with initial velocity, mathbf{v}, that is perpendicular to the impact parameter, the distance of closest approach, to the field star whose gravitational field will affect the original orbit. Using Newton’s laws, the change in the subject star’s velocity, delta mathbf{v}, is approximately equal to the acceleration at the impact parameter, multiplied by the time duration of the acceleration.The relaxation time can be thought as the time it takes for delta mathbf{v} to equal mathbf{v}, or the time it takes for the small deviations in velocity to equal the star’s initial velocity. The number of “half-diameter” crossings for an average star to relax in a stellar system of N objects is approximately {t_text{relax} over t_text{Ï‚}} = N^{text{relax}} backsimeq frac{0.123(N-1)}{ln (N-1)} gg 1from a more rigorous calculation than the above mean free time estimates for strong deflection.The answer makes sense because there is no relaxation for a single body or 2-body system. A better approximation of the ratio of timescales is left.frac{N’}{ln sqrt{1+ N’^2}}right|_{N’=0.123(N-2)}, hence the relaxation time for 3-body, 4-body, 5-body, 7-body, 10-body, ..., 42-body, 72-body, 140-body, 210-body, 550-body are about 16, 8, 6, 4, 3, ..., 3, 4, 6, 8, 16 crossings. There is no relaxation for an isolated binary, and the relaxation is the fastest for a 16-body system; it takes about 2.5 crossings for orbits to scatter each other. A system with N sim 10^2 - 10^{10} have much smoother potential, typically takes sim ln N’ approx (2-20) weak encounters to build a strong deflection to change orbital energy significantly.

Relation between friction and relaxation

Clearly that the dynamical friction of a black hole is much faster than the relaxation time by roughly a factor M_odot / M_bullet , but these two are very similar for a cluster of black holes,
N^text{fric} ={t_text{fric} over t_text{Ï‚}} rightarrow {t_text{relax} over t_text{Ï‚}} = N^text{relax} sim {(N-1) over 10-100}, ~ text{when}~ {M_bullet rightarrow m leftarrow M_odot}.
For a star cluster or galaxy cluster with, say, N=10^3, ~ R=mathrm{1 pc-10^5 pc}, ~ V=mathrm{1 km/s - 10^3 km/s }, we have t_{text{relax}} sim 100 t_text{Ï‚}approx 100 mathrm{Myr} -10 mathrm{Gyr} . Hence encounters of members in these stellar or galaxy clusters are significant during the typical 10 Gyr lifetime.On the other hand, typical galaxy with, say, N=10^6 - 10^{11} stars, would have a crossing time t_text{Ï‚} sim {1 mathrm{kpc} - 100 mathrm{kpc} over 1 mathrm{km/s} - 100 mathrm{km/s}} sim 100 mathrm{Myr} and their relaxation time is much longer than the age of the Universe. This justifies modelling galaxy potentials with mathematically smooth functions, neglecting two-body encounters throughout the lifetime of typical galaxies. And inside such a typical galaxy the dynamical friction and accretion on stellar black holes over a 10-Gyr Hubble time change the black hole’s velocity and mass by only an insignificant fraction
Delta sim {M_bullet over 0.1 N M_odot} {t over t_text{Ï‚}} le {M_bullet over 0.1% N M_odot}
if the black hole makes up less than 0.1% of the total galaxy mass N M_odot sim 10^{6-11}M_odot. Especially when M_bullet sim M_odot , we see that a typical star never experiences an encounter, hence stays on its orbit in a smooth galaxy potential.The dynamical friction or relaxation time identifies collisionless vs. collisional particle systems. Dynamics on timescales much less than the relaxation time is effectively collisionless because typical star will deviate from its initial orbit size by a tiny fraction t/t_{text{relax}} ll 1 . They are also identified as systems where subject stars interact with a smooth gravitational potential as opposed to the sum of point-mass potentials. The accumulated effects of two-body relaxation in a galaxy can lead to what is known as mass segregation, where more massive stars gather near the center of clusters, while the less massive ones are pushed towards the outer parts of the cluster.

A Spherical-Cow Summary of Continuity Eq. in Collisional and Collisionless Processes

Having gone through the details of the rather complex interactions of particles in a gravitational system, it is always helpful to zoom out and extract some generic theme, at an affordable price of rigour, so carry on with a lighter load.First important concept is “gravity balancing motion” near the perturber and for the background as a whole
text{Perturber Virial} approx {GM_bullet over s_bullet} approx V_text{cir}^2 approx langle V rangle^2 approx overline{langle V^2 rangle} approx sigma^2 approx left({R over t_text{Ï‚}}right)^2 approx c_text{Ï‚}^2 approx {G (N m) over R} approx text{Background Virial},
by consistently omitting all factors of unity 4pi , pi, ln text{Λ} etc for clarity, approximating the combined mass M_bullet + m approx M_bullet and being ambiguous whether the geometry of the system is a thin/thick gas/stellar disk or a (non)-uniform stellar/dark sphere with or without a boundary, and about the subtle distinctions among the kinetic energies from the local Circular rotation speed V_text{cir}, radial infall speed langle V rangle , globally isotropic or anisotropic random motion sigma in one or three directions, or the (non)-uniform isotropic Sound speed c_text{ς} to emphasize of the logic behind the order of magnitude of the friction time scale.Second we can recap very loosely summarise the various processes so far of collisional and collisionless gas/star or dark matter by Spherical cow style Continuity Equation on any generic quantity Q of the system:
{d Qover dt} approx {pm Q over ({l over c_text{Ï‚} }) }, ~text{Q being mass M, energy E, momentum (M V), Phase density f, size R, density} {N m over {4pi over 3} R^3}...,
where the pm sign is generally negative except for the (accreting) mass M, and the Mean free path l = c_text{Ï‚} t_text{fric} or the friction time t_text{fric} can be due to direct molecular viscosity from a physical collision Cross section, or due to gravitational scattering (bending/focusing/Sling shot) of particles; generally the influenced area is the greatest of the competing processes of Bondi accretion, Tidal disruption, and Loss cone capture,
s^2 approx maxleft[text{Bondi radius}~ s_bullet, text{Tidal radius}~s_text{Hill}, text{physical size}~ s_text{Loss cone}right]^2.
E.g., in case Q is the perturber’s mass Q = M_bullet , then we can estimate the Dynamical friction time via the (gas/star) Accretion rate
begin{align}dot{M}_bullet=& {M_bullet over t_text{fric} } approx int_0^{s^2} d(text{area}) ~(text{background mean flux}) approx s^2 (rho c_text{Ï‚})
approx & frac{text{Perturber influenced cross section}~(s^2)}{text{background system cross section}~(R^2) } times frac{text{background mass}~(N m)}{text{crossing time}~t_text{Ï‚} approx {R over c_text{Ï‚}} approx {1 over sqrt{G (N m) over R^3} sim sqrt{G rho} sim kappa } } approx & {G M_bullet over G t_text{Ï‚} } {G M_bullet over G (Nm) } approx (rho c_text{Ï‚}) left({G M_bullet over c_text{Ï‚}^2 }right)^2 ,~~text{if consider only gravitationally focusing,} approx & {M_bullet over N t_text{Ï‚} },~~text{if for a light perturber} M_bullet rightarrow m = M_odot rightarrow & 0, ~~text{if practically collisionless}~~N rightarrow infty,end{align} where we have applied the relations motion-balancing-gravity.In the limit the perturber is just 1 of the N background particle, M_bullet rightarrow m , this friction time is identified with the (gravitational) Relaxation time. Again all Coulomb logarithm etc are suppressed without changing the estimations from these qualitative equations.For the rest of Stellar dynamics, we will consistently work on precise calculations through primarily Worked Examples, by neglecting gravitational friction and relaxation of the perturber, working in the limit N rightarrow infty as approximated true in most galaxies on the 14Gyrs Hubble time scale, even though this is sometimes violated for some clusters of stars or clusters of galaxies.of the cluster.A concise 1-page summary of some main equations in Stellar dynamics and Accretion disc physics are shown here, where one attempts to be more rigorous on the qualitative equations above.(File:GAPbrief.pdf|thumb|Stellar dynamics Key concepts and equations)

Connections to statistical mechanics and plasma physics

The statistical nature of stellar dynamics originates from the application of the kinetic theory of gases to stellar systems by physicists such as James Jeans in the early 20th century. The Jeans equations, which describe the time evolution of a system of stars in a gravitational field, are analogous to Euler’s equations for an ideal fluid, and were derived from the collisionless Boltzmann equation. This was originally developed by Ludwig Boltzmann to describe the non-equilibrium behavior of a thermodynamic system. Similarly to statistical mechanics, stellar dynamics make use of distribution functions that encapsulate the information of a stellar system in a probabilistic manner. The single particle phase-space distribution function, f(mathbf{x},mathbf{v},t), is defined in a way such thatf(mathbf{x},mathbf{v},t) , dmathbf{x} , dmathbf{v} = dN where dN/N represents the probability of finding a given star with position mathbf{x} around a differential volume dmathbf{x} and velocity text{v} around a differential velocity space volume dmathbf{v}. The distribution function is normalized (sometimes) such that integrating it over all positions and velocities will equal N, the total number of bodies of the system. For collisional systems, Liouville’s theorem is applied to study the microstate of a stellar system, and is also commonly used to study the different statistical ensembles of statistical mechanics.

Convention and notation in case of a thermal distribution

In most of stellar dynamics literature, it is convenient to adopt the convention that the particle mass is unity in solar mass unit M_odot, hence a particle’s momentum and velocity are identical, i.e.,
mathbf{p} = m mathbf{v} = mathbf{v}, ~ m=1, ~ N_text{total} = M_text{total},


{dM over dx^3 dv^3} = f(mathbf{x},mathbf{v},t) = f(mathbf{x},mathbf{p},t) equiv {dN over dx^3 dp^3}
For example, the thermal velocity distribution of air molecules (of typically 15 times the proton mass per molecule) in a room of constant temperature T_0 sim mathrm{300K} would have a Maxwell distribution
f^text{Max}(x,y,z,m V_x,m V_y,m V_z) = {1 over (2pi hbar)^3} {1 over expleft({E(x,y,z,p_x,p_y,p_z) - mu over kT_0}right) + 1 }
f^text{Max} sim {1 over (2pi hbar/m)^3} e^{mu over kT_0 } e^ {-E over msigma_1^2},
where the energy per unit mass E/m = Phi(x,y,z) + (V_x^2 + V_y^2 + V_z^2)/2, where Phi(x,y,z) equiv g_0 z = 0and sigma_1 =sqrt{k T_0/m} sim mathrm{0.3km/s} is the width of the velocity Maxwell distribution, identical in each direction and everywhere in the room, and the normalisation constant e^{mu over kT_0} (assume the chemical potential mu sim (msigma_1^2) lnleft[n_0 left({sqrt{2pi}hbar over m sigma_1}right)^3right] ll 0 such that the Fermi-Dirac distribution reduces to a Maxwell velocity distribution) is fixed by the constant gas number density n_0 = n(x,y,0) at the floor level, where n(x,y,0) = !! int_{-infty}^infty m dV_x !! int_{-infty}^infty m dV_y !! int_{-infty}^infty m dV_z f(x,y,0,mV_x,mV_y,mV_z)
n approx {(2pi)^{3/2} (msigma_1)^3 over (2pi hbar)^3} e^{mu over m sigma_1^2}.

The CBE

In plasma physics, the collisionless Boltzmann equation is referred to as the Vlasov equation, which is used to study the time evolution of a plasma’s distribution function.The Boltzmann equation is often written more generally with the Liouville operator {mathcal{L}} as{mathcal{L}} f(t,mathbf{x},mathbf{p}) = {f^text{Max}_text{fit} - f(t,mathbf{x},mathbf{p}) over t_text{relax}},
{mathcal{L}} equiv frac{partial}{partial t} + frac{mathbf{p}}{m} cdot nabla + mathbf{F}cdotfrac{partial}{partial mathbf{p}},.
where mathbf{F} equiv mathbf{dot{p}} =-m nabla Phi is the gravitational force and f^text{Max}_text{fit} is the Maxwell (equipartition) distribution (to fit the same density, same mean and rms velocity as f(t,mathbf{x},mathbf{p})). The equation means the non-Gaussianity will decay on a (relaxation) time scale of t_text{relax} , and the system will ultimately relaxes to a Maxwell (equipartition) distribution.Whereas Jeans applied the collisionless Boltzmann equation, along with Poisson’s equation, to a system of stars interacting via the long range force of gravity, Anatoly Vlasov applied Boltzmann’s equation with Maxwell’s equations to a system of particles interacting via the Coulomb Force.JOURNAL, Henon, M, June 21, 1982, Vlasov Equation?, Astronomy and Astrophysics, 114, 1, 211–212, 1982A&A...114..211H, Both approaches separate themselves from the kinetic theory of gases by introducing long-range forces to study the long term evolution of a many particle system. In addition to the Vlasov equation, the concept of Landau damping in plasmas was applied to gravitational systems by Donald Lynden-Bell to describe the effects of damping in spherical stellar systems.JOURNAL, Lynden-Bell, Donald, 1962, The stability and vibrations of a gas of stars, Monthly Notices of the Royal Astronomical Society, 124, 4, 279–296, 10.1093/mnras/124.4.279, 1962MNRAS.124..279L, A nice property of f(t,x,v) is that many other dynamical quantities can be formed by its moments, e.g., the total mass, local density, pressure, and mean velocity. Applying the collisionless Boltzmann equation, these moments are then related by various forms of continuity equations, of which most notable are the Jeans equations and Virial theorem.

Probability-weighted moments and hydrostatic equilibrium

Jeans computed the weighted velocity of the Boltzmann Equation after integrating over velocity space
{1 over rho_p } int! left{mathbf{v}_p {d [f_p m_p]over dt} - langle{mathbf{v}}rangle_p {d [f_p m_p]over dt}right} d^3mathbf{v} = 0, and obtain the Momentum (Jeans) Eqs. of a ^population (e.g., gas, stars, dark matter):
begin{align}overbrace{ left({partial over partial t}+sum_{j=1}^{3} langle{v_j^p}rangle {partial over partial x_j}right) langle{v_i^p}rangle}^{dot{langle{v}rangle}_i^p} &underbrace{=}_{EoM} overbrace{-partial Phi(t,mathbf{x})over partial x_i}^{g_isim O(-GM/R^2)} ~~underbrace{-}^text{pressure}_text{balance}~~sum_{j=1}^{3} {partial over rho^p partial x_j}overbrace{[underbrace{rho^p(t,mathbf{x})}_{int_infty!!!!m_p f_p d^3mathbf{v}} underbrace{sigma_{ji}^p(t,mathbf{x})}_{O(c_s^2)}]}^{intlimits_infty!! dmathbf{v}^3 (mathbf{v}_j-langle{v}rangle^p_j) (mathbf{v}_i-langle{v}rangle^p_i)m_pf_p }- {underbrace{langle{v_i^p}rangle overbrace{[dot{m}_p/m_p]}^{1/t|^text{fric}_{text{visc}~m_p=M_text{gas}}}}_text{snow.plough}}, The general version of Jeans equation, involving (3 x 3) velocity moments is cumbersome. It only becomes useful or solvable if we could drop some of these moments, epecially drop the off-diagonal cross terms for systems of high symmetry, and also drop net rotation or net inflow speed everywhere.The isotropic version is also called Hydrostatic equilibrium equation where balancing pressure gradient with gravity; the isotropic version works for axisymmetric disks as well, after replacing the derivative dr with vertical coordinate dz. It means that we could measure the gravity (of dark matter) by observing the gradients of the velocity dispersion and the number density of stars.

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.JOURNAL, Einstein, Albert, 2002, A Simple Application of the Newtonian Law of Gravitation to Star Clusters,alberteinstein.info/vufind1/images/einstein/5-166.tr.pdf, The Collected Papers of Albert Einstein, 7, 230–233, Princeton University Press, JOURNAL, Zwicky, Fritz, 2009, Republication of: The redshift of extragalactic nebulae, General Relativity and Gravitation, 41, 1, 207–224, 10.1007/s10714-008-0707-4, 2009GReGr..41..207Z, 119979381, 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.BOOK, Astrophysics for Physicists, Choudhuri, Arnab Rai, Cambridge University Press, 2010, 978-0-521-81553-6, New York, 213–214, 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 mergers. 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 begin{align} Phi(R,z) & ={G M_0 over 2z_0}
left[2sinh^{-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} where z_0, R_0 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 M_0 because
Phi(R,z)
rightarrow {G M_0 over 2z_0} (2 Q_{-} - Q_{-} -Q_{+}) = -{G M_0 over R} , when we take the large radii limit
R rightarrow infty, ~|z| ge z_0, , so that Q = Q_{-}=Q_{+}-{2z_0 over R} = {|z| + (R_0 - z_0) over R} rightarrow 0.
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 M_0 , and that of a uniform-Needle mass distribution:
Phi_{KM}(R,z) = -{G M_0 over sqrt{ R^2 + (|z|+R_0)^2}}, ~~ z_0=0,
Phi_{PT}(R,z) = -{G M_0 over sqrt{R^2+z^2}} , ~~ z_0=R_0=0,
Phi_{UN}^{R_0=0}(R, z) = {G M_0 over 2z_0} left[2sinh^{-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].

A worked example of gravity vector field in a thick disk

First consider the vertical gravity at the boundary,
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,
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, partial_{|z|} (2 Q) - partial_{|z|} Q_{-} = partial_{|z|} left(Q_{+} - frac{2z_0}{R}right) = {1 over R} is continuous. Apply Gauss’s theorem by integrating the vertical force over the entire disk upper and lower boundaries, we have
2 int_0^infty (2 pi R dR) |g_z(R,z_0)| = 4 pi G M_0,
confirming that M_0 takes the meaning of the total disk mass.The vertical gravity drops with -g_z rightarrow G M_0 z (1+R_0/z_0)/R^3 at large radii, which is enhanced over the vertical gravity of a point mass G M_0 z/R^3 due to the self-gravity of the thick disk.

Density of a thick disk from Poisson Equation

Insert in the cylindrical Poisson eq.
rho(R,z) ={partial_z partial_z Phi over 4 pi G} + {partial_R (Rpartial_R Phi) over 4 pi G R} = { M_0 R_0/z_0 over 4pi (R^2+R_0^2)^{3/2}} H(z_0-|z|),
which drops with radius, and is zero beyond |z| > z_0 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 2 z_0 , we find the surface density and the total mass as
Sigma(R) = (2 z_0)rho(R,0), ~~ M_0 = int_0^infty (2pi R dR) Sigma(R).
This confirms that the absence of extra razor thin discs at the boundaries. In the limit, z_0 rightarrow 0 , this thick disc potential reduces to that of a razor-thin Kuzmin disk, for which we can verify {|g_z (R,0+)| over 2pi G} rightarrow Sigma(R) rightarrow {M_0 R_0 over 2pi (R^2+R_0^2)^{3/2}} .

Oscillation frequencies in a thick disk

To find the vertical and radial oscillation frequencies, we do a Taylor expansion of potential around the midplane.
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
and we find the circular speed V_text{cir} and the vertical and radial epicycle frequencies to be given by
(R omega)^2 equiv V^2_text{cir} = left[{(1+R_0/z_0) G M_0over sqrt{R^2+(R_0+z_0)^2} } - {(R_0/z_0) G M_0 over sqrt{R^2+R_0^2}} right],
nu^2 = {G M_0 (R_0/z_0 + 1) over (R^2+(R_0+z_0)^2)^{3/2}},
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}}.
Interestingly the rotation curve V_text{cir} is solid-body-like near the centre R ll R_0 , and is Keplerian far away.At large radii three frequencies satisfy
left.left[omega, nu, kappa, sqrt{4 pi G rho}right]right|_{Rto infty} to [1,1+R_0/z_0,1, R_0/z_0]^{1over 2} sqrt{G M_0over R^3} .
E.g., in the case that R to infty and R_0 / z_0 = 3, the oscillations omega: nu: kappa = 1: 2 : 1 forms a resonance.In the case that R_0 =0 , the density is zero everywhere except uniform needle between |z| le z_0 along the z-axis.If we further require z_0=0, then we recover a well-known property for closed ellipse orbits in point mass potential,
omega: nu: kappa = 1: 1 : 1 .

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
f(mathbf{x},mathbf{v},t) = {dN over dx^3 dv^3} le {6 over (2pi hbar/m)^3}, ~~~
where the Fermi-Dirac statistics says there are at most 6 flavours of neutrinos within a volume dx^3 and a velocity volume dv^3 = (dp/m)^3 = [(2pihbar/dx)/m]^3,.Let’s approximate the distribution is at maximum, i.e.,
f(x, y, z, V_x, V_y, V_z) = {6 over (2pi hbar/m)^3} q^{alpha over 2}, ~~ 0 le q(E)={Phi_{max}- E over V_0^2/2} le 1,
where 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} such that E_{min}, E_{max} , 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 berho(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), which reduces torho(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} Bleft(1+{alpha over 2}, {3 over 2}right) over (2pi hbar/m)^3} Take the simple case alpha to 0 , and estimate the density at the centre r=0 with an escape speed V_0 , we have
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]}.
Clearly eV-scale neutrinos with m_{eV} sim 0.1-1 is too light to make up the 100–10000 over-density in galaxies with escape velocity V_{200} equiv V/(mathrm{200km/s}) sim 0.1-3.4 , while neutrinos in clusters with V sim mathrm{2000 km/s} could make up 100-1000 times cosmic background density.By the way the freeze-out cosmic neutrinos in your room have a non-thermal random momentum sim {(mathrm{2.7 K}) k over c} sim (1~mathrm{eV}/c^2) (mathrm{70 km/s}) , 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 rho_0 and potential Phi(r)
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, ~~ Omegaequiv 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}where V_e(r) = V_0 sqrt{1-{r^2 over r_0^2}} =sqrt{2Phi(r_0)-2Phi(r)} is the speed to escape to the edge r_0.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 Omega with 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}Loosely speaking our goal is to put stars on a weighted distribution of orbits with various energies fleft(I_x(x,dot{x}), I_y(y,dot{y}), I_z(z,dot{z}right) = DF(mathbf{r},mathbf{V}), 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 f(mathbf{x},mathbf{v}) 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 (r,theta,phi) and its velocity components (V_r,V_theta,V_phi) is
f(r,theta,varphi,V_r,V_theta,V_varphi) = {C_0 over V_0^3} sqrt{V_0^2 over 2Q},
where C_0 = 2pi^{-2} rho_0 is a normalisation constant, which has the dimension of (mass) density. And we define a (positive enthalpy-like dimension text{km}^2/text{s}^2 ) Quantity
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|}, 0right]_max .
Clearly anti-clockwise rotating stars with J_z le 0,~~ Q=0 are excluded.It is easy to see in spherical coordinates that
J^2 = r^2 V_t^2 = r^2 (V_theta^2+V_varphi^2),


J_z = V_varphi r sintheta,


E = {V_r^2+V_t^2 over 2} + Phi(r), ~ V_t equiv sqrt{V_theta^2+V_varphi^2}
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
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^2right)}, 0 right]_max, which implies |V_r| le V_e(r), and |V_theta|, V_varphi between zero and V_0 .
To verify the above E, ~J_z being constants of motion in our spherical potential, we note
dE/dt = {partial Eover partial t} + mathbf{v} {partial E over partial mathbf{x}} + (mathbf{-nabla Phi}) {partial E over partial mathbf{v}}


dE/dt

{partial Phiover partial t} + mathbf{v} {partial Phi over partial mathbf{x}} + (mathbf{-nabla Phi}) mathbf{v} {partial Phiover partial t} 0 for any “steady state” potential.

dJ_z/dt

{partial J_zover partial t} + {partial J_z over partial mathbf{x}} cdot mathbf{v} - (mathbf{nabla Phi}) cdot {partial J_z over partial mathbf{v}}, which reduces to 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) {yover R}{partial Phi(R,z) over partial R}right]

0 around the z-axis of any axisymmetric potential, where Rsqrt{x^2+y^2} .

Likewise the x and y components of the angular momentum are also conserved for a spherical potential. Hence dJ/dt =0 .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
dE[mathbf{x},mathbf{v}]/dt = dJ[mathbf{x},mathbf{v}]/dt= dJ_z[mathbf{x},mathbf{v}]/dt =0 .
Hence using the chain rule, we have
{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 ,
i.e., {d over dt} f= f’(Q) {d Q[mathbf{x},mathbf{v}]over dt} =0 , so that CBE is satisfied, i.e., our
f(mathbf{x},mathbf{v}) = f(E[mathbf{x},mathbf{v}],J[mathbf{x},mathbf{v}],J_z[mathbf{x},mathbf{v}])
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,
f(|mathbf{r}|,V_r,V_theta,V_varphi) =
{C_0 over V_0^3} left.{text{H}(1-x) over left(1-x^2right)^{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}, once we input the expression for the earlier potential Phi(r) inside r le r_0 , or even better the speed to “escape from r to the edge” r_0 of a uniform sphere V_e(r)=V_0 sqrt{1-{r^2 over r_0^2}}. Clearly the factor {V_e(|mathbf{r}|) over sqrt{2Q} } = sqrt{max[0,{1 over 1-q}]} in the DF (distribution function) is well-defined only if Q ge 0 rightarrow qle 1 , which implies a narrow range on radius 0 le |mathbf{r}| V_0 > V_e(r) , from the distribution function (DF, i.e., phase space density).In fact, the positivity carves the ( V_varphi ge 0 ) left-half of an ellipsoid in the [V_r, V_theta, V_varphi] velocity space (“velocity ellipsoid“),
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, where (u_r,u_theta,u_varphi) is (V_r, V_theta,V_varphi) rescaled by the function V_e(r)=V_0 sqrt{1-r^2/r_0^2} or V_0 respectively.The velocity ellipsoid (in this case) has rotational symmetry around the r axis or V_r axis. It is more squashed (in this case) away from the radial direction, hence more tangentially anisotropic because everywhere V_e(r) < V_0 , except at the origin, where the ellipsoid looks isotropic. Now we compute the moments of the phase space.
E.g., the resulting density (moment) isbegin{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} (2pi u^2 du)} = rho_0end{align}is indeed a spherical (angle-independent) and uniform (radius-independent) density inside the edge, where the normalisation constant C_0 =2 pi^{-2} rho_0 .The streaming velocity is computed as the weighted mean of the velocity vectorbegin{align}
langlemathbf{V}rangle (mathbf{x}) & equiv
{int f dmathbf{V}^3 mathbf{V} over int f dmathbf{V}^3 } & = {1 over rho} int f dmathbf{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 (2pi 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 3pi}right] = overline{mathbf{V}(mathbf{x})}, end{align}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 (r,theta) plane.Incidentally, the angular momentum global average of this flat-rotation sphere is
overline{mathbf{r} times langlemathbf{V}rangle} = int_0^{r_0} {(rho 4pi r^2 dr) over M_0} [0,0,r langle V_varphirangle] = [0, 0, {3 r_0over 4} overline{V_varphi}].
Note global average of centre of mass does not change, so overline{mathbf{V}_i(mathbf{x})} =0 due to global momentum conservation in each rectangular direction i=x,y,z, and this does not contradict the global non-zero rotation.Likewise thanks to the symmetry of f(r,theta,varphi,V_r,V_theta,V_varphi) = f(r,theta,pm varphi, pm V_r, pm V_theta,V_varphi) , we have
langlemathbf{(pm V_r) V_varphi}rangle =0 , ~ langle mathbf{(pm V_theta) V_varphi}rangle =0 , ~
langlemathbf{(pm V_r) V_theta}rangle =0 everywhere}.Likewise the rms velocity in the rotation direction is computed by a weighted mean as follows, E.g.,begin{align}langlemathbf{V}_varphi^2rangle(|mathbf{x}|)&equiv {int f dmathbf{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 { (2pi 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 { (2pi u^2 du) (1 - u^2 )^{-1/2} } } & =langlemathbf{V}_theta^2rangle(|mathbf{x}|), end{align}Here
langle V_t^2 rangle = langle V_theta^2 + V_varphi^2rangle =0.5V_0^2.
Likewise
langlemathbf{V}_r^2rangle(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 { (2pi u^2 du) (1 - u^2 )^{-1/2} } } = left({V_0 over 2} sqrt{1-{r^2 over r_0^2}} right)^2.
So the pressure tensor or dispersion tensor is
begin{align}
sigma^2_{ij}(mathbf{r})= & {P_{ij}(mathbf{r}) over rho(mathbf{r})}
=&langlemathbf{V}_imathbf{V}_jrangle- langlemathbf{V}_iranglelanglemathbf{V}_jrangle

& begin{bmatrix}

left[1-({r over r_0})^2right]left({V_0over 2}right)^2 & 0 & 0
0 & left({V_0over 2}right)^2 & 0
0 & 0 & left[1- ({8 over 3 pi})^2right]left({V_0over 2}right)^2
end{bmatrix}
end{align}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 {8 over 3 pi} = 0.8488 in the random velocity spread in the azimuthal direction. Among the diagonal dispersion tensor moments, sigma_theta equiv sqrt{sigma^2_{thetatheta}} = 0.5V_0 is the biggest among the three at all radii, and sigma_varphi equiv sqrt{sigma^2_{varphivarphi}} ge sigma_r equiv sqrt{sigma^2_{rr}} only near the edge between 0.8488 r_0 le r le r_0 .The larger tangential kinetic energy than that of radial motion seen in the diagonal dispersions is often phrased by an anisotropy parameter
beta(r) equiv
1 - { 0.5langle {mathbf V_t}^2(|mathbf{r}|) rangle over langle{mathbf V_r}^2rangle(|mathbf{r}|) } = 1 - {langle {mathbf V_theta}^2(|mathbf{r}|) rangle over langle{mathbf V_r}^2rangle(|mathbf{r}|) } = - {r^2 over r_0^2 - r^2} le 0; 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 isbegin{align}{2K over M_0} & = overline{langle V^2rangle} 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}which balances the potential energy per unit mass of the uniform sphere, inside which M propto r^3 propto x^3 .The average Virial per unit mass can be computed from averaging its local value mathbf{r} cdot (-mathbf{nabla} Phi), which yieldsbegin{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 dmathbf{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} 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
begin{align}
{E_text{pot} over M_0} & = overline{langle {Phi over 2}rangle} & = M_0^{-1} int_{x>0}^{x 0, and partly due to anisotropic pressure
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.25Omega^2 r - {0.1801 V_0^2 over r} = pm , end{align}so 0.2643V_0 = sigma_varphi < sigma_r =0.5V_0 at the very centre, but the two balance at radius r=0.8488r_0, and then reverse to 0.2643 V_0 = sigma_varphi > sigma_r =0 at the very edge.Now we can verify thatbegin{align}{partial langle V_r rangle over partial t} & = (-sum_{i=x,y,z} V_i partial_i langle V_rrangle ) - cancel{ langle V_rrangle 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 -2bar{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 -2sigma_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} 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 beta ne 0 ) rotational (aka langle V_varphirangle ne 0) axisymmetric ( partial_varphi Phi(mathbf{x},t) =0 ) sphere (aka partial_theta n(mathbf{x},t) =0) after much coordinate manipulations of the dispersion tensor; similar equation of motion can be obtained for the two tangential direction, e.g., {partial langle V_varphirangle over partial t} , which are useful in modelling ocean currents on the rotating earth surface or angular momentum transfer in accretion disks, where the frictional term -{ langle V_varphirangle over t_text{fric} } is important. The fact that the l.h.s. {partial V_r over partial t} =0 means that the force is balanced on the r.h.s. for this uniform (aka nabla_mathbf{x} m n(mathbf{x},t) =0) spherical model of a galaxy (cluster) to stay in a steady state (aka time-independent equilibrium {partial n(mathbf{x},t) over partial t} = 0 everywhere) statically (aka with zero flow langle mathbf{V}(mathbf{x},t) rangle =0 everywhere). Note systems like accretion disk can have a steady net radial inflow langle mathbf{V}(mathbf{x}) rangle < 0 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 z=pm z_0 . To find the peak of the pressure, we note that
P(R,z) = int^{z_0}_z partial_zPhi rho(R) dz = rho(R) [Phi(R,z_0) - Phi(R,z)].
So the fluid temperature per unit mass, i.e., the 1-dimensional velocity dispersion squared would be
sigma^2(R,z) = {P(R,z) over rho(R)}, ~~ |z| le z_0


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}.
Along the rotational z-axis,
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) }
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 } } ,
which is clearly the highest at the centre and zero at the boundaries z=pm z_0 . Both the pressure and the dispersion peak at the midplane z=0 . In fact the hottest and densest point is the centre, where 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} .

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 langle{v_j}rangle =0 frictionless system t_text{fric} rightarrow infty, hence the local velocity speed
sigma_j^2 (r) = langle{v_j^2}rangle(r) - underbrace{langle{v_j}rangle^2(r)}_{=0} = { intlimits_infty!! dv_r dv_theta dv_varphi ({v}_j-overbrace{langle{v}rangle^p_j}^{=0})^2 f_p over
intlimits_infty!! dv_r dv_theta dv_varphi f_p},
everywhere for each of the three directions ~_j =~ _r, ~_theta, ~_varphi.
One can project the phase space into these moments, which is easily if in a highly spherical system, which admits conservations of energy E = 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.,
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}

which matches the expectation from the Virial theorem overline{r partial_r Phi} = overline{v_text{cir}^2} = overline{G M over r}= overline{langle v_t^2 rangle} , or in other words, the overline{text{global average}} kinetic energy of an equilibrium equals the average kinetic energy on circular orbits with purely transverse motion.

See also

Further reading

References

{{reflist}}{{Astronomy subfields}}{{Star}}

- content above as imported from Wikipedia
- "stellar dynamics" does not exist on GetWiki (yet)
- time: 8:55am EDT - Wed, May 22 2024
[ this remote article is provided by Wikipedia ]
LATEST EDITS [ see all ]
GETWIKI 21 MAY 2024
GETWIKI 09 JUL 2019
Eastern Philosophy
History of Philosophy
GETWIKI 09 MAY 2016
GETWIKI 18 OCT 2015
M.R.M. Parrott
Biographies
GETWIKI 20 AUG 2014
CONNECT