Past exam of the mathematics course of the University of Cambridge / 2017 / iii / Paper 320 (Galactic Astronomy and Dynamics)

The stellar relaxation time is the time on which accumulated discrete gravitational encounters change a typical star's velocity by an amount comparable to its original velocity. It concerns two-body relaxation, rather than the much faster orbital evolution in a smooth gravitational potential. A collisionless stellar system requires this time to greatly exceed the duration being studied.

In the straight-line impulse approximation, an encounter at relative speed $u$ has transverse acceleration $Gmb/(b^2+u^2{t}^2{)}^{3/2}$. ThusA small-angle gravitational encounter needs $\Delta v_{\perp }\ll u$. The transition to order-one deflections is therefore $b\sim Gm/u^2$, ignoring factors such as the equal-mass relative deflection. For the spherical estimate $u\sim v$,This is the lower cutoff of the weak-scattering estimate; closer encounters actually occur and require strong-scattering treatment.

The number of encounters in time $T$ with impact parameters in $[b,b+db]$ is $2\pi bdbnuT$. Independent random transverse directions make mean kicks cancel while their variances add. ConsequentlyEach logarithmic interval contributes equally, giving the Coulomb logarithm in stellar dynamics. Defining a deflection time by variance $v^2$ and setting $u=v$ gives $v^3/(8\pi G^2{m}^{2}n\log \Lambda )$. The paper instead uses an order-one normalization three quarters of this estimate:Both have the same physical scaling. The specified characteristic speed, the word “comparable”, the velocity-distribution average and strong-encounter cutoff do not fix that numerical coefficient uniquely. The crude equal-speed impulse calculation must not be claimed to determine $3/(32\pi )$ exactly.

For a self-gravitating, approximately virialized system, the virial theorem gives $v^2\sim GNm/R$. With $n=3N/(4\pi R^3)$, $\Lambda \sim N$ and the stellar crossing time $t_{\mathrm{dyn}}=R/v$, substitution into the paper's convention givesAn externally dominated gravitational potential or a different structural constant changes this substitution. For an N-body simulation lasting $T=K{t}_{\mathrm{dyn}}$, a tolerable fractional velocity-squared diffusion $\epsilon \ll 1$ requires $T/t_{\mathrm{rel}}\lesssim \epsilon$, henceThere is no unique smallest $N$ without $K$, the error tolerance and the force prescription. For scale, equality at $t_{\mathrm{rel}}/t_{\mathrm{dyn}}=100$ occurs near $N=7.1\times 10^3$; $N=10^5$ gives about $1.1\times 10^3$ crossing times. A calculation lasting 100 crossing times therefore needs substantially more than the first threshold to have negligible relaxation. Gravitational softening can raise the effective cutoff and reduce artificial scattering, but it also sets the spatial force resolution.

Put $u=\lambda v$ and let $n$ now mean number per area. During time $T$, the two strips of impact parameters between $b$ and $b+db$, on opposite sides of the trajectory, have combined area $2uTdb$. ThusThe same straight-line impulse approximation gives a kick $2Gm/(bu)$. Summing its squared size over these encounters givesassuming $b_{\max }\gg b_{\min }$. The relevant change is now the small random speed $u$, not the full circular speed $v$. Setting $D_v{t}_{\mathrm{rel}}\sim u^2$ yieldsThe razor-thin disk relaxation integral is dominated by nearby encounters, rather than equally by logarithmic intervals. Because kicks cease to be weak near the cutoff, its numerical coefficient is an estimate.

Here the weak-deflection cutoff must use the relative speed: $b_{\min }\sim Gm/u^2=Gm/({\lambda }^2{v}^2)$, which is larger than the spherical cutoff. Thus $t_{\mathrm{rel}}\sim u/(Gmn)$. Using $n\sim N/R^2$ and $v^2\sim GNm/R$ from the virial theorem givesFor example $n=N/(\pi R^2)$ gives the order-one coefficient $\pi /8$ in this last substitution; the paper's approximate equality suppresses such structural factors. The cancellation of $N$ is the central result.

An unsoftened point-particle system confined to an infinitesimally thin plane is therefore not collisionless over orbital times for $\lambda \sim 0.1$ in this encounter model. At fixed total mass, increasing $N$ decreases $m$ and the cutoff together; the enhanced importance of close planar encounters cancels the usual benefit of increasing particle number. This is a singular zero-thickness limit. A fixed nonzero gravitational softening length or physical thickness restores a different $N$ dependence, so the conclusion is not a prohibition on collisionless numerical models of disks.

