NBPhO 2025

Overview

The problem strings together one piece of textbook physics — vertical oscillations in a self-gravitating slab — with one observational trick that turned a phase-space photograph into a measurement of the local dark-matter density. Three ideas do the heavy lifting:

1. Gauss’s law in 1D. For a horizontally infinite, mirror-symmetric mass distribution $\rho(z)=\rho(-z)$, the only field component is vertical, and Poisson’s equation collapses to $d^2\Phi/dz^2 = 4\pi G\rho(z)$. The vertical force $a_z=-d\Phi/dz$ is therefore set entirely by the enclosed surface mass between the midplane and $z$.
2. Anharmonic oscillator with amplitude-dependent period. Because $\rho$ falls with $|z|$, the effective spring weakens at larger amplitudes, so a high-flying star takes longer to come back than a low-flying one. The problem pins this down with a “matching parabola” prescription: at amplitude $z_m$, replace $\Phi$ by $kz^2$ with $k=\Phi(z_m)/z_m^2$, giving $$T(z_m)\;=\;\frac{2\pi}{\sqrt{2k}}\;=\;\pi z_m\sqrt{\tfrac{2}{\Phi(z_m)}}.$$
3. Phase mixing and the spiral. Stars perturbed at the same time but with different periods drift apart in phase. After time $t$, the phase difference between two stars is $(\omega_1-\omega_2)\,t$, which the spiral reads off directly: each crossing of an axis advances the spiral phase by $\pi/2$. This both gives $\Phi(z)$ at the turning points (where energy equals $\Phi$) and times the perturbation.

A small parameter that tells you why the phase-spiral method works at all is the fractional spread of periods,

$$\frac{T_\text{outer}-T_\text{inner}}{T_\text{outer}}\;\sim\;\frac{1}{2},$$

across the observed range of amplitudes. The factor of two between fastest and slowest oscillators is exactly what produces a tightly wound, easily readable spiral on a $\sim$ Gyr timescale.

A useful unit conversion throughout:

$$1\,\text{kpc}/(\text{km}/\text{s}) \;=\; 0.978\,\text{Gyr}\;\approx\; 978\,\text{Myr}.$$

Combined with $G = 4.30\times 10^{-3}\,\text{pc}\,M_\odot^{-1}(\text{km/s})^2$, all numerical work below stays in $\{$kpc, km/s, $M_\odot/\text{pc}^3\}$ with a single conversion at the end.

---

Part (i) — Acceleration in a uniform slab

For a horizontally infinite mass sheet, the field has only a $z$-component, $\vec g = g_z(z)\hat z$, with $g_z(0)=0$ by reflection symmetry. Apply Gauss’s law to a “pillbox” of cross-section $A$ extending from $-z$ to $+z$:

$$\oint \vec g\cdot d\vec A \;=\; -4\pi G M_\text{enc}.$$

The end caps contribute $g_z(z)\,A - g_z(-z)\,A = 2g_z(z)\,A$ (using $g_z(-z)=-g_z(z)$ by symmetry), and $M_\text{enc}=\rho_0\,(2zA)$. Hence

$$2g_z(z)\,A \;=\; -4\pi G\,\rho_0\,(2zA),$$ $$\boxed{\,a_z \;=\; g_z(z)\;=\;-\,4\pi G\rho_0\,z\,}.$$

The acceleration grows linearly with $z$ and points back toward the midplane — exactly the form of a Hooke’s-law restoring force, so the motion is harmonic.

Educational remark

The same result follows from $\nabla^2\Phi=4\pi G\rho$: in 1D with $\rho=\rho_0$ const, $\Phi(z)=2\pi G\rho_0 z^2$ (taking $\Phi(0)=0$). The factor $4\pi$, rather than $2\pi$ as one might guess by analogy with a single sheet, comes from including matter on both sides of the test point.

---

Part (ii) — Oscillation period

The equation of motion $\ddot z = -4\pi G\rho_0\,z$ is that of a SHO with

$$\omega^2 \;=\; 4\pi G\rho_0,$$

so a star released from rest at $z=a$ executes $z(t)=a\cos\omega t$ with

$$\boxed{\,T\;=\;\frac{2\pi}{\omega}\;=\;\sqrt{\frac{\pi}{G\rho_0}}\,}.$$

Notice that $T$ is independent of the amplitude $a$ — a feature special to a strictly harmonic potential. As soon as $\rho$ varies with $z$, the period acquires an amplitude dependence, and that is what produces the spiral in part (iii).

Numerical value

For the canonical local density $\rho_0 \approx 0.1\,M_\odot/\text{pc}^3$,

$$T \;=\; \sqrt{\frac{\pi}{4.30\times 10^{-3}\cdot 0.1}}\;\text{pc}/(\text{km/s}) \;\approx\; 85.5\;\text{pc}/(\text{km/s}) \;\approx\; 84\,\text{Myr}.$$

This $\sim 80\,\text{Myr}$ vertical “epicyclic” period is the standard timescale for a star bobbing through the disc in the Solar neighbourhood.

---

Part (iii) — Reading $\Phi(z)$ from the phase spiral

Method

Each star moves on a closed orbit in the $(z,v_z)$ plane satisfying conservation of energy,

$$E \;=\; \tfrac{1}{2}v_z^2 + \Phi(z).$$

The plotted spiral is a continuous curve on which energy varies monotonically from inner to outer end. Two kinds of axis-crossings give us direct probes of $E$:

• At a crossing of $v_z=0$ at $z=z_m$: that star is at its turning point, so $E = \Phi(z_m)$. We learn $\Phi$ if we know $E$.
• At a crossing of $z=0$ at $v_z=v_m$: $E = \tfrac{1}{2}v_m^2$ directly, no $\Phi$ needed.

The trick is that the two kinds of crossings interleave along the spiral. So we can read $E$ off the $z=0$ crossings (where it is unambiguous), then linearly interpolate $E$ along the spiral to a neighbouring $v_z=0$ crossing, and finally identify that interpolated value with $\Phi(z_m)$.

Concretely, label the crossings consecutively along the spiral as $n=0,1,2,\ldots$ (alternating type). For each $v_z=0$ crossing at $z=z_m^{(n)}$, the two adjacent $z=0$ crossings have $|v_z|=v_-$ and $v_+$. Linear interpolation of $E$ in the crossing index gives

$$\boxed{\,\Phi(z_m^{(n)}) \;\approx\; \tfrac{1}{2}\!\left[\tfrac{1}{2}v_-^2 + \tfrac{1}{2}v_+^2\right] \;=\; \tfrac{1}{4}\!\left(v_-^2 + v_+^2\right).\,}$$

Reading the spiral

Tracing from the inner end of the spiral outward, the axis crossings appear in alternating order. From the figure I read (signs alternating between consecutive crossings of the same type):

order	type	value	$E=\tfrac12 v_m^2$ $[(\text{km/s})^2]$
0	$v_z=0$	$\lvert z_m\rvert=0.10\,\text{kpc}$	—
1	$z=0$	$\lvert v_m\rvert=8\,\text{km/s}$	$32$
2	$v_z=0$	$\lvert z_m\rvert=0.18\,\text{kpc}$	—
3	$z=0$	$\lvert v_m\rvert=18\,\text{km/s}$	$162$
4	$v_z=0$	$\lvert z_m\rvert=0.32\,\text{kpc}$	—
5	$z=0$	$\lvert v_m\rvert=22\,\text{km/s}$	$242$
6	$v_z=0$	$\lvert z_m\rvert=0.50\,\text{kpc}$	—
7	$z=0$	$\lvert v_m\rvert=30\,\text{km/s}$	$450$
8	$v_z=0$	$\lvert z_m\rvert=0.70\,\text{kpc}$	—
9	$z=0$	$\lvert v_m\rvert=36\,\text{km/s}$	$648$
10	$v_z=0$	$\lvert z_m\rvert=1.00\,\text{kpc}$	—
11	$z=0$	$\lvert v_m\rvert=44\,\text{km/s}$	$968$
12	$v_z=0$	$\lvert z_m\rvert=1.40\,\text{kpc}$	—
13	$z=0$	$\lvert v_m\rvert=55\,\text{km/s}$	$1513$

(Readings are quantized to the nearest gridline; their precision sets the dominant uncertainty in everything that follows.)

Computing $\Phi(z_m)$

Apply $\Phi(z_m^{(n)}) = \tfrac14(v_-^2 + v_+^2)$ to each $v_z=0$ crossing whose two neighbours are known:

$$\begin{aligned} \Phi(0.18\,\text{kpc}) &\;=\; \tfrac{1}{4}(8^2 + 18^2) \;=\; 97 \;(\text{km/s})^2,\\ \Phi(0.32\,\text{kpc}) &\;=\; \tfrac{1}{4}(18^2 + 22^2) \;=\; 202,\\ \Phi(0.50\,\text{kpc}) &\;=\; \tfrac{1}{4}(22^2 + 30^2) \;=\; 346,\\ \Phi(0.70\,\text{kpc}) &\;=\; \tfrac{1}{4}(30^2 + 36^2) \;=\; 549,\\ \Phi(1.00\,\text{kpc}) &\;=\; \tfrac{1}{4}(36^2 + 44^2) \;=\; 808,\\ \Phi(1.40\,\text{kpc}) &\;=\; \tfrac{1}{4}(44^2 + 55^2) \;=\; 1240. \end{aligned}$$

(The innermost crossing at $|z_m|=0.10\,\text{kpc}$ has no inner neighbour — the spiral simply ends — so it is dropped as unreliable.)

Plot

A schematic of $\Phi(z)$ from these six points:

 Φ [(km/s)²]
 1500 ┤                                            ●   z=1.40
 1250 ┤
 1000 ┤                                  ●            z=1.00
  750 ┤                          ●                    z=0.70
  500 ┤                  ●                            z=0.50
  250 ┤          ●                                    z=0.32
      ┤   ●                                           z=0.18
    0 ┼──┼─────┼─────┼─────┼─────┼─────┼─────┼──── z [kpc]
      0  0.2   0.4   0.6   0.8   1.0   1.2   1.4

The curve is approximately parabolic for small $z$ (constant-density region) and bends toward a roughly linear regime at large $z$ (where the slab acts like a thin sheet, $\Phi\propto z$). Quantifying this transition is exactly the next two parts.

Educational remark

A natural worry is that linear interpolation of $E$ in the crossing index is arbitrary. It is justified by the matching-parabola picture: in the harmonic limit, the spiral is a uniformly-pitched logarithmic-like winding in which each axis crossing is equally spaced in phase angle (one per quarter-period). Energies of nearby orbits differ by an amount that, in the same limit, is also linear in phase index. Both linearities ride on the same approximation, so they are consistent.

A consistency check: $\Phi(z)/z^2$ is the local “matching $k$”, and on the assumption $\rho(z)$ decreases with $|z|$ it should decrease monotonically with $z_m$. Computing:

$$\frac{\Phi}{z^2}\;:\;\;3000,\;1972,\;1384,\;1120,\;808,\;633\;\;[(\text{km/s})^2/\text{kpc}^2],$$

indeed monotonically decreasing from $z=0.18$ to $z=1.40$. This is one of the cleanest confirmations from the data alone that the disc density genuinely falls off with height.

---

Part (iv) — Density near the midplane

Within the constant-density region $|z|\leq 0.3\,\text{kpc}$, the potential is parabolic,

$$\Phi(z) \;=\; 2\pi G\rho_0\,z^2,$$

so any single $\Phi(z)$ measurement inside that region directly yields $\rho_0$. Using the most reliable inner point at $z=0.18\,\text{kpc} = 180\,\text{pc}$:

$$\rho_0 \;=\; \frac{\Phi(0.18)}{2\pi G\,(180\,\text{pc})^2} \;=\; \frac{97}{2\pi\cdot 4.30\times 10^{-3}\cdot (180)^2}\,M_\odot/\text{pc}^3.$$

Computing the denominator: $2\pi\cdot 4.30\times 10^{-3}\cdot 32400 = 875.5$, so

$$\boxed{\,\rho_0 \;\approx\; \frac{97}{875.5}\,M_\odot/\text{pc}^3 \;\approx\; 0.11\,M_\odot/\text{pc}^3\,}.$$

Cross-check

A useful sanity check — almost circular but worthwhile — is that part (ii)‘s formula gives the corresponding period:

$$T \;=\; \sqrt{\pi/(G\rho_0)} \;\approx\; 81\,\text{Myr},$$

which matches the period one would extract from the innermost spiral orbits ($T(0.18)\approx 79\,\text{Myr}$, computed in part (vi) below). The two routes agree, as they must.

Educational remark