For finite-thickness disk relaxation, assume a well-separated range $b_{\min }\ll z_0\ll R$ and approximately virialized self-gravity. The three factors have distinct origins. Replacing the characteristic encounter speed by $u=\lambda v$ shortens the diffusion time by ${\lambda }^3$: the variance-production rate has one inverse power of speed, and the target random velocity squared has two powers. Concentrating the same number of stars into volume of order $R^2{z}_0$, rather than $R^3$, raises the number volume density by $R/z_0$ and shortens the time by $z_0/R$. Finally, only separations below the thickness sample a three-dimensional encounter geometry. The logarithmic upper cutoff becomes $z_0$, and therefore the inverse-logarithm time acquires $\log (R/b_{\min })/\log (z_0/b_{\min })$. More distant encounters have the planar geometry and provide a convergent, non-logarithmic correction. A smaller Coulomb logarithm in stellar dynamics actually partly offsets the first two shortenings.

Thus, using the same cutoff for this comparison,For actual point particles the disk cutoff is $Gm/(\lambda v{)}^2$, whereas the spherical one is $Gm/v^2$; one must use the appropriate cutoff in each logarithm rather than silently identifying them. For a softened calculation it is instead of order the larger of the weak-deflection scale and the softening length. The formula is not valid when $z_0$ approaches that cutoff.

Write $h=z_0/R$ and $q=K/\epsilon$, with simulation duration $K{t}_{\mathrm{dyn}}$ and allowed fractional diffusion $\epsilon$. Combining the paper's spherical normalization and its comparison formula gives the leading criterionThe spherical comparison logarithm cancels when its numerical convention and a common comparison cutoff are used; order-one shape factors remain approximate. In an unsoftened disk with $v^2\sim GNm/R$, $z_0/b_{\min }\sim h{\lambda }^{2}N$. Hence the implicit requirement is $N/\log (h{\lambda }^{2}N)\gtrsim 8q/({\lambda }^{3}h)$, in the regime $h{\lambda }^{2}N\gg 1$. For $\lambda =h=0.1$, the stellar relaxation time is of order $10^{-4}$ of the spherical value apart from logarithms; $N\sim 10^8$ gives a time of order 100 crossing times in this estimate. Negligible relaxation over that duration requires an additional margin. A disk generally needs far more particles than a round galaxy for the same relaxation tolerance; no universal minimum follows without duration, thickness and softening. Collective spiral or bending responses are a separate source of evolution even in a collisionless stellar system.

Assume a time-independent spherical relative potential $\psi (r)$, with acceleration $\nabla \psi$. The relative energy and specific angular momentum areAlong an orbit, $\dot{E}=\nabla \psi \cdot \mathbf{v}-\mathbf{v}\cdot \nabla \psi =0$ and $\dot{\mathbf{\ell }}=\mathbf{r}\times \nabla \psi =0$, since a spherical gradient is radial. Thus both $E$ and $L$ are integrals of motion. Choose the escape level $\psi =0$ and take the bound galactic distribution function to vanish for $E\le 0$.

For an isotropic $f(E)$, integrate over a spherical coordinate system in velocity space:Differentiating gives ${\rho }^{\prime }(\psi )=2\pi \sqrt{2}{\int }_0^{\psi }f(E)/\sqrt{\psi -E}dE$. Convolve once more with the reciprocal-square-root kernel. Reversing the order of integration, under the usual integrability assumptions, usesTherefore ${\int }_0^{E_0}{\rho }^{\prime }(\psi )/\sqrt{E_0-\psi }d\psi =2\sqrt{2}{\pi }^2{\int }_0^{E_0}f(E)dE$. Differentiating proves Eddington inversion:Here ${\pi }^2$ is outside the square root, as in the original PDF; the TeX transcription wrongly places it inside. The assumptions include a locally integrable bound distribution and sufficient density regularity, with no extra unbound or boundary population. An inverted expression must additionally be nonnegative to be physical. Since the velocity measure and $f$ are invariant under all velocity rotations, $⟨v_r^2⟩=⟨v_{\theta }^2⟩=⟨v_{\varphi }^2⟩=⟨v^2⟩/3$ whenever these moments exist; mixed moments vanish.

For the half-anisotropic distribution function, write the polar velocity angle from the radial direction as $\eta$. Then $L=rv\sin \eta$, andThe apparent $1/L$ singularity is integrable because it cancels the factor $v\sin \eta$ in the velocity measure. Integrating $\eta \in [0,\pi ]$ and $\phi \in [0,2\pi ]$ givesThe density relation requires $r\rho \to 0$ at the escape boundary and $f_E\ge 0$ for a physical distribution. Consider a monotonically decreasing $\psi$. The largest allowed radius at fixed binding energy is the radial-orbit limit $r_E$ with $\psi (r_E)=E$; a nonzero-$L$ orbit normally has a smaller apocenter because of its centrifugal term. Changing from $\psi$ to $r$ givesThe numerator in the original PDF instead uses $\rho +d\rho /d\psi$. That is a genuine printed error: the chain rule for $d(r\rho )/d\psi$ requires $\rho +rd\rho /dr$. It is not legitimate to prove the printed expression as written. A concrete counterexample is the unit-parameter Hernquist model at $r=1$: ${\psi }^{\prime }=-1/4$, ${\rho }^{\prime }=-(5/2)\rho$ and $d\rho /d\psi =10\rho$. The printed quotient is $-44\rho$, whereas the correct $g$ is $6\rho >0$. The former would give a negative distribution.

The angular weight after cancellation is uniform in $\eta$, so the averages of ${\cos }^2\eta$ and ${\sin }^2\eta$ are both $1/2$. The azimuthal angle splits the tangential term equally. Hence $⟨v_r^2⟩=⟨v_{\theta }^2+v_{\varphi }^2⟩$ and each individual tangential moment is half the radial one. The velocity-anisotropy parameter is consequentlyThis is a radially biased constant-anisotropy distribution function, not an isotropic one or purely radial motion. It does not by itself prove dynamical stability; it preferentially weights small angular momentum while remaining integrable.

For the Hernquist model, use the Poisson equation for Newtonian gravity with the relative-potential sign, ${\nabla }^2\psi =-4\pi G\rho$. Since ${\psi }^{\prime }=-GM/(r+a{)}^2$,Its enclosed mass is $M{r}^2/(r+a{)}^2$, which tends to $M$. Expressing the augmented density in the relative potential gives $r\rho =a{\psi }^3/(2\pi G^3{M}^2)$. The preceding derivative therefore giveswith zero bound-model population outside the allowed energy domain. It is nonnegative and reproduces the density by direct integration. The Hernquist model has a central density cusp, so the formulas for the local $1/L$ distribution are interpreted at $r>0$. As a further check, the radial second moment obtained from the same integral is $\psi /4$ and each tangential second moment is $\psi /8$, consistent with $\beta =1/2$.

On a length scale $\ell$, pressure communicates support over a sound-crossing time $\ell /v_s$, whereas self-gravity acts over a time of order $(G{\rho }_0{)}^{-1/2}$. When the former is longer, gravity can amplify an overdensity before pressure smooths it. Thus Jeans instability occurs at scales of order $\ell \gtrsim v_s/\sqrt{G{\rho }_0}$; the precise numerical factor depends on whether scale means wavelength or inverse wavenumber.

Take a uniform static background, small perturbations, an inviscid barotropic fluid, and $\delta p=v_s^2\delta \rho$ with $v_s^2=(dp/d\rho {)}_0>0$. An infinite Newtonian homogeneous density cannot itself have zero background gravitational acceleration while satisfying the unmodified Poisson equation for Newtonian gravity. Specify the Jeans swindle: subtract the homogeneous field, or regard an external prescription as balancing it, and apply Poisson's equation only to perturbations. The linearized continuity equation and Euler equations for an inviscid fluid areTaking a time derivative of the first equation and substituting the second yields ${\partial }_t^2\delta \rho =v_s^2{\nabla }^2\delta \rho +4\pi G{\rho }_0\delta \rho$. A plane wave $e^{i(\mathbf{k}\cdot \mathbf{x}-\omega t)}$ therefore givesFor $0<k<k_J$ there is exponential growth at rate $\sqrt{4\pi G{\rho }_0-v_s^2{k}^2}$. At $k_J$ the density mode is neutral to this linear order, and above it there are acoustic oscillations. The critical wavelength is $2\pi /k_J=\sqrt{\pi }{v}_s/\sqrt{G{\rho }_0}$. Transverse velocity perturbations decouple as zero-frequency vortical modes in this idealization.