The number $\rho_0\sim 0.1\,M_\odot/\text{pc}^3$ is the so-called Oort limit, the dynamically-determined total mass density at the Solar position. Star counts give the visible density at $\sim 0.08\,M_\odot/\text{pc}^3$; the small ($\sim 20\%$) excess is what dynamicists have argued for a century is “missing matter”. Phase-space methods like this one make that excess visible in a single picture.

---

Part (v) — Local dark-matter density

Beyond $z\geq 0.7\,\text{kpc}$, the problem instructs us to take $\rho(z)\approx\rho_\text{DM}=\text{const}$. Poisson’s equation then integrates twice:

$$\Phi(z) \;-\; \Phi(0.7) \;=\; g_0\,(z-0.7)\;+\;2\pi G\rho_\text{DM}\,(z-0.7)^2,\qquad z\geq 0.7,$$

where $g_0\equiv g(0.7)$ is the field at the boundary (sourced by all the visible matter inside $|z|<0.7$). Two unknowns, $g_0$ and $\rho_\text{DM}$, demand two equations from the data; use the two outer measurements at $z=1.00$ and $z=1.40\,\text{kpc}$:

$$\begin{aligned} z=1.00\,\text{kpc}\;:\quad &808 - 549 \;=\; 259 \;=\; 0.30\,g_0 \;+\; 0.09\,q,\\ z=1.40\,\text{kpc}\;:\quad &1240 - 549 \;=\; 691 \;=\; 0.70\,g_0 \;+\; 0.49\,q, \end{aligned}$$

where $q\equiv 2\pi G\rho_\text{DM}$ in $(\text{km/s})^2/\text{kpc}^2$.

Eliminating $g_0$: from the first, $g_0 = (259 - 0.09\,q)/0.30$. Substituting,

$$691 \;=\; 0.70\cdot\frac{259-0.09q}{0.30} + 0.49q \;=\; 604.3 - 0.21q + 0.49q \;=\; 604.3 + 0.28q,$$

so $0.28q = 86.7$ and

$$q \;=\; 309\,(\text{km/s})^2/\text{kpc}^2.$$

Convert $q$ to $\rho_\text{DM}$:

$$\rho_\text{DM} \;=\; \frac{q}{2\pi G} \;=\; \frac{309\times 10^{-6}\,(\text{km/s})^2/\text{pc}^2}{2\pi\cdot 4.30\times 10^{-3}\,\text{pc}\cdot(\text{km/s})^2/M_\odot} \;=\; \frac{309\times 10^{-6}}{0.0270}\,M_\odot/\text{pc}^3,$$ $$\boxed{\,\rho_\text{DM} \;\approx\; 0.011\,M_\odot/\text{pc}^3\,}.$$

Sanity checks

• $g_0 = (259 - 0.09\cdot 309)/0.30 \approx 770\,(\text{km/s})^2/\text{kpc}$. Comparing with $g_0\approx 2\pi G\Sigma_\text{vis}$ for a thin disc: $\Sigma_\text{vis}\approx g_0/(2\pi G)\approx 28\,M_\odot/\text{pc}^2$, in the right ballpark for the Solar-neighbourhood visible surface density inside $\pm 0.7\,\text{kpc}$.
• $\rho_\text{DM}/\rho_0 \approx 0.011/0.11 = 10\%$. Of the total matter near the midplane, roughly $90\%$ is baryons (gas + stars) and $10\%$ is dark matter, consistent with independent disc-dynamics estimates.
• $0.011\,M_\odot/\text{pc}^3 \approx 0.4\,\text{GeV/cm}^3$, the standard “particle-physics” value used in direct-detection experiments.

The fact that all three numbers fall out at once from a single phase-space photograph is the achievement of this method.

Educational remark

The simplification ”$\rho_\text{DM}$ constant for $|z|\geq 0.7\,\text{kpc}$” exploits a clean separation of scales: dark-matter haloes have scale-heights of tens of kpc, while the stellar disc thickness is $\sim 0.3\,\text{kpc}$. Within $\pm 1.5\,\text{kpc}$, the halo density varies by less than a percent while the stellar density falls by more than an order of magnitude. So the gradient of the gravitational acceleration far from the disc is, to good accuracy, set entirely by the dark matter — exactly the leverage point we used.

---

Part (vi) — Time since the perturbation

Setting up phase mixing

At the moment of perturbation (time $t$ ago), the bar-shaped overdensity put stars of all the relevant amplitudes at nearly the same orbital phase. Since then, each star has advanced in phase at its own angular frequency $\omega(z_m)$, so its phase today is