For the attractive Yukawa potential, assume a real screening parameter $\alpha \ge 0$. Its Green function obeys $({\nabla }^2-{\alpha }^2)(-Gm{e}^{-\alpha r}/r)=4\pi Gm{\delta }^{(3)}(\mathbf{r})$, including the delta source. Thus $\delta {\varphi }_{\mathbf{k}}=-4\pi G\delta {\rho }_{\mathbf{k}}/(k^2+{\alpha }^2)$. The same fluid elimination gives Jeans instability with Yukawa gravity:There exists a growing nonzero-wavenumber mode precisely whenEquality gives no exponentially growing nonzero mode. If screening exceeds this threshold, the finite-range gravitational response is too weak to overcome pressure at any wavelength. Otherwise all sufficiently long but nonconstant waves remain unstable, with a longer minimum unstable wavelength than in Newtonian gravity. For $\alpha >0$ the long-wave gravitational contribution also scales as $k^2$, so the growth rate tends to zero as $k\to 0$; screening does not produce a second positive lower cutoff in wavenumber. A homogeneous screened background can have the constant gravitational potential $-4\pi G{\rho }_0/{\alpha }^2$, with zero force, so no Newtonian background subtraction is needed for nonzero screening. The limits $k\to 0$ and $\alpha \to 0$ are nonuniform. An exact $k=0$ density change is not one of the growing finite-wavelength modes. In the pressureless edge case $v_s=0$, every finite $k>0$ is unstable for positive density and finite screening.

Use the same specified homogeneous Newtonian background and Jeans swindle, but describe stars by a mass-normalized galactic distribution function. The Collisionless Boltzmann equation is ${\partial }_{t}f+\mathbf{v}\cdot {\nabla }_{x}f-\nabla \varphi \cdot {\nabla }_{v}f=0$. Let $f=f_0(\mathbf{v})+f_1{e}^{i(\mathbf{k}\cdot \mathbf{x}-\omega t)}$ and similarly perturb the gravitational potential. To first order,Eliminating $f_1$ and the nonzero gravitational potential amplitude gives the collisionless Jeans dispersion relationFor growing modes $\mathrm{Im}\omega >0$ there is no real-velocity pole. Real and damped frequencies require a causal contour prescription or analytic continuation; one must not simply integrate through a pole without specifying it.

Take $b>0$ in the Cauchy velocity distribution. Its normalization is $\int f_0{d}^{3}v={\rho }_0$. Orient $\mathbf{k}$ along the $u=v_z$ axis and integrate over the perpendicular velocities:At the marginal mode $\omega =0$, $F^{\prime }(u)/u=-2{\rho }_{0}b/[\pi (u^2+b^2{)}^2]$, with its finite limiting value at zero. The provided Gamma function integral at power two givesConsequentlyTo check that this really separates growing modes, put $c=\omega /k$ in the upper half-plane. Residue integration, or integration by parts followed by the Cauchy resolvent, gives $\int F^{\prime }(u)/(u-c)du={\rho }_0/(c+ib{)}^2$. The upper-half-plane solution has $\omega =i(\sqrt{4\pi G{\rho }_0}-kb)$, which grows exactly for $0<k<k_J$. At the threshold it approaches the marginal mode continuously. Although $b$ plays the threshold role of a velocity scale, this distribution's second moment diverges: the radial integrand for $\int v^2{f}_0{d}^{3}v$ tends to a nonzero constant at large speed. It is therefore incorrect to identify $b^2$ with a finite velocity dispersion or to infer this result by substituting an rms speed into a fluid formula.

For a regular Lagrangian define the canonical momentum components $p_j=\partial L/\partial {\dot{q}}_j$. Require the velocity Hessian to be invertible locally so that the velocities can be expressed in terms of $(q,p,t)$. The Legendre transform in mechanics definesThe explicit time argument cannot generally be omitted when the given Lagrangian depends on time. Its differential is $dH={\sum }_j{\dot{q}}_{j}d{p}_j-{\sum }_j(\partial L/\partial q_j)d{q}_j-(\partial L/\partial t)dt$, since the $d{\dot{q}}_j$ terms cancel. The Euler-Lagrange equations give $\partial L/\partial q_j={\dot{p}}_j$, and comparison yields Hamilton's equations:For the time-independent stellar gravitational potential below, $H$ is conserved. Its sign convention differs from the binding energy in question 2: here $E=H=v^2/2-\psi$ is the ordinary specific orbital energy.

The prescribed nondegenerate prolate spheroidal coordinates require $0<\alpha <\beta$; equality collapses the coordinate family. Their physical velocity components are $v_j=h_j{\dot{q}}_j$, and $h_{\varphi }=R$. Per unit stellar mass the Lagrangian is $L=\frac{1}{2}{\sum }_j{h}_j^2{\dot{q}}_j^2+\psi$, so $p_j=h_j^2{\dot{q}}_j=h_j{v}_j$. Its Hamiltonian isThis uses acceleration $+\nabla \psi$, as specified, rather than importing the opposite gravitational potential convention.

For the Stäckel potential, use the Hamilton-Jacobi equation with principal function $S=W_{\lambda }(\lambda )+W_{\mu }(\mu )+l\varphi -Et$. Axisymmetry makes $p_{\varphi }=l=R{v}_{\varphi }$ constant. Write $d=\lambda -\mu$ and $A(\tau )=(\tau +\alpha )(\tau +\beta )$. Multiplying $H=E$ by $d$ and usingseparates it as $F(\lambda )-E\lambda =F(\mu )-E\mu$, where $F(\tau )=2A(\tau )(W_{\tau }^{\prime }{)}^2-(\alpha -\beta )l^2/[2(\tau +\alpha )]-G(\tau )$. Each side depends on its own coordinate only and therefore equals a constant, chosen as $-I$. Here $G(\tau )$ is the generating function of the gravitational potential, not the gravitational constant. HenceThis proves the conserved axisymmetric Stäckel third integral. Eliminating $E$ between the two expressions and using $p_j=h_j{v}_j$ givesAt regular generic phase-space points the three separation constants $E,l,I$ are independent; special circular or degenerate orbits can reduce their local independence. Thus the model is an integrable Hamiltonian system.

For the intended regular, phase-mixed steady model, Jeans theorem permits a galactic distribution function $f=F(E,l,I)$. At any fixed point, $E$ and $I$ depend on $v_{\lambda }$ and $v_{\mu }$ only through their squares, and $l$ depends only on $v_{\varphi }$. Thus $f$ is even separately in $v_{\lambda }$ and $v_{\mu }$, whether or not it is even in $v_{\varphi }$. Odd integration proves velocity-ellipsoid alignment in a Stäckel potential:The argument allows azimuthal streaming: it does not require zero mean $v_{\varphi }$. The means in the two meridional directions vanish, so the off-diagonal velocity covariance also vanishes. Second moments must exist for this conclusion.

There is a hypothesis behind the printed universal wording. On a generic nonresonant regular invariant torus, a continuous steady distribution is constant along a dense orbit and hence depends only on the torus integrals; this justifies the phase-mixed form. Bare stationarity alone is not sufficient in a degenerate resonant gravitational potential with additional integrals. For example $H=\frac{1}{2}(v^2+{\omega }^2{r}^2)$ is a Stäckel Hamiltonian: choose $G(\tau )=-\frac{1}{2}{\omega }^2[{\tau }^2+(\alpha +\beta )\tau ]$, up to a constant vacuum shift. The additional conserved quantity $K_{xy}=v_x{v}_y+{\omega }^{2}xy$ gives a positive normalizable steady tracer $f\propto \exp (-aH+ϵK_{xy})$ for $a>|ϵ|>0$. Its local $⟨v_x{v}_y⟩=ϵ/(a^2-{ϵ}^2)\ne 0$. At azimuth zero and a regular point where the $\lambda$ direction has a radial component, this gives nonzero $⟨v_{\lambda }{v}_{\varphi }⟩$. This tracer need not itself source the harmonic gravitational potential. It exhibits why the stronger, ordinary phase-mixed three-integral assumption is necessary instead of claiming that every imaginable steady distribution has the asserted alignment.

The geometry follows by multiplying the coordinate-defining equation and evaluating at $\tau =-\alpha ,-\beta$:Constant $\lambda$ surfaces are prolate spheroids, with equatorial semiaxis $\sqrt{\lambda +\alpha }$ and polar semiaxis $\sqrt{\lambda +\beta }$. Constant interior $\mu$ surfaces are two-sheeted hyperboloids. Their common foci lie at $z=\pm \sqrt{\beta -\alpha }$. The orbit volume in an axisymmetric Stäckel potential is determined byIts allowed coordinate intervals have turning boundaries or coordinate endpoints. A typical bound orbit with $l\ne 0$ oscillates between inner and outer spheroids, oscillates vertically between hyperboloidal limits, and circulates about the symmetry axis with $\dot{\varphi }=l/R^2$. The equatorial fold $\mu =-\beta$ joins the positive and negative $z$ branches; it is not necessarily a physical turning point in $z$.

For three incommensurate frequencies, the regular trajectory densely samples its projected allowed volume, with nonuniform occupation determined by its speeds. It does not follow one coordinate surface, and it need not close. Resonant orbits may be closed or lie in a smaller projected set; $l=0$ confines azimuthal motion to a meridional plane (with the usual axis coordinate qualification), and circular or equatorial limits are also lower-dimensional. Smooth Stäckel separation removes the generic chaotic behavior of nonintegrable potentials; it does not assert that every orbit fills a three-dimensional volume.