$$\varphi(z_m,\,t) \;=\; \varphi_\text{init}\;+\;\omega(z_m)\,t \pmod{2\pi}.$$

Two stars at different amplitudes therefore differ in phase by

$$\Delta\varphi \;=\; \big[\omega(z_{m,1})-\omega(z_{m,2})\big]\,t,$$

without any modular ambiguity if we count windings, which is exactly what the spiral lets us do: each axis crossing along the spiral advances $\Delta\varphi$ by $\pi/2$, equivalently $\pi$ between consecutive $v_z=0$ crossings.

Periods from the matching parabola

Using the $\Phi(z_m)$ values from part (iii) and $T(z_m)=\pi z_m\sqrt{2/\Phi(z_m)}$, then converting via $1\,\text{kpc/(km/s)}=978\,\text{Myr}$:

$z_m$ [kpc]	$\Phi$ $[(\text{km/s})^2]$	$T$ [Myr]	$\omega$ [Myr$^{-1}$]
0.18	97	79.4	0.0791
0.32	202	97.8	0.0643
0.50	346	116.7	0.0539
0.70	549	129.8	0.0484
1.00	808	152.8	0.0411
1.40	1240	172.7	0.0364

Solving for $t$

Pick the two amplitudes that are best-determined and most widely separated in $\omega$. Between the $v_z=0$ crossings at $z_m=0.18$ and $z_m=1.40\,\text{kpc}$, the spiral makes 5 half-turns (counting consecutive $v_z=0$ crossings on alternating sides), so

$$\Delta\varphi \;=\; 5\pi.$$

Equating to $(\omega_\text{in}-\omega_\text{out})\,t$:

$$t \;=\; \frac{5\pi}{\omega(0.18)-\omega(1.40)} \;=\; \frac{5\pi}{0.0791 - 0.0364}\,\text{Myr} \;=\; \frac{15.71}{0.0427}\,\text{Myr} \;\approx\; 370\,\text{Myr}.$$

Cross-checks with other pairs

The same procedure with different pairs gives:

inner $z_m$	outer $z_m$	half-turns $N$	$t = N\pi/\Delta\omega$ [Myr]
0.18	1.00	4	330
0.32	1.00	3	406
0.32	1.40	4	451
0.50	1.00	2	491
0.50	1.40	3	539

The values cluster around $\sim 400\,\text{Myr}$ but spread by a factor of $\sim 1.5$, signalling that the matching-parabola formula and our crossing readings together have $\sim 30\%$ uncertainty. A reasonable compromise:

$$\boxed{\,t \;\approx\; 4\times 10^{2}\,\text{Myr}\,}.$$

This places the disturbance during the late Cenozoic from Earth’s standpoint and is consistent with the passage of the Sagittarius dwarf galaxy through the Milky Way disc inferred from independent stellar-stream evidence.

Educational remark

The spread among the entries above is physically meaningful, not just numerical noise: the matching-parabola period systematically under-estimates the true period for non-parabolic potentials (in the limit of constant gravitational force, by about $27\%$), and the under-estimate is largest for the most anharmonic — i.e. outermost — orbits. So pairs that include $z_m=1.4\,\text{kpc}$ tend to overestimate $t$. The original Antoja et al. (2018) analysis, using the full nonlinear period, narrowed the answer to the range $300$–$900\,\text{Myr}$, with subsequent refinements (including Guo et al. 2024) settling near $\sim 500\,\text{Myr}$.

---

Summary of results

Part	Answer
(i)	$a_z = -4\pi G\rho_0 z$
(ii)	$T = \sqrt{\pi/(G\rho_0)} \approx 84\,\text{Myr}$ for $\rho_0=0.1\,M_\odot/\text{pc}^3$
(iii)	$\Phi(0.18)\approx 97,\;\Phi(0.32)\approx 202,\;\Phi(0.50)\approx 346,\;\Phi(0.70)\approx 549,\;\Phi(1.00)\approx 808,\;\Phi(1.40)\approx 1240$ in $(\text{km/s})^2$
(iv)	$\rho_0 \approx 0.11\,M_\odot/\text{pc}^3$
(v)	$\rho_\text{DM} \approx 0.011\,M_\odot/\text{pc}^3$
(vi)	$t \approx 4\times 10^{2}\,\text{Myr}$