I usually write about physics here, but I felt I had to say something about what has recently been reported at Cornell University.
Reading about the alleged sexual assault and everything that has followed has honestly been disturbing. I don’t know all the details, and I don’t want to pretend that I do. The case is still being investigated, so I won’t speculate about what happened beyond what has been publicly reported.
My thoughts are with the survivor and with everyone who has been affected by this. I also think about people who may read about this and be reminded of things they have experienced themselves.
Sexual violence has no place anywhere. No student should have to feel unsafe at their university or in their own community.
We must stand up against it. We have to listen to survivors, support them, and speak up when something is wrong. We cannot just look away because something happened to people we don’t know.
I hope the investigation is fair, thorough, and transparent, and that everyone affected is treated with the dignity and support they deserve.
Here is a number that should bother you. In the solar wind near the Earth’s orbit the proton density is about $5$ per cubic centimetre and the temperature is around $10^{5}$ K. Work out how far a proton travels before Coulomb collisions turn it appreciably off course and you get something like $10^{13}$ cm, most of the way to an astronomical unit. Yet stream-interaction regions, shocks and flux ropes in that wind are routinely modelled with fluid equations, often on scales comparable to this or smaller. The hot gas that fills a galaxy cluster is worse. There the mean free path is some twenty kiloparsecs, comparable to the turbulent eddies and cold fronts that people model with magnetohydrodynamics every day. By the standard we apply to air or water, these plasmas have no business behaving like fluids.1The numbers in this article are representative, order-of-magnitude values. A Coulomb mean free path depends on the species, the temperature and the relaxation process used to define it. I use the thermal estimate common in the cluster literature, $\lambda_{\rm mfp}\approx23\,{\rm kpc}\,(T/10^{8}\,{\rm K})^{2}(n_e/10^{-3}\,{\rm cm^{-3}})^{-1}$ (Sarazin 1986), with the Coulomb logarithm adjusted for each environment. That gives about $10^{13}$ cm for the solar wind (with $L=1$ AU), $2\times10^{12}$ cm for the warm ionised medium ($L=100$ pc) and $20$ kpc for the cluster ($L=200$ kpc). Other conventions give values a few times smaller. In a real application $L$ should be the relevant gradient scale, not the size of the whole system. The NRL Plasma Formulary (2023) lists species-dependent collision rates and plasma scales.
They behave like fluids anyway. Magnetohydrodynamics predicts the structure of the solar wind, the propagation of Alfvén waves through it, the shapes of cold fronts in clusters, the collimation of astrophysical jets and the stability of accretion discs, and it does so well enough that its failures are informative rather than catastrophic. The electromagnetic field obviously has something to do with this. The question is what, exactly, the field supplies, and what else is needed before a closed set of fluid equations can be justified.
Coulomb mean free path $\lambda_{\rm mfp}$ (filled) and characteristic scale $L$ (open) on a logarithmic axis. The Knudsen number $\mathrm{Kn}=\lambda_{\rm mfp}/L$ is about $0.7$ in the solar wind and $0.1$ in the cluster, against $7\times10^{-9}$ in the warm ionised medium.
Start with why a neutral gas needs collisions at all, since the requirement is easy to take for granted. A fluid description assigns a density, a velocity and a pressure to each point and claims that their evolution closes on itself. That is a strong claim. The underlying object is a distribution function $f(\mathbf{x},\mathbf{v},t)$ in six-dimensional phase space, and to shrink it down to a handful of local fields you need some way of writing the higher velocity moments in terms of the ones you keep. Collisions provide it. They push the local distribution toward a Maxwellian on a timescale $\nu^{-1}$. If that is short compared with the time over which the flow changes, and the mean free path $\lambda$ is short compared with the scale $L$ over which gradients vary, the distribution stays close to Maxwellian everywhere and the fluid equations close. The controlling parameter is the Knudsen number $$\mathrm{Kn}=\frac{\lambda}{L},$$ and the local collisional fluid limit needs $\mathrm{Kn}\ll1$, along with collisions that are fast compared with the dynamics. For air at room temperature and atmospheric pressure $\lambda$ is about $70$ nanometres, so a flow that varies over $10$ cm has $\mathrm{Kn}\sim7\times10^{-7}$. Nobody worries about this for a river. People do worry about it in microfluidics and in the thin upper atmosphere.2The transport coefficients of a neutral gas come out of the same expansion: viscosity and thermal conduction are first-order corrections to the Maxwellian. Kinematic viscosity and thermal diffusivity both scale as the mean free path times the thermal speed; dynamic viscosity and conductivity carry extra factors of density and heat capacity. Once $\mathrm{Kn}$ is no longer small the Chapman–Enskog expansion stops being controlled and you need kinetic theory.
In the solar wind and the cluster, $\mathrm{Kn}$ is $0.1$ or larger. What rescues the situation is that a plasma can organise its motion collectively, through the electromagnetic field, instead of through binary encounters. To see how, it helps to build the relevant scales from scratch, because the relations between them do most of the explaining.
The most basic fact about a plasma is that it hides its own charges. Drop a test charge $q$ into a neutral plasma of density $n$ and electron temperature $T$. The electrons are light and rearrange quickly, so let them reach thermal equilibrium in the resulting potential while the much heavier ions stay put as a uniform neutralising background. Away from the test charge, Poisson’s equation reads $$\nabla^{2}\phi=-4\pi\rho_q=-4\pi e\left(n_i-n_e\right),$$ with $n_i=n_0$ fixed and the electrons Boltzmann-distributed, $n_e=n_0\exp(e\phi/k_BT)$. For potentials small compared with the thermal energy, $|e\phi|\ll k_BT$, expand the exponential to first order, so that $n_i-n_e\simeq-n_0e\phi/k_BT$ and $$\nabla^{2}\phi=\frac{4\pi n_0e^{2}}{k_BT}\,\phi\equiv\frac{\phi}{\lambda_D^{2}}.$$ This is a screening equation, not Laplace’s equation, and its spherically symmetric solution around a point charge is $$\phi(r)=\frac{q}{r}\,e^{-r/\lambda_D},\qquad \lambda_D=\sqrt{\frac{k_BT}{4\pi n_0e^{2}}}.$$ The Coulomb potential dies off exponentially beyond $\lambda_D$, the electron Debye length.3I use Gaussian cgs units throughout. The symbol $n$ is the electron number density, equal to the proton density in a quasineutral hydrogen plasma; $\rho_q$ is charge density, while $\rho$ in the fluid equations below is mass density. At the origin the screening equation has the point source $-4\pi q\delta^{(3)}(\mathbf r)$, which fixes the prefactor of the Yukawa potential. If the ions are allowed to respond too, they add their own screening term, and in general $\lambda_D^{-2}=\sum_s 4\pi n_sq_s^{2}/k_BT_s$. For singly charged ions at the electron temperature this gives $\lambda_D=\sqrt{k_BT/8\pi ne^{2}}$, smaller by $\sqrt2$. Which one applies depends on the timescale: faster than the ion plasma period, the ions haven’t had time to move. I stick with the electron Debye length because it is the one tied to the plasma frequency and to the classical Coulomb cutoff below. The literature uses both, so check which one a quoted number assumes.
Debye shielding. The bare Coulomb potential (dashed) and the screened potential (solid), with the same normalisation. The two marked points show that at $r=\lambda_D$ the screened potential has fallen to $1/e$ of the Coulomb value. The cutoff is exponential, not a hard edge.
In practical units, $\lambda_D\simeq743\,\sqrt{T[\mathrm{eV}]/n[\mathrm{cm^{-3}}]}$ cm. The warm ionised component of the interstellar medium, with $n\simeq0.3$ and $T\simeq8000$ K, has a Debye length of about $11$ metres. The solar wind at $1$ AU comes out at about $10$ metres, and a coronal plasma with $n\simeq5\times10^8\,\mathrm{cm^{-3}}$ and $T\simeq10^6$ K at about $3$ millimetres. Set against the systems these plasmas fill, those numbers are absurdly small. On scales much larger than $\lambda_D$, and for dynamics slow compared with the plasma’s response, quasineutrality holds to excellent accuracy. That is the first thing a plasma does that a neutral gas cannot: it polices its own charge separation.
How many particles take part in that policing? Count the number inside a Debye sphere, $$N_D=n\cdot\frac{4\pi}{3}\lambda_D^{3}.$$ For the warm interstellar medium this is about $2\times10^{9}$, for the solar wind at $1$ AU about $2\times10^{10}$, and for an intracluster plasma with $n\simeq10^{-3}\,\mathrm{cm^{-3}}$ and $T\simeq10^8$ K about $4\times10^{16}$. These are enormous numbers, and they are the defining feature of the plasmas in question. Each particle interacts, weakly, with billions of others at once through their overlapping Coulomb fields, and no single neighbour matters much. A plasma with $N_D\gg1$ is called weakly coupled, and the name is precise: each pairwise interaction is weak, while the collective ones are not.
The second thing a plasma does is oscillate. Take a slab of plasma and shift all the electrons by a small distance $x$ relative to the ions, which stay put because they are nearly two thousand times heavier. The shift leaves a surface charge density $\sigma=nex$ on one face and the opposite on the other, and the field between them, $E=4\pi nex$, pulls the electrons back. Newton’s law for an electron reads $$m_e\ddot{x}=-eE=-4\pi ne^{2}x,$$ a harmonic oscillator with $$\omega_p=\sqrt{\frac{4\pi ne^{2}}{m_e}}.$$ In practical units $\omega_p\simeq5.6\times10^{4}\sqrt{n[\mathrm{cm^{-3}}]}$ radians per second, so in the warm interstellar medium the plasma frequency is about $5$ kHz and the period about $0.2$ ms. Any charge separation is undone on a timescale of order $\omega_p^{-1}$. In this idealised model the electrons overshoot and keep oscillating; it takes damping, which we have left out, to make them settle.
The plasma frequency. Shifting the electrons (dashed) a distance $x$ relative to the ions (shaded) exposes surface charges $\pm nex$. The field $E=4\pi nex$ between them pulls the electrons back, and they oscillate at $\omega_p$.
These two quantities are not independent, and the connection between them is the first real insight. The electron thermal speed is $v_{\rm th,e}=\sqrt{k_BT/m_e}$, and dividing it by the plasma frequency gives $$\frac{v_{\rm th,e}}{\omega_p}=\sqrt{\frac{k_BT}{m_e}}\cdot\sqrt{\frac{m_e}{4\pi ne^{2}}}=\sqrt{\frac{k_BT}{4\pi ne^{2}}}=\lambda_D.$$ The Debye length is the distance a thermal electron covers in one response time $\omega_p^{-1}$ (not in a full period, $2\pi/\omega_p$). Screening and oscillation are the same balance between thermal motion and electrostatic restoring force, one seen in space and the other in time.
Now for collisions, which in a plasma are stranger than they look. Two charges scatter through the Coulomb interaction, and the impact parameter that gives a large deflection is set by equating the potential energy at closest approach to the thermal kinetic energy, $$\frac{e^{2}}{b_{90}}\sim k_BT\qquad\Longrightarrow\qquad b_{90}\sim\frac{e^{2}}{k_BT}.$$ In the warm interstellar medium that is about $2\times10^{-7}$ cm, a tiny target. If large-angle encounters were the whole story, the cross-section would be $\pi b_{90}^{2}$ and the mean free path would be gigantic. But the Coulomb force is long-ranged, and a particle passing at $b\gg b_{90}$ still picks up a small deflection, of order $b_{90}/b$. These small kicks are far more numerous, since the number of encounters with impact parameter near $b$ grows as $b\,db$, and they add up as a random walk in angle. Summing the mean-square deflection rate over impact parameters gives, up to numerical factors, $$\frac{d\langle\Delta\theta^{2}\rangle}{dt}\sim nv\int\left(\frac{b_{90}}{b}\right)^{2}b\,db=nvb_{90}^{2}\ln\left(\frac{b_{\max}}{b_{\min}}\right).$$ The integral diverges logarithmically at both ends. It is cut off below at $b_{90}$, where the small-angle approximation fails, and above at $\lambda_D$, where screening switches the interaction off. The effective cross-section is therefore $\sigma\sim b_{90}^{2}\ln\Lambda$ with $$\Lambda=\frac{\lambda_D}{b_{90}},$$ and the mean free path is $$\lambda_{\rm mfp}=\frac{1}{n\sigma}\sim\frac{(k_BT)^{2}}{ne^{4}\ln\Lambda}\propto\frac{T^{2}}{n\ln\Lambda}.$$ The logarithm changes slowly, so the scaling is usually quoted as $T^2/n$. Numerical prefactors differ between treatments by up to an order of magnitude depending on which species and which relaxation process you follow, so trust the scaling more than the coefficient.4The standard results are the Spitzer collision times, which distinguish electron–electron, electron–ion and ion–ion relaxation and give different coefficients for momentum transfer, energy exchange and temperature equilibration. The classical estimate $b_{\min}\sim b_{90}$ also has to be replaced by the de Broglie wavelength when that is larger, as it is for electrons in the hottest cluster gas. Notice, too, that the large-$b$ divergence is cut off by screening, so the Debye length enters the collision rate itself. Even in principle, collisions in a plasma are not independent of its collective behaviour.
Coulomb scattering. Attractive encounters, with the impact parameter $b$ measured from the incoming asymptote. At $b\sim b_{90}$ the deflection is large; grazing encounters give $|\Delta\theta|\sim b_{90}/b$, whatever the sign of the charges. Adding up their mean-square deflections gives the Coulomb logarithm, $\ln\Lambda\simeq\ln(\lambda_D/b_{90})$.
The $T^{2}/n$ scaling settles half of the paradox we started with. Hot, thin plasmas have long mean free paths because both factors push the same way. Whether a given plasma counts as collisionless still depends on the length and time scales in the problem, and a great many astrophysical problems land on the collisionless side. What the scaling doesn’t explain is why those plasmas behave like fluids.
The explanation starts from a relation that is easy to miss. Put the electron Debye length into the classical estimate for $\Lambda$ and use $k_BT=4\pi ne^{2}\lambda_D^{2}$: $$\Lambda=\frac{\lambda_D}{b_{90}}\sim\frac{\lambda_Dk_BT}{e^{2}}=\frac{\lambda_D\cdot4\pi ne^{2}\lambda_D^{2}}{e^{2}}=4\pi n\lambda_D^{3}=3N_D.$$ So the Coulomb logarithm is, give or take the choice of cutoffs, the logarithm of the number of particles in a Debye sphere. That ties together two facts we had been treating separately. Large $N_D$ means that individual Coulomb interactions are weak and that the collective field is smooth and well defined. It does not stop collisions mattering over long enough times. Whether a particular problem is collisional is decided by $\nu/\omega$ and $\lambda_{\rm mfp}/L$, not by $N_D$. The astrophysical plasmas here are deep in the weakly coupled regime, with $\ln\Lambda$ between about $20$ and $40$.
So Coulomb relaxation is slow, yet the plasma still reacts to charge separation within $\omega_p^{-1}$ and over $\lambda_D$. That is a kind of local coupling, but it only constrains the charge density. It says nothing about how momentum gets around. For that we need the magnetic field.
A nonrelativistic charged particle in a uniform magnetic field $\mathbf{B}$, with nothing else acting on it, moves freely along the field and gyrates around it with angular frequency $$\Omega_c=\frac{|q|B}{mc}$$ and radius $r_g=v_\perp/\Omega_c$. The sign of $q$ sets the sense of rotation. In practical units the proton gyrofrequency is $\Omega_{cp}\simeq9.6\times10^{3}\,B[\mathrm{G}]$ radians per second, and the electron’s is larger by the mass ratio. In the warm interstellar medium, with $B\simeq5\,\mu$G, a thermal proton moving at $v_{\rm th,p}=\sqrt{k_BT/m_p}$ gyrates on a radius of about $170$ km.
Gyromotion. A helical orbit in a uniform field $\mathbf B$, seen side-on. The marked particle sits one gyroradius $r_g$ from its guiding centre. Motion along the field is unconstrained.
Compare that with the collisional mean free path in the same gas, about $2\times10^{7}$ km. The gyroradius is smaller by a factor of $10^{5}$. In the solar wind, with a field of about $5$ nT, the ratio is around a million. In a cluster with a microgauss field, a thermal proton gyrates on a radius of about $10^{5}$ km while travelling some $20$ kpc between collisions, a ratio of $10^{12}$ to $10^{13}$. Across the field, its orbit is the size of a large planet. Along the field, it can cross a good fraction of the cluster core before anything deflects it. The Lorentz force ties each particle to its field line far more tightly than collisions ever could, but it does nothing to stop the particle streaming along that line or its guiding centre drifting.
It is tempting to declare victory here: across the field, gyration does the job that collisions do in a neutral gas, so the fluid equations follow. That is too quick, and it pays to separate what magnetisation gives you from what it doesn’t.
What it gives you is scale separation. When the fields change slowly compared with the gyrofrequency and over distances large compared with $r_g$, the fast gyration can be averaged away, leaving guiding centres that drift across the field and stream along it. A collision then moves a particle across the field by about $r_g$ rather than $\lambda_{\rm mfp}$, although the resulting transport still depends on the scattering rate and the field geometry. That is a big simplification, and it underpins a family of reduced kinetic theories that work even without collisions.
What it doesn’t give you is closure. A fluid theory needs the higher moments written in terms of the ones it keeps, and gyration supplies no such rule. Nothing confines particles along the field, so the parallel distribution can be far from Maxwellian, carry a heat flux that no local formula captures, and trade energy with waves through Landau damping. The pressure can become anisotropic, because $p_\perp$ and $p_\parallel$ respond differently to compression. At finite $k_\perp r_g$ there are finite-Larmor-radius corrections that ordinary MHD leaves out. Scale separation gets you a gyro-averaged kinetic theory. Getting from there to single-fluid MHD takes further, separately justified approximations for the pressure, the heat flux and the moments above them. The last part of this article is about how a plasma can arrange some of those conditions for itself.5For a systematic treatment see Schekochihin et al. (2009), Astrophysical gyrokinetics. In their low-frequency, anisotropic ordering, the Alfvénic part of the turbulent cascade obeys reduced-MHD equations even at collisionless scales, while the compressive fluctuations need a kinetic description. Whether a fluid reduction is valid can depend on which fluctuations you are following.
Put the scales together and you get a hierarchy, which is the real subject of this article. For the warm interstellar medium, $$\lambda_D\sim11\ \mathrm{m}\ \ll\ r_g\sim170\ \mathrm{km}\ \ll\ \lambda_{\rm mfp}\sim2\times10^{7}\ \mathrm{km}\ \ll\ L\sim100\ \mathrm{pc},$$ spanning some seventeen orders of magnitude from the Debye length to the system size. The ordering matters more than the individual numbers. Below $\lambda_D$, charge separation has to be kept. Between $\lambda_D$ and $r_g$, low-frequency fluctuations are quasineutral but still resolve the gyration. On perpendicular scales well above $r_g$, and at frequencies well below $\Omega_c$, gyration averages out, though guiding centres still stream along the field. Only when the parallel scales are also much longer than $\lambda_{\rm mfp}$, and the evolution is slow compared with collisions, does a local collisional closure become available. So the hierarchy shows that a gyro-averaged description becomes useful far below the collisional mean free path, provided the timescales cooperate too. The step to a single fluid needs an argument about the distribution function, and the ordering of scales doesn’t provide one.
Hierarchy of scales in the warm ionised interstellar medium, on a logarithmic length axis. Gyro-averaging (heavy bar) needs both spatial and temporal separation from the gyromotion. A local collisional closure also needs parallel gradients longer than $\lambda_{\rm mfp}$ and evolution slower than collisions. The ordering alone does not give single-fluid MHD.
With the hierarchy in hand, it’s worth writing down exactly what magnetohydrodynamics assumes, because the electromagnetic orderings and the fluid closure are separate ingredients. The equations are the usual continuity and momentum equations with a Lorentz force, plus an induction equation for the field: $$\begin{aligned}\frac{\partial\rho}{\partial t}+\nabla!\cdot(\rho\mathbf v)&=0,\ \rho\frac{D\mathbf v}{Dt}&=-\nabla p+\frac{1}{c}\mathbf J\times\mathbf B,\ \frac{\partial\mathbf B}{\partial t}&=\nabla\times(\mathbf v\times\mathbf B)+\eta\nabla^2\mathbf B.\end{aligned}$$ Here $D/Dt=\partial_t+\mathbf v\cdot\nabla$, $\rho$ is the mass density, $p$ the total scalar thermal pressure and $\eta$ a uniform magnetic diffusivity, with $\nabla\cdot\mathbf B=0$ and $\mathbf J=(c/4\pi)\nabla\times\mathbf B$. Several orderings sit behind these equations. Scales much larger than $\lambda_D$ and frequencies much lower than $\omega_p$ give quasineutrality. Scales larger than the gyroradii and frequencies below the gyrofrequencies allow gyro-averaging. Slow, nonrelativistic dynamics let us drop the displacement current. None of this gives $\mathbf E=-\mathbf v\times\mathbf B/c$. That ideal Ohm’s law is a separate approximation, which needs the resistive, Hall, electron-pressure and electron-inertia terms of the generalised Ohm’s law to be small. Keeping only the resistive term gives $\mathbf E+\mathbf v\times\mathbf B/c=(4\pi\eta/c^2)\mathbf J$, which is the form used above. None of these orderings requires $\lambda_{\rm mfp}\ll L$, and none of them closes the equations either. Moments of the kinetic equation form an infinite hierarchy in which each moment depends on the next, so scalar-pressure MHD needs an energy equation and some justified treatment of the heat flux and viscous stress. The simplest non-dissipative choice is $D(p/\rho^\gamma)/Dt=0$, with $\gamma=5/3$ for a monatomic gas. Collisions can justify it, and so, sometimes, can effective scattering or a restricted dynamical regime. A collisionless plasma is under no obligation to obey it.6Once resistivity is kept, the energy equation must include Ohmic heating, and other transport or heating terms where they matter; the adiabatic law is the ideal limit. Dropping the Hall and electron-inertia terms brings in the ion and electron inertial lengths, so the four scales in the figure above are not a complete checklist for MHD. See Fitzpatrick’s derivation of the MHD equations, and, for the Hall extension, Acheritogaray et al. (2011).
The induction equation carries most of the physics that makes magnetised plasmas behave as they do. The ratio of the two terms on its right-hand side is the magnetic Reynolds number, $$R_m=\frac{vL}{\eta},$$ where $\eta$ is the magnetic diffusivity, proportional to the electrical resistivity. In a fully ionised plasma the Spitzer value scales as $T_e^{-3/2}$, apart from the Coulomb logarithm. For the warm ionised medium at $T\simeq8000$ K this gives $\eta\sim10^{7}\ \mathrm{cm^{2}\,s^{-1}}$, and with $v\sim10\ \mathrm{km\,s^{-1}}$ and $L\sim1$ pc, $R_m$ comes out around $10^{17}$. If the other nonideal terms are small too, flux freezing is an excellent approximation on large scales: the magnetic flux through any surface moving with the fluid is conserved, and fluid elements that start on the same field line stay on the same field line. Thin current sheets are the exception, since their local $R_m$ can be far smaller. That one property accounts for much of what magnetic fields do in astrophysics, from the winding up of galactic fields by differential rotation to the launching of jets.
$R_m=\dfrac{vL}{\eta}\gg1$
Flux freezing. In ideal MHD, fluid elements $A$ and $B$ that start on the same field line stay on it however the flow deforms the line. Large $R_m$ makes resistive slippage negligible on the scale considered; the other nonideal terms have to be small too.
It also means the field is not a fixed backdrop. Currents in the plasma reshape it, and it pushes back. Tension along the field lines resists bending, the magnetic pressure $B^{2}/8\pi$ resists compression across them, and the ratio of thermal to magnetic pressure is $$\beta=\frac{8\pi p}{B^{2}},\qquad p=p_e+p_i.$$ For a hydrogen plasma with equal electron and ion temperatures, $p\simeq2nk_BT$. The warm ionised example then has $\beta\simeq0.7$. In the solar wind $\beta$ is typically of order unity, though it varies a lot. In the cluster it is in the hundreds, so magnetic pressure hardly matters to the overall force balance. The field still matters, though. It steers heat transport, sets the thresholds for the pressure-anisotropy instabilities described below, and shapes local structures.
So far this reads as a success story, so it is worth being clear where it stops. Gyration limits motion across the field, but a uniform field does nothing to stop particles streaming along it. Many of the important departures from MHD come from that asymmetry, though not all of them. Nonuniform fields reflect particles through the mirror force, and Hall effects, electron inertia and finite-Larmor-radius corrections impose limits of their own.
Heat transport is the clearest case. Start with the classical collisional estimate, which makes the geometry easy to see. An electron travels a distance of order $\lambda_{{\rm mfp},e}=v_{{\rm th},e}\tau_e$ along the field in one deflection time $\tau_e$, while each scattering shifts its guiding centre across the field by about $r_{g,e}$. In the strongly magnetised limit $\Omega_{ce}\tau_e\gg1$ the ratio of conductivities is, up to numerical factors, $$\frac{\kappa_{e\perp}}{\kappa_{e\parallel}}\sim\frac{1}{(\Omega_{ce}\tau_e)^2}=\left(\frac{r_{g,e}}{\lambda_{{\rm mfp},e}}\right)^2\ll1.$$ These are electron quantities. When the electron and ion temperatures are similar, electrons carry most of the parallel heat flux, so the proton gyroradius from earlier doesn’t belong in this estimate. In the collisional regime the parallel conductivity is the Spitzer–Braginskii value, as long as the electron mean free path is short compared with the parallel temperature-gradient scale. In a weakly collisional plasma such as the intracluster medium, the heat flux can instead become nonlocal or be capped by kinetic instabilities, and tangled fields and turbulence change how heat moves across the mean field. So the classical ratio shows how extreme the microscopic anisotropy is; it is not a suppression factor to plug into a cluster model. What survives is that the field geometry decides where heat can go. A cold-front model with isotropic conduction misses this, and the observed sharpness of cold fronts has led to models in which the field is draped along the front.7For the classical coefficients and where they apply, see Fitzpatrick, Classical Closure Scheme. For kinetic limits on the parallel electron heat flux, see Roberg-Clark et al. (2018), Wave generation and heat flux suppression in astrophysical plasma systems. For draping at cold fronts, see Lyutikov (2006), Magnetic draping of merging cores and radio bubbles in clusters of galaxies.
Classical anisotropic electron heat transport. Along the field, an electron moves about $\lambda_{{\rm mfp},e}$ between scatterings. Across it (right, $\mathbf B$ into the page), each scattering shifts the guiding centre by about $r_{g,e}$: the solid and dashed orbits meet at the point where the electron scattered. The ratio holds for strongly magnetised, collisional transport.
A subtler failure is that the pressure stops being a single number. A scalar pressure is justified when collisions keep the velocity distribution nearly isotropic. Without them there is no reason for motion along $\mathbf{B}$ and motion across it to share a temperature. The natural description has two pressures, $p_\parallel$ and $p_\perp$, and the adiabatic invariants of gyromotion tell us how each one evolves. The magnetic moment $\mu=m v_\perp^{2}/2B$ is conserved when changes are slow, so strengthening the field raises $p_\perp$, while squeezing plasma along a flux tube raises $p_\parallel$. Following both through gives the double-adiabatic, or CGL, relations, $$\frac{d}{dt}\left(\frac{p_\perp}{\rho B}\right)=0,\qquad \frac{d}{dt}\left(\frac{p_\parallel B^{2}}{\rho^{3}}\right)=0,$$ which replace the single adiabatic law of ordinary gas dynamics. They come from taking moments of the kinetic equation in the long-wavelength, low-frequency limit, ignoring collisions, and closing the system by setting the heat fluxes to zero. That last step is the weak point.8The closure is due to Chew, Goldberger and Low (1956). Its problem is that zero heat flux is not generally justified in a collisionless plasma, where particles stream along field lines and carry energy with them. There is no universal local closure for that transport. Landau-fluid closures, which mimic collisionless phase mixing and Landau damping in the higher moments, are a useful compromise in the regimes they are built for; see Hammett and Perkins (1990).
What happens next is the most interesting part of the story, and it brings back a version of the fluid picture from an unexpected direction. Pressure anisotropy cannot grow without limit, because a plasma with $p_\parallel$ well above $p_\perp$ is unstable. Particles streaming along a curved field line exert a centrifugal force on it that helps the bend grow, and once the pressure anisotropy beats the magnetic tension the field buckles. In the long-wavelength fluid description the threshold is $$p_\parallel-p_\perp>\frac{B^{2}}{4\pi}.$$ This is the firehose instability, named for the way an unsecured hose whips about. The opposite anisotropy is unstable too. When $p_\perp$ exceeds $p_\parallel$ by enough, small dips in the field strength trap particles and deepen. That is the mirror instability, with a threshold for a simple gyrotropic distribution of approximately $$\frac{p_\perp}{p_\parallel}-1>\frac{1}{\beta_\perp},$$ where $\beta_\perp=8\pi p_\perp/B^2$. Precise kinetic thresholds depend on the species and the shape of the distribution. How fast either instability grows depends on how far past threshold the plasma has been pushed, but both can respond much faster than the large-scale flow that drives them. They settle the anisotropy in different ways. Firehose fluctuations scatter particles in pitch angle. Mirror modes at first regulate the anisotropy by trapping particles as the fluctuations grow, with strong scattering appearing later at the sharp edges of the mirror structures.
Where that scattering happens, it acts like an effective collision. It breaks conservation of the magnetic moment and reshuffles pitch angles at a rate set by the fluctuations, not by Coulomb encounters. Trapping offers a second way to cap the anisotropy without randomising every orbit straight away. The plasma is, in effect, responding to its own distribution function: driving pushes it past a stability boundary, fluctuations grow, and their feedback holds it near marginal stability. At high $\beta$ those boundaries allow only a small fractional difference between $p_\perp$ and $p_\parallel$, although even that small difference can rival the magnetic tension. It isn’t exact isotropy, and it doesn’t fix the heat flux or the higher moments. But it is one important reason a weakly collisional plasma can behave more like a fluid than Coulomb collisions alone would suggest. Where this effective scattering dominates, the corresponding mean free path can be far shorter than the Coulomb one. A full MHD closure still needs its own justification.9The difference between firehose scattering and mirror trapping is shown in Kunz, Schekochihin and Stone (2014), Firehose and Mirror Instabilities in a Collisionless Shearing Plasma. A similar self-regulating loop appears in cosmic-ray transport: cosmic rays streaming along the field excite waves that scatter them. If wave damping is weak, their bulk drift is held near the Alfvén speed; stronger damping lets them stream faster. The loops are analogous, not identical; see Thomas and Pfrommer (2019), Cosmic-ray hydrodynamics.
Pressure-anisotropy regulation (schematic). Isotropy sits at $p_\perp/p_\parallel=1$. Pushing the plasma past the firehose boundary, $p_\parallel-p_\perp>B^2/4\pi$, or the mirror boundary, $p_\perp/p_\parallel-1>1/\beta_\perp$, excites fluctuations that pull it back toward marginal stability. The width of the stable band shrinks as $\beta$ grows.
Two other places show different limits of the simplest fluid picture. The first is magnetic reconnection. Ideal MHD preserves magnetic connectivity, yet solar flares, magnetospheric substorms and accretion-disc coronae all involve changes in field topology. In collisionless reconnection, current sheets thin down to ion or electron kinetic scales, where the terms dropped from ideal MHD become important and the field is no longer frozen in. The large-scale flow can stay fluid-like while a small nonideal region decides how the connectivity changes. Fluid theory isn’t abandoned entirely, though: resistive MHD reconnects field lines perfectly well in a collisional plasma, without the layer reaching kinetic scales.10For the resistive-fluid picture, see Fitzpatrick’s plasma-physics notes, chapter Magnetic Reconnection.
The second is high-energy particles. The gyroradius grows with perpendicular momentum, $r_g=p_\perp c/(|q|B)$, so an ultrarelativistic proton at $10^{15}$ eV in a $1\,\mu$G field has a gyroradius of about a parsec: $$r_g\simeq1.1\,\mathrm{pc}\left(\frac{E}{1\,\mathrm{PeV}}\right)\left(\frac{B}{1\,\mu\mathrm G}\right)^{-1}\sin\alpha,$$ where $\alpha$ is the pitch angle between the momentum and the field. At $10^{20}$ eV the same formula gives around $100$ kpc, several times the diameter of the Galactic disc. Particles like these cannot be folded into the thermal fluid. Their transport has to be treated kinetically, or through a cosmic-ray fluid or diffusion model with its own scattering closure. The same plasma can therefore be a fluid for its thermal particles and a scattering medium for its suprathermal ones. There is no contradiction in that. The fluid approximation is always a statement about a particular population on a particular scale.11The coefficient follows from $pc\simeq E$ and $r_g=p_\perp c/(eB)$ in Gaussian units. Cosmic-ray fluid models with an explicit scattering closure are developed, for example, by Thomas and Pfrommer (2019).
Underneath all of this sits $N_D$, the number of particles in a Debye sphere. When it is large, pair interactions are weak and collective behaviour is well defined, and $\ln\Lambda\simeq\ln(3N_D)$ links the two without settling whether any given problem is collisional. Space is nearly empty, but that doesn’t make each particle a free agent. Dilution weakens binary collisions while leaving the collective field free to organise the plasma, as long as the plasma’s own scales stay small compared with the system. Whether that organised motion can be reduced to a fluid then comes down to scale separation, and to how the kinetic degrees of freedom that remain are kept in check.
References and Footnotes
1
The numbers in this article are representative, order-of-magnitude values. A Coulomb mean free path depends on the species, the temperature and the relaxation process used to define it. I use the thermal estimate common in the cluster literature, $\lambda_{\rm mfp}\approx23\,{\rm kpc}\,(T/10^{8}\,{\rm K})^{2}(n_e/10^{-3}\,{\rm cm^{-3}})^{-1}$ (Sarazin 1986), with the Coulomb logarithm adjusted for each environment. That gives about $10^{13}$ cm for the solar wind (with $L=1$ AU), $2\times10^{12}$ cm for the warm ionised medium ($L=100$ pc) and $20$ kpc for the cluster ($L=200$ kpc). Other conventions give values a few times smaller. In a real application $L$ should be the relevant gradient scale, not the size of the whole system. The NRL Plasma Formulary (2023) lists species-dependent collision rates and plasma scales. ↩︎
2
The transport coefficients of a neutral gas come out of the same expansion: viscosity and thermal conduction are first-order corrections to the Maxwellian. Kinematic viscosity and thermal diffusivity both scale as the mean free path times the thermal speed; dynamic viscosity and conductivity carry extra factors of density and heat capacity. Once $\mathrm{Kn}$ is no longer small the Chapman–Enskog expansion stops being controlled and you need kinetic theory. ↩︎
3
I use Gaussian cgs units throughout. The symbol $n$ is the electron number density, equal to the proton density in a quasineutral hydrogen plasma; $\rho_q$ is charge density, while $\rho$ in the fluid equations below is mass density. At the origin the screening equation has the point source $-4\pi q\delta^{(3)}(\mathbf r)$, which fixes the prefactor of the Yukawa potential. If the ions are allowed to respond too, they add their own screening term, and in general $\lambda_D^{-2}=\sum_s 4\pi n_sq_s^{2}/k_BT_s$. For singly charged ions at the electron temperature this gives $\lambda_D=\sqrt{k_BT/8\pi ne^{2}}$, smaller by $\sqrt2$. Which one applies depends on the timescale: faster than the ion plasma period, the ions haven’t had time to move. I stick with the electron Debye length because it is the one tied to the plasma frequency and to the classical Coulomb cutoff below. The literature uses both, so check which one a quoted number assumes. ↩︎
4
The standard results are the Spitzer collision times, which distinguish electron–electron, electron–ion and ion–ion relaxation and give different coefficients for momentum transfer, energy exchange and temperature equilibration. The classical estimate $b_{\min}\sim b_{90}$ also has to be replaced by the de Broglie wavelength when that is larger, as it is for electrons in the hottest cluster gas. Notice, too, that the large-$b$ divergence is cut off by screening, so the Debye length enters the collision rate itself. Even in principle, collisions in a plasma are not independent of its collective behaviour. ↩︎
5
For a systematic treatment see Schekochihin et al. (2009), Astrophysical gyrokinetics. In their low-frequency, anisotropic ordering, the Alfvénic part of the turbulent cascade obeys reduced-MHD equations even at collisionless scales, while the compressive fluctuations need a kinetic description. Whether a fluid reduction is valid can depend on which fluctuations you are following. ↩︎
6
Once resistivity is kept, the energy equation must include Ohmic heating, and other transport or heating terms where they matter; the adiabatic law is the ideal limit. Dropping the Hall and electron-inertia terms brings in the ion and electron inertial lengths, so the four scales in the figure above are not a complete checklist for MHD. See Fitzpatrick’s derivation of the MHD equations, and, for the Hall extension, Acheritogaray et al. (2011). ↩︎
The closure is due to Chew, Goldberger and Low (1956). Its problem is that zero heat flux is not generally justified in a collisionless plasma, where particles stream along field lines and carry energy with them. There is no universal local closure for that transport. Landau-fluid closures, which mimic collisionless phase mixing and Landau damping in the higher moments, are a useful compromise in the regimes they are built for; see Hammett and Perkins (1990). ↩︎
For the resistive-fluid picture, see Fitzpatrick’s plasma-physics notes, chapter Magnetic Reconnection. ↩︎
11
The coefficient follows from $pc\simeq E$ and $r_g=p_\perp c/(eB)$ in Gaussian units. Cosmic-ray fluid models with an explicit scattering closure are developed, for example, by Thomas and Pfrommer (2019). ↩︎
Newtonian mechanics contains two masses with no obvious reason to be related. One appears in the law of gravitation, $$\mathbf{F}=-\frac{GMm_g}{r^{2}}\hat{\mathbf r},$$ where $m_g$ measures how strongly a body couples to a gravitational field. The other appears in the second law, $$m_i\mathbf a=\mathbf F,$$ where $m_i$ measures resistance to being accelerated by any force at all and has nothing to do with gravity in particular. Combining them gives $$\mathbf a=-\frac{GM}{r^{2}}\frac{m_g}{m_i}\hat{\mathbf r},$$ so the acceleration of a falling body depends on the ratio $m_g/m_i$. Nothing in Newtonian mechanics requires that ratio to be the same for different materials. Experimentally, though, it’s universal to very high precision, so we set $m_i=m_g$ and write $$\mathbf a=-\frac{GM}{r^{2}}\hat{\mathbf r}.$$
The question worth asking is why gravity should be the only interaction whose trajectories don’t care what the test body is made of. For electromagnetism, $$\mathbf a=\frac{q}{m}\mathbf E,$$ so the motion depends explicitly on properties of the particle. A gravitational acceleration apparently doesn’t, and the consequence is that a free-fall trajectory is fixed entirely by an initial position and an initial velocity, with no other information about the falling body required.
That’s a strange property for a force law. It’s a much more ordinary property for geometry: once a starting point and a starting direction are given on a sphere, the corresponding great circle is fixed without asking what is moving along it. This doesn’t prove that gravity is geometric, but it suggests that the universality of free fall may be telling us something about spacetime itself.
The useful next question is how much of gravity can be removed by changing the observer. If the geometric suggestion is pointing in the right direction, then some of what we call a gravitational field may turn out to depend on how spacetime is being described, while some part of it should survive every such change. The way to find out is to start with the simplest possible experiment: a falling object inside an accelerating laboratory.
Objects of different mass and different composition, released together with non-gravitational forces removed, follow the same trajectory and arrive together. Newtonian mechanics contains two logically independent masses, the inertial mass in $m_i\mathbf a=\mathbf F$ and the gravitational mass in the force law, and this records the experimental result that their ratio is the same for every body tested. The acceleration carries no information about the falling object.
The way to test it is to see how much of gravity can be removed by a change of description. Put a laboratory in deep space, far from any mass, and accelerate it with proper acceleration $a$ using a rocket. Someone inside releases a ball; no force acts on it, so relative to the distant stars it moves at constant velocity while the floor accelerates upward to meet it, and to the occupant the ball appears to accelerate toward the floor at $a$. Every ball behaves identically, since the ball is doing nothing and the floor is doing all the work. Now stand the same laboratory on a planet where the gravitational acceleration happens to equal $a$. A released ball accelerates toward the floor at $a$, again independently of mass or composition, this time because of gravity. The two situations have completely different descriptions, one with an inertial ball and an accelerating lab, the other with a static lab and a ball being pulled, and no experiment with falling bodies performed inside a windowless box distinguishes them. Universality of free fall is what guarantees this; if different balls fell differently, dropping two of them would settle the question at once.1This is a reconstruction of the logic rather than the history. Einstein worked at the problem for most of a decade with several wrong turns, and the question here is whether we, knowing modern physics, can reconstruct why this structure makes sense.
A released body has the same local motion in an upward-accelerating laboratory in empty space and in a stationary laboratory with a downward gravitational field of equal magnitude. The equivalence is local; tidal effects distinguish the two over a sufficiently large region.
Running the argument in reverse suggests that the effect of a gravitational field can be removed by accelerating appropriately, which for a field pointing downward means letting the laboratory fall. Cut the cable on an elevator near the Earth’s surface, and the occupant and everything with them accelerates downward at the same rate, so relative to the elevator nothing accelerates. A released ball stays where it’s put, two released balls hover next to each other, and a ball given a gentle push drifts across the elevator in a straight line at constant speed. Mechanics experiments performed inside return the answers of special relativity with no gravitational term in them. It’s tempting to say the gravitational field has been cancelled, and that phrasing causes trouble later, so it’s worth fixing the language now. What has happened is that at a chosen event a freely falling observer can adopt a frame in which the non-gravitational laws take their special-relativistic form and the gravitational acceleration doesn’t appear. Whether anything else survives is a separate question, and the answer to it is the whole point of what follows.
The laboratory, its occupant, and the loose objects all share nearly the same local free-fall acceleration, so nothing accelerates relative to anything else inside the box. Tidal effects appear only when the laboratory is large enough for the gravitational field to vary appreciably across it.
The statement that this holds for all of physics rather than only for mechanics is the equivalence principle, and the version needed here says that in a sufficiently small region of spacetime around any event, local non-gravitational experiments give the same results as in an inertial frame of special relativity with no gravitational field present. The word local is carrying most of the weight in that sentence. Take the freely falling elevator and make it very large, say a thousand kilometres tall, still falling toward the Earth. Release two balls, one near the ceiling and one near the floor. The lower ball sits where $GM/r^{2}$ is larger, so it accelerates slightly more strongly than the elevator average and drifts toward the floor, while the upper ball drifts toward the ceiling. Release two balls side by side instead. Both accelerate toward the Earth’s centre, so their acceleration vectors aren’t quite parallel and the balls drift toward each other. Inside a large enough freely falling laboratory, something is left over.
Local therefore means small enough that this residue sits below whatever precision you care about, and the comparison is quantitative rather than rhetorical. Near the Earth’s surface the tidal gradient is $2GM_\oplus/R_\oplus^{3}\approx3.1\times10^{-6}\ \mathrm{s^{-2}}$, so across a laboratory two metres tall the relative acceleration between top and bottom is about $6\times10^{-6}\ \mathrm{m\,s^{-2}}$, or roughly six parts in $10^{7}$ of $g$. An experiment that can’t resolve that level sees a two-metre freely falling lab as an inertial frame; one that can does not, and needs a smaller lab. There’s no absolute answer to how small counts as local, only an answer relative to the curvature scale and the precision of the measurement. The region also has to be small in time, since an elevator falling for an hour has moved a long way through a field that changes along the way.
So what exactly is the residue? Consider a single freely falling particle first. Construct a small frame that falls along with it, and by the equivalence principle that frame is inertial in a neighbourhood of each event on the worldline. In fact one can do better than event by event: a single frame can be built that stays locally inertial all along the worldline, so an observer riding with the particle finds no gravitational acceleration anywhere in their own history.2Coordinates of this kind are called Fermi normal coordinates. They put the metric in Minkowski form and set the connection coefficients to zero at every point of the central worldline, with corrections appearing only at second order in the distance away from it. Both of those objects are defined later in this post. What curvature obstructs is not the maintenance of such a frame along the curve, but its extension to a finite region surrounding the curve. Nothing local along the way distinguishes a gravitational field from empty space, and the accurate statement is that curvature can’t be detected by measurements confined to an ideally pointlike freely falling frame. A real observer with a gradiometer, or an interferometer with arms of finite length, measures tidal gravity perfectly well; what’s needed is a comparison between points at some finite separation, and the cleanest version of that comparison uses two nearby worldlines.
Release two particles a short distance apart in the freely falling elevator and watch the vector between them. Each particle is individually weightless and each has a local inertial frame in which nothing happens, and yet the separation changes. Two particles released side by side at the same altitude both fall toward the Earth’s centre, and since that centre is a single point their acceleration vectors converge, so the horizontal separation shrinks. Two released one above the other fall along the same radius, the lower one accelerates more strongly, and the vertical separation grows. A small spherical cloud of freely falling particles is therefore stretched radially and squeezed transversely. What no change of coordinates can do is make this covariant relative acceleration disappear. One can even choose coordinates in which both particles keep fixed spatial labels, but then their changing physical separation is carried by the metric rather than by their coordinate positions. It would be too quick to say that all observers agree about the distance between them, since coordinate distance and simultaneity are not invariant and different observers assign different components to the separation vector; what’s invariant is the relative acceleration of the two worldlines described covariantly, and that cannot be transformed away. Making that description precise is the business of the rest of this post.
A local inertial frame can be carried along a single freely falling worldline. Curvature can reveal itself through the relative acceleration of neighbouring freely falling worldlines, governed by the evolution of their separation vector.
In Newtonian language the residue has a clean form. Write the field as the gradient of a potential, $\mathbf{g}=-\nabla\Phi$, so a particle obeys $\ddot{x}^{i}=-\partial_{i}\Phi$. Take two particles, one at $x^{i}$ and one at $x^{i}+\xi^{i}$ with $\xi$ small, subtract their equations of motion and expand the gradient to first order in $\xi$: $$\ddot{\xi}^{i}=-\partial_{i}\partial_{j}\Phi\,\xi^{j}.$$ The relative acceleration of two nearby freely falling particles is governed by the Hessian of the potential. Not the potential itself, which can be shifted by a constant with no physical consequence, and not its gradient, which is exactly what falling absorbs. With $\Phi=-GM/r$ and $\partial_{i}\partial_{j}\Phi=GM(\delta_{ij}/r^{3}-3x_{i}x_{j}/r^{5})$, the eigenvalue along the radial direction is $-2GM/r^{3}$ and in each transverse direction $+GM/r^{3}$, so the radial separation grows and the transverse ones shrink. The three eigenvalues sum to zero, so a small cloud keeps its volume to leading order while its shape distorts, which is just $\nabla^{2}\Phi=0$ holding in the vacuum outside the Earth.
That equation is the target, and it’s also the obstacle. Everything in it is written in a fixed background space with a universal time, and it treats the components of a separation as though they meant something on their own. But the lesson of the falling elevator is that the acceleration of a body depends on who’s describing it: the occupant and the observer on the ground disagree about whether the ball accelerates, and neither of them is making a mistake. If part of gravity is an artefact of the description, then a claim that some other part is not an artefact has to be made in a language where the choice of description can’t influence the answer. Otherwise we have no way of telling a residue from a bookkeeping accident. So the next thing we need isn’t more physics, it’s a way of writing physics in which relabelling events changes nothing.
Suppose Alice describes the falling elevator using coordinates $x^{\mu}$ and Bob uses different coordinates $x’^{\mu}$. Nothing physical has changed; they’ve relabelled events. A coordinate system is a perfectly real and useful thing, a rule for attaching four numbers to each event, and Alice can measure everything she needs with it. What has no invariant meaning is the labels themselves and the components taken with respect to them: the statement that a clock reads $t=5$, or that a vector has components $(3,4)$, is a statement about a labelling scheme, whereas the proper time along a worldline is something a clock actually measures and comes out the same whoever computes it. The laws we’re looking for therefore can’t depend on which labels were used, and this becomes a working requirement: we want statements whose truth doesn’t depend on the coordinate system.
To make that requirement usable we need to be more careful about what a vector is, because the whole difficulty ahead is hiding in a place where the usual picture of an arrow is too loose. Fix an event $P$. Take all smooth curves through $P$ and, for each one, its velocity at $P$. These velocities add and scale like arrows, and together they form a four-dimensional vector space attached to $P$ and to no other point, written $T_{P}M$ and called the tangent space at $P$. A vector at $P$ is an element of $T_{P}M$, and that’s all it is; it isn’t an object stretching from one place to another. Choosing coordinates gives a basis $\mathbf{e}_{\mu}$ of $T_{P}M$, namely the velocities of the four curves along which one coordinate varies and the rest are held fixed, and the components $V^{\mu}$ are the numbers appearing in $\mathbf{V}=V^{\mu}\mathbf{e}_{\mu}$. A covector at $P$ is a linear map from $T_{P}M$ to the real numbers, which in plain terms is a machine that eats one vector and hands back a single number, in a way that respects adding and scaling. The gradient of a function is the standard example, since it takes a velocity and returns the rate at which the function changes along that curve, and its components carry a lower index. A tensor at $P$ is a multilinear map that takes some vectors and some covectors and returns a number, one index per slot; multilinear just means linear in each slot separately, so doubling what you feed into one slot doubles the answer and the other slots don’t notice. That’s the whole content of index notation: an index is a slot, an upper index is a slot for a covector and a lower index a slot for a vector, and the components are the numbers the map returns when fed basis elements. Notice what hasn’t been defined. Nothing so far relates $T_{P}M$ to $T_{Q}M$ at a different event, and that omission is about to become the central problem.
The transformation law follows from this. Under a change of coordinates the basis of $T_{P}M$ changes, so the components change to compensate, giving $$V’^{\mu}=\frac{\partial x’^{\mu}}{\partial x^{\nu}}V^{\nu},\qquad \omega’_{\mu}=\frac{\partial x^{\nu}}{\partial x’^{\mu}}\omega_{\nu},$$ for a vector and a covector, and one such factor per index for anything with more slots, as in $T’^{\mu}{}_{\nu}=\frac{\partial x’^{\mu}}{\partial x^{\alpha}}\frac{\partial x^{\beta}}{\partial x’^{\nu}}T^{\alpha}{}_{\beta}$. This is what makes tensors worth the notation. If $A^{\mu}{}_{\nu}=B^{\mu}{}_{\nu}$ holds in Alice’s system, then multiplying both sides by the same transformation factors gives $A’^{\mu}{}_{\nu}=B’^{\mu}{}_{\nu}$ in Bob’s, so a tensor equation is true in every coordinate system or in none. In particular, if every component of a tensor vanishes for Alice then every component vanishes for Bob, since each of Bob’s components is a linear combination of Alice’s. A tensor equation is therefore a statement about spacetime rather than about labels, which is exactly the kind of statement we were told to look for.
A concrete example makes the danger clear, and it will come back later. Draw an arrow on a flat tabletop. In one set of Cartesian axes it has components $(3,4)$, and its length follows from $dl^{2}=dx^{2}+dy^{2}$ as $5$. Rotate the axes to lie along the arrow and the components become $(5,0)$. Two different pairs of numbers, one arrow, and a length that never moved. Now describe the same tabletop in polar coordinates. The line element reads $$dl^{2}=dr^{2}+r^{2}d\theta^{2},$$ which looks nothing like $dx^{2}+dy^{2}$: there’s a coefficient that varies from place to place, as though the geometry were doing something at large $r$ that it wasn’t doing near the origin. The table is still a plane. What’s changed is that the polar basis vectors are genuinely different objects at different points, unlike the Cartesian ones: $\mathbf{e}_{r}$ points away from the origin, so on opposite sides of the origin it points in opposite directions, and $\mathbf{e}_{\theta}$ gets longer as $r$ grows, because a step $d\theta$ carries you a distance $r\,d\theta$. Take the unit radial vector field as a test case. Its components with respect to the local polar basis are $(1,0)$ at every point of the plane, and yet the arrows at two different points plainly point in different directions. Three things have to be held apart from here on: the geometric vector, which lives in the tangent space at one event; its components, which are just numbers; and the basis, which is what turns the numbers back into the vector. Only the first of the three is coordinate-independent, and equal components in a position-dependent basis don’t mean equal vectors. Hold on to that line element too, because we’ll meet it again in a more dangerous disguise.
Some notation, then. An event is a point of spacetime labelled by four coordinates $x^{\mu}$ with $\mu$ running over $0,1,2,3$, where $x^{0}=ct$, and Greek indices run over all four values while Latin indices $i,j,k$ run over the three spatial ones.3Two smaller conventions, used without further comment below. An index appearing once raised and once lowered is summed over, so $A_{\mu}B^{\mu}$ means $\sum_{\mu=0}^{3}A_{\mu}B^{\mu}$. And I keep $c$ explicit throughout rather than setting it to one, so that every formula can be checked against a Newtonian limit. In special relativity the geometry is carried by the interval between nearby events, $$ds^{2}=\eta_{\mu\nu}\,dx^{\mu}dx^{\nu}=-c^{2}dt^{2}+dx^{2}+dy^{2}+dz^{2},$$ with $\eta_{\mu\nu}$ the Minkowski metric, a constant array with diagonal $(-1,+1,+1,+1)$. That sign pattern is the signature, and I use $(-,+,+,+)$ everywhere below, which matters because many signs flip under the other convention.
In the language just set up, the object defining that interval is a tensor at each event. The metric at $P$ is a symmetric bilinear form $$g_{P}:T_{P}M\times T_{P}M\longrightarrow\mathbb{R},$$ a map taking two vectors at $P$ and returning a number, linear in each argument, symmetric under exchanging them, and non-degenerate, meaning the only vector orthogonal to everything is the zero vector. Its components are its values on the basis, $g_{\mu\nu}=g_{P}(\mathbf{e}_{\mu},\mathbf{e}_{\nu})$, and symmetry leaves ten independent ones at each event, with matrix inverse $g^{\mu\nu}$ defined by $g^{\mu\alpha}g_{\alpha\nu}=\delta^{\mu}_{\nu}$. The interval itself is what that form returns when fed a displacement twice over, which is a number rather than a tensor. Special relativity is the case where coordinates exist making the components $\eta_{\mu\nu}$ everywhere at once. Dropping that requirement and letting them depend on position, $$ds^{2}=g_{\mu\nu}(x)\,dx^{\mu}dx^{\nu},$$ costs nothing in interpretation, because every physical use of the metric is local anyway. A clock carried along a timelike worldline with four-velocity $u$ accumulates proper time at a rate fixed by $g(u,u)=-c^{2}$; the length of a spatial separation is $\sqrt{g(\xi,\xi)}$ for a spacelike vector $\xi$ lying in an observer’s local rest space; two vectors are orthogonal when $g(A,B)=0$; and a direction is null, so a light ray can go that way, when $g(k,k)=0$, which is what fixes the light cones and with them the causal structure. Calling $g_{\mu\nu}$ a gravitational potential isn’t exactly wrong, but it undersells it: it’s the object that tells you what every clock, ruler and light signal in a small neighbourhood will do.
With the metric available, the equivalence principle turns into a statement about coordinates. Around any event $P$ we can choose coordinates in which $$g_{\mu\nu}(P)=\eta_{\mu\nu},\qquad \partial_{\alpha}g_{\mu\nu}(P)=0.$$ The first condition says that the geometry at $P$ is Minkowskian, which it always can be made to be, since a symmetric matrix can be brought to diagonal form with entries $\pm1$ and the signature fixes which. The second says the metric components are stationary at $P$, so to first order in the distance from $P$ the geometry looks flat. I want to resist saying that the first derivatives of the metric are the gravitational acceleration, because that isn’t quite what they are, and the precise version has to wait until we have the object they actually build. The obvious next question is whether we can go further and remove the second derivatives too. We’ll be in a position to answer that properly later, so I’ll leave it standing.
What we can’t do yet is compare anything at two different events, and this is where the tangent-space definition starts to bite. A vector at $P$ belongs to $T_{P}M$ and a vector at a neighbouring event $Q$ belongs to $T_{Q}M$, and these are two different vector spaces. Subtracting one from the other is not a defined operation. But subtraction is exactly what differentiation is: to say how a vector field changes near $P$ we have to take the vector at $Q$, compare it with the vector at $P$, and divide by the separation. The components, being numbers, can of course be subtracted, and that’s what $\partial_{\mu}V^{\nu}$ does, but the answer has no meaning on its own because it mixes a change in the vector with a change in the basis. The polar example already showed this: the radial field has constant components and is not a constant field.
The failure shows up algebraically as well, and it’s worth seeing both faces of it. If $\partial_{\nu}V^{\mu}$ were a tensor it would transform with two factors of the coordinate change. Differentiate the transformation law and see what happens. Using the chain rule for the outer derivative and then the product rule, $$\partial’_{\alpha}V’^{\mu}=\frac{\partial x^{\beta}}{\partial x’^{\alpha}}\partial_{\beta}\left(\frac{\partial x’^{\mu}}{\partial x^{\nu}}V^{\nu}\right)=\frac{\partial x^{\beta}}{\partial x’^{\alpha}}\frac{\partial x’^{\mu}}{\partial x^{\nu}}\,\partial_{\beta}V^{\nu}+\frac{\partial x^{\beta}}{\partial x’^{\alpha}}\frac{\partial^{2}x’^{\mu}}{\partial x^{\beta}\partial x^{\nu}}\,V^{\nu}.$$ The first term is the tensor transformation law we wanted. The second is an obstruction, and it contains second derivatives of the coordinate change, so it vanishes only for linear relabellings. A rotation of Cartesian axes is linear and does no damage, which is why nobody meets this problem in elementary vector calculus; a change to polar coordinates is not linear, and it does. The two faces are the same fact: the obstruction term is precisely the failure to have a canonical way of comparing $T_{P}M$ with $T_{Q}M$.
So the logic is forced. Differentiation requires comparing neighbouring vectors; no comparison exists by default; therefore one has to be supplied as extra structure. That structure is called a connection, and supplying it means saying how the basis vectors change from point to point, $$\nabla_{\nu}\mathbf{e}_{\mu}=\Gamma^{\rho}_{\ \nu\mu}\,\mathbf{e}_{\rho},$$ which declares what it means to carry a vector from one tangent space to the next without changing it. Given that rule, differentiating a field $\mathbf{V}=V^{\mu}\mathbf{e}_{\mu}$ is just the product rule, $$\nabla_{\nu}\mathbf{V}=(\partial_{\nu}V^{\rho})\mathbf{e}_{\rho}+V^{\mu}\nabla_{\nu}\mathbf{e}_{\mu}=\left(\partial_{\nu}V^{\rho}+\Gamma^{\rho}_{\ \nu\mu}V^{\mu}\right)\mathbf{e}_{\rho},$$ so that in components $$\nabla_{\mu}V^{\nu}=\partial_{\mu}V^{\nu}+\Gamma^{\nu}_{\ \mu\lambda}V^{\lambda}.$$ The first term asks how the components changed, the second corrects for the fact that the basis changed too, and the sum is ordinary differentiation with the change of basis subtracted off. For a lower index the correction comes with the opposite sign, $\nabla_{\mu}V_{\nu}=\partial_{\mu}V_{\nu}-\Gamma^{\lambda}_{\ \mu\nu}V_{\lambda}$, and a general tensor picks up one $\Gamma$ term per slot.
The coefficients themselves can’t be tensors, and the reason is structural rather than accidental. For $\nabla_{\mu}V^{\nu}$ to transform correctly, the $\Gamma$ term has to cancel the obstruction found above, so the $\Gamma$ must transform with an inhomogeneous piece of their own, $$\Gamma’^{\rho}_{\ \mu\nu}=\frac{\partial x’^{\rho}}{\partial x^{\lambda}}\frac{\partial x^{\alpha}}{\partial x’^{\mu}}\frac{\partial x^{\beta}}{\partial x’^{\nu}}\Gamma^{\lambda}_{\ \alpha\beta}+\frac{\partial x’^{\rho}}{\partial x^{\lambda}}\frac{\partial^{2}x^{\lambda}}{\partial x’^{\mu}\partial x’^{\nu}},$$ where the first term is what a tensor would do and the second is the repair. The second derivatives of the coordinate change appear here with exactly the structure needed to cancel the ones in the obstruction. So a connection is not a tensor, and an equation saying $\Gamma=0$ in some coordinates carries no invariant information at all. That will matter quite a lot in a moment.
Nothing so far says which connection to take, and in fact any set of coefficients with the right transformation behaviour defines a legitimate covariant derivative. General relativity makes two additional demands, and it’s worth stating them separately because they’re independent and neither is a theorem. The first is metric compatibility, $$\nabla_{\rho}g_{\mu\nu}=0,$$ which says that transporting two vectors along a curve preserves their inner product, so lengths and angles are carried unchanged. The second is vanishing torsion, $$\Gamma^{\rho}_{\ \mu\nu}=\Gamma^{\rho}_{\ \nu\mu},$$ symmetry in the two lower indices, which says that the failure of an infinitesimal parallelogram to close is not part of the structure. Metric compatibility doesn’t imply the symmetry, and theories keeping one without the other are perfectly consistent.4Einstein-Cartan theory is the standard example, where torsion is sourced by intrinsic spin and vanishes outside matter, so the vacuum predictions coincide with general relativity. These are choices that define standard general relativity, not mathematical necessities, and together they pin the connection down completely.
The derivation is short. Expanding metric compatibility with one $\Gamma$ per index gives $\partial_{\rho}g_{\mu\nu}=\Gamma_{\nu\rho\mu}+\Gamma_{\mu\rho\nu}$, where $\Gamma_{\alpha\beta\gamma}\equiv g_{\alpha\lambda}\Gamma^{\lambda}_{\ \beta\gamma}$. Write it three times with the indices permuted, $$\partial_{\mu}g_{\nu\rho}=\Gamma_{\rho\mu\nu}+\Gamma_{\nu\mu\rho},\qquad \partial_{\nu}g_{\rho\mu}=\Gamma_{\mu\nu\rho}+\Gamma_{\rho\nu\mu},\qquad \partial_{\rho}g_{\mu\nu}=\Gamma_{\nu\rho\mu}+\Gamma_{\mu\rho\nu},$$ then add the first two and subtract the third. Symmetry in the last two indices makes $\Gamma_{\nu\mu\rho}$ cancel against $\Gamma_{\nu\rho\mu}$ and $\Gamma_{\mu\nu\rho}$ against $\Gamma_{\mu\rho\nu}$, leaving $2\Gamma_{\rho\mu\nu}$. Raising the first index, $$\Gamma^{\rho}_{\ \mu\nu}=\tfrac12 g^{\rho\sigma}\left(\partial_{\mu}g_{\sigma\nu}+\partial_{\nu}g_{\sigma\mu}-\partial_{\sigma}g_{\mu\nu}\right),$$ the Christoffel symbols of the Levi-Civita connection. Every term is a first derivative of the metric, which is the fact we’ll use twice over in what follows.
It’s worth stopping to see what has and hasn’t been achieved. We set out to find the part of gravity that free fall can’t remove, and we haven’t found it. What we’ve built is the machinery for comparing directions at different events, which we needed because without it we couldn’t write down a single coordinate-independent statement about how anything varies from place to place. The connection is that machinery, and by construction it isn’t a tensor, so it can’t itself be the invariant we’re after. What it does give us is a way of saying what free motion means, and that’s the next thing to check.
A particle moves freely if its velocity is carried along unchanged by the transport the geometry provides, which is to say $$\nabla_{u}u=0,\qquad u^{\mu}=\frac{dx^{\mu}}{d\tau},\qquad g_{\mu\nu}u^{\mu}u^{\nu}=-c^{2}.$$ Writing it in components, $(\nabla_{u}u)^{\mu}=u^{\nu}\nabla_{\nu}u^{\mu}=u^{\nu}(\partial_{\nu}u^{\mu}+\Gamma^{\mu}_{\ \nu\lambda}u^{\lambda})$. The first piece simplifies because of what $u^{\nu}\partial_{\nu}$ means. The four-velocity is tangent to the worldline, so contracting it with a gradient is the directional derivative along that curve, and for any quantity defined along the worldline the chain rule gives $u^{\nu}\partial_{\nu}f=\frac{dx^{\nu}}{d\tau}\frac{\partial f}{\partial x^{\nu}}=\frac{df}{d\tau}$. Applied to $u^{\mu}$ this turns $u^{\nu}\partial_{\nu}u^{\mu}$ into $du^{\mu}/d\tau=d^{2}x^{\mu}/d\tau^{2}$, and the whole condition becomes $$\frac{d^{2}x^{\mu}}{d\tau^{2}}+\Gamma^{\mu}_{\ \alpha\beta}\frac{dx^{\alpha}}{d\tau}\frac{dx^{\beta}}{d\tau}=0,$$ the geodesic equation. A curve satisfying it parallel-transports its own tangent vector, which is the closest available statement of going straight, and when the metric components are constant every $\Gamma$ vanishes and it reduces to $d^{2}x^{\mu}/d\tau^{2}=0$.
The same equation follows from a completely different demand, and the agreement is worth noting because it isn’t obvious. In special relativity, among timelike paths joining two events the straight one has the largest proper time; that’s the twin paradox as a variational statement. Extremising $\tau=\frac{1}{c}\int\sqrt{-g_{\mu\nu}dx^{\mu}dx^{\nu}}$ is awkward because of the square root. It’s convenient to use an affine parameter $\lambda$ along the path and extremise instead $$S=\frac12\int g_{\mu\nu}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda}\,d\lambda,$$ the same trick as using $\tfrac12 mv^{2}$ rather than a path length in ordinary mechanics. This gives the same curves, and for a timelike geodesic we may afterwards take $\lambda=\tau$, up to an affine rescaling. Run the Euler-Lagrange equations on that integrand and the derivatives of the metric assemble themselves, without being asked, into precisely the combination that defines $\Gamma^{\sigma}_{\ \mu\nu}$, and the geodesic equation comes back out.5The steps, for anyone checking: with $L=\tfrac12 g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$ and a dot meaning $d/d\lambda$, we have $\partial L/\partial\dot{x}^{\rho}=g_{\rho\nu}\dot{x}^{\nu}$, whose derivative along the curve is $\partial_{\mu}g_{\rho\nu}\dot{x}^{\mu}\dot{x}^{\nu}+g_{\rho\nu}\ddot{x}^{\nu}$, set against $\partial L/\partial x^{\rho}=\tfrac12\partial_{\rho}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$. The first term is contracted with something symmetric in $\mu\nu$ and may be symmetrised, and multiplying through by $g^{\sigma\rho}$ produces the Christoffel combination. So a geodesic has two independent characterisations: it parallel-transports its own tangent, and it extremises proper time. That they agree is not automatic. The first notion needs only a connection and the second needs only a metric, and they coincide because we chose the connection to be metric-compatible and torsion-free; for a general connection the autoparallels and the extremal-proper-time curves are different families.6Timelike geodesics extremise proper time locally and are genuinely maximal for sufficiently nearby events, but over long paths families of geodesics can refocus at conjugate points, beyond which a geodesic is no longer the longest path between its endpoints. I say extremise rather than maximise throughout for that reason.
Now go back to the freely falling frame and put the pieces together. In those coordinates $g_{\mu\nu}(P)=\eta_{\mu\nu}$ and $\partial_{\alpha}g_{\mu\nu}(P)=0$, and every term in the Levi-Civita formula is a first derivative of the metric, so $$\Gamma^{\rho}_{\ \mu\nu}(P)=\tfrac12 g^{\rho\sigma}(0+0-0)=0.$$ This is the precise statement of what free fall accomplishes, and it’s worth spelling out the chain rather than compressing it. The first derivatives of the metric determine the connection; the connection coefficients are what appear in the geodesic equation; and there they supply the entire coordinate acceleration of a freely moving particle, since $d^{2}x^{\mu}/d\tau^{2}=-\Gamma^{\mu}_{\ \alpha\beta}\dot{x}^{\alpha}\dot{x}^{\beta}$. Kill the $\Gamma$ at $P$ and a free particle has no coordinate acceleration at $P$, which is the occupant of the elevator watching the ball hang in the air. The first derivatives of the metric aren’t themselves the gravitational field; they’re the ingredients of the object that shows up as acceleration once you ask how a free particle moves.
At this point everything appears to line up. The geodesic equation rearranges to acceleration equals minus a $\Gamma$ term, sitting in exactly the slot where $-\nabla\Phi$ sits in Newton’s law, and the connection coefficients are built from first derivatives of the metric, just as the Newtonian field is the first derivative of the potential. On that basis you’ll see the $\Gamma$ called the gravitational field, and it’s hard to see what else they could be.
So compute them for a space with no gravity in it at all. Take the flat tabletop from earlier in polar coordinates, $$dl^{2}=dr^{2}+r^{2}d\theta^{2},$$ so $g_{rr}=1$, $g_{\theta\theta}=r^{2}$, and the inverse components are $g^{rr}=1$, $g^{\theta\theta}=1/r^{2}$. The only nonvanishing derivative is $\partial_{r}g_{\theta\theta}=2r$, and feeding that into the Levi-Civita formula gives $$\Gamma^{r}_{\ \theta\theta}=-\tfrac12 g^{rr}\partial_{r}g_{\theta\theta}=-r,\qquad \Gamma^{\theta}_{\ r\theta}=\Gamma^{\theta}_{\ \theta r}=\tfrac12 g^{\theta\theta}\partial_{r}g_{\theta\theta}=\frac1r,$$ with every other component zero. The space is a plane. There’s no gravitational field anywhere, no matter anywhere, nothing curved about it, and the connection coefficients are not zero. The geodesic equation becomes $\ddot{r}-r\dot{\theta}^{2}=0$ together with an equation for $\theta$, whose solutions are straight lines written in polar form, and that $-r\dot{\theta}^{2}$ is the centrifugal term from first-year mechanics, appearing here as a connection coefficient rather than a force. Nonzero Christoffel symbols do not imply curvature. They record that the coordinate basis rotates and stretches from point to point, which is a fact about the grid, and the warning attached to this same line element earlier was exactly this.
The same thing happens in spacetime. An observer with constant proper acceleration in flat Minkowski spacetime can use coordinates adapted to their own motion, called Rindler coordinates, in which the metric components aren’t $\eta_{\mu\nu}$, the connection coefficients are nonzero, and objects released from rest appear to accelerate away. Locally that observer could reasonably call what they feel a gravitational field, while the spacetime is exactly flat. Put that beside the calculation two paragraphs above and both directions are now established. In flat spacetime we can have $\Gamma^{\rho}_{\ \mu\nu}\neq0$, and in curved spacetime we can have $\Gamma^{\rho}_{\ \mu\nu}(P)=0$ at any event we like. A quantity that can be made nonzero where there’s no curvature and zero where there is curvature cannot be the invariant measure of curvature, and this is not an argument about intuition; it follows from the inhomogeneous transformation law, since an object with a non-tensorial piece can always be shifted by a suitable relabelling. Whatever survives free fall isn’t the connection.
Which brings back the question left standing earlier: how far can a change of coordinates go? Counting the freedom gives a rough sense of where it runs out, and the count is supporting intuition rather than proof, so it’s worth doing quickly and not leaning on it too hard. A coordinate change near $P$ is specified by its Taylor coefficients: the first derivatives $\partial x’^{\mu}/\partial x^{\alpha}$ form a $4\times4$ array with $16$ entries; the second derivatives are symmetric in their two lower indices, giving $10$ pairs for each of $4$ values of $\mu$, so $40$; the third derivatives are symmetric in three lower indices, giving $20$ triples for each of $4$, so $80$. What we’d like removed is the metric and its derivatives: $g_{\mu\nu}$ has $10$ independent components, $\partial_{\alpha}g_{\mu\nu}$ has $40$, and $\partial_{\alpha}\partial_{\beta}g_{\mu\nu}$, symmetric in $\alpha\beta$ and in $\mu\nu$ separately, has $100$. Setting $g_{\mu\nu}=\eta_{\mu\nu}$ imposes $10$ conditions on $16$ parameters and leaves $6$ over, which are the dimension of the Lorentz group and the residual freedom to boost and rotate a local inertial frame that special relativity says should be there. Setting $\partial_{\alpha}g_{\mu\nu}=0$ imposes $40$ conditions with exactly $40$ parameters. That matching count is consistent with the existence of locally inertial coordinates, although the count itself is not a proof of it. Setting $\partial_{\alpha}\partial_{\beta}g_{\mu\nu}=0$ would impose $100$ conditions with only $80$ parameters, and it fails, short by $20$.
A suggestive count of the coordinate freedom available at an event. A symmetric pair of four-dimensional indices has $10$ independent choices, while a symmetric triple has $20$. The count is consistent with the fact that the metric and its first derivatives can be put into locally inertial form, while second-order information remains. The shortfall of $20$ matches the number of independent components of the Riemann tensor in four dimensions, although this counting argument is not itself a proof.
Take that as a hint about where to look rather than as a result. The tally treats coordinate changes as though they acted on the metric derivatives one parameter at a time, which they don’t, and it proves nothing about what the leftover information is. Two things it does suggest. First, the coordinate freedom available to remove the metric’s first derivatives is just sufficient, while at second derivatives something is left that no relabelling reaches, which is the same order at which the Newtonian residue lived. Second, the size of the shortfall is consistent with the fact that the Riemann tensor in four dimensions has twenty independent components, which is a remark worth filing and not an argument for its existence or structure. The distinction to keep hold of is that the raw second derivatives $\partial_{\alpha}\partial_{\beta}g_{\mu\nu}$ are themselves coordinate-dependent, and individual components can be changed at will by a different choice of labels. What can’t be changed is a particular combination of them together with products of first derivatives, and finding that combination is a separate job the count can’t do for us.
Each small neighbourhood can be described by an approximately flat local frame, but curvature prevents those frames from fitting together into one globally flat description. The surface is drawn embedded only for illustration; the curvature being discussed is intrinsic.
One picture is worth clearing away before we start. The bowling ball on a stretched rubber sheet explains the marble’s motion by appealing to real gravity pulling it into the dip, so it uses gravity to explain gravity; it shows two spatial dimensions, so what it depicts is spatial curvature rather than the curvature of spacetime; and it hides the component that actually matters, since we’ll find in the weak-field limit that what makes a slowly moving body fall comes from $g_{00}$, which is a feature of the temporal part of the geometry that no curved spatial surface can display. We’ll work intrinsically instead, with measurements made inside the geometry and no reference to a higher-dimensional space for it to sit in.
The intrinsic notion we need is already in hand, because a connection is exactly a rule for carrying a vector from one tangent space to the next. Carry one along a curve, keeping it parallel by that rule at every step, and you have parallel transport. Now the question that matters: if two different curves run from $P$ to $Q$, do they deliver the same vector at $Q$? Equivalently, transport a vector around a closed loop and ask whether it returns to itself. In a simply connected flat region with the ordinary rule it does. On a sphere it doesn’t. Stand on the equator with an arrow pointing north, carry it a quarter of the way around the equator keeping it parallel to itself, then up a meridian to the pole, then back down the next meridian to where you started, and the arrow comes home rotated by a right angle. Nothing twisted it. The rotation is measurable from within the surface, with no reference to any space the sphere might sit in, and its size depends on the loop: for a sphere of radius $R$ the angle equals the enclosed area divided by $R^{2}$, so shrinking the loop shrinks the effect in proportion to the area. This path-dependence is holonomy. For an infinitesimal loop, its leading behaviour is controlled locally by the curvature tensor, which is the object we’re about to build, and that’s the route we’ll take: holonomy is built from the connection, which isn’t a tensor, and yet it’s stated as a comparison between two transports, which is a question with an observer-independent answer. It also explains the failure of the coordinate count in geometric terms: each small patch carries a good flat frame, and carrying that frame around a loop doesn’t return it to itself, so the patches can’t be assembled into one flat frame over a finite region.
Parallel transport around the closed path A → B → C → A returns the vector to its starting point rotated by $\pi/2$. For a sphere, the rotation equals the enclosed area divided by the square of the sphere’s radius.
To get a formula we make the loop infinitesimal. Transport around a small rectangle with sides along the $\mu$ and $\nu$ directions means going out one way and back the other, and comparing the two routes is the same as differentiating first in one direction and then the other and asking whether the order matters. So compute the commutator of two covariant derivatives. Write $W_{\nu}{}^{\rho}\equiv\nabla_{\nu}V^{\rho}$, which has one upper and one lower index, so differentiating it covariantly needs one $\Gamma$ for each: $$\nabla_{\mu}\nabla_{\nu}V^{\rho}=\partial_{\mu}W_{\nu}{}^{\rho}+\Gamma^{\rho}_{\ \mu\lambda}W_{\nu}{}^{\lambda}-\Gamma^{\lambda}_{\ \mu\nu}W_{\lambda}{}^{\rho}.$$ Substituting $W_{\nu}{}^{\rho}=\partial_{\nu}V^{\rho}+\Gamma^{\rho}_{\ \nu\lambda}V^{\lambda}$ and expanding, $$\nabla_{\mu}\nabla_{\nu}V^{\rho}=\partial_{\mu}\partial_{\nu}V^{\rho}+(\partial_{\mu}\Gamma^{\rho}_{\ \nu\lambda})V^{\lambda}+\Gamma^{\rho}_{\ \nu\lambda}\partial_{\mu}V^{\lambda}+\Gamma^{\rho}_{\ \mu\lambda}\partial_{\nu}V^{\lambda}+\Gamma^{\rho}_{\ \mu\lambda}\Gamma^{\lambda}_{\ \nu\sigma}V^{\sigma}-\Gamma^{\lambda}_{\ \mu\nu}\left(\partial_{\lambda}V^{\rho}+\Gamma^{\rho}_{\ \lambda\sigma}V^{\sigma}\right).$$ Now write the same thing with $\mu$ and $\nu$ exchanged and subtract, term by term. The second derivatives go first: $\partial_{\mu}\partial_{\nu}V^{\rho}-\partial_{\nu}\partial_{\mu}V^{\rho}=0$, since partial derivatives commute. The terms carrying a first derivative of $V$ go next: the pair $\Gamma^{\rho}_{\ \nu\lambda}\partial_{\mu}V^{\lambda}+\Gamma^{\rho}_{\ \mu\lambda}\partial_{\nu}V^{\lambda}$ is already symmetric under exchanging $\mu$ and $\nu$, so it cancels against itself. The last bracket cancels too, because $\Gamma^{\lambda}_{\ \mu\nu}$ is symmetric in its lower indices, which is where the torsion-free assumption earns its keep.7With torsion the antisymmetric part of $\Gamma^{\lambda}_{\ \mu\nu}$ survives here and contributes an extra term $-T^{\lambda}{}_{\mu\nu}\nabla_{\lambda}V^{\rho}$ to the commutator. Dropping it is a consequence of the choice made earlier, not an algebraic accident. Everything containing a derivative of $V$ has now gone, which is the important structural fact: what’s left acts on $V$ itself, with no derivatives, so it’s a linear map on vectors rather than a differential operator. What’s left is $$\left(\nabla_{\mu}\nabla_{\nu}-\nabla_{\nu}\nabla_{\mu}\right)V^{\rho}=\left(\partial_{\mu}\Gamma^{\rho}_{\ \nu\sigma}-\partial_{\nu}\Gamma^{\rho}_{\ \mu\sigma}+\Gamma^{\rho}_{\ \mu\lambda}\Gamma^{\lambda}_{\ \nu\sigma}-\Gamma^{\rho}_{\ \nu\lambda}\Gamma^{\lambda}_{\ \mu\sigma}\right)V^{\sigma},$$ and we give the bracket a name, $$R^{\rho}_{\ \sigma\mu\nu}=\partial_{\mu}\Gamma^{\rho}_{\ \nu\sigma}-\partial_{\nu}\Gamma^{\rho}_{\ \mu\sigma}+\Gamma^{\rho}_{\ \mu\lambda}\Gamma^{\lambda}_{\ \nu\sigma}-\Gamma^{\rho}_{\ \nu\lambda}\Gamma^{\lambda}_{\ \mu\sigma},$$ the Riemann curvature tensor. It’s worth being clear that this is not a formula someone wrote down and then found a use for. It’s the coefficient left behind when covariant derivatives fail to commute, and the expression is whatever the cancellation leaves.
The indices now say something definite. The left side is linear in $V$, so $\sigma$ labels the component of the vector fed in and $\rho$ the component of the vector that comes back out, exactly as for any linear map on the tangent space. The pair $\mu\nu$ is antisymmetric, since exchanging them changes the sign of the commutator, and an antisymmetric pair of directions is an oriented two-plane: it specifies the infinitesimal patch the transport was taken around and which way round it was traversed. So for each oriented two-plane at an event, $R^{\rho}_{\ \sigma\mu\nu}$ is a linear map taking a vector to the change it suffers when carried around that patch. Four indices is the least this can be done with, and it also explains why the object is a tensor even though the $\Gamma$ are not. The left side is manifestly a tensor: it’s built from covariant derivatives, each of which was constructed to transform correctly, and the inhomogeneous pieces of the $\Gamma$ cancel there rather than surviving. Since the left side is a tensor and $V^{\sigma}$ is arbitrary, the coefficient has to be one too. This is the identification we were promised, although the usual shorthand for it needs one qualification. People say the curvature is a combination of second derivatives of the metric, and it’s worth seeing why that’s only half the story. Each $\Gamma$ already carries one derivative of the metric, so the definition above is schematically $$R\sim\partial^{2}g+(\partial g)^{2},$$ with both kinds of term present. The quadratic piece isn’t decoration; it’s part of what makes the non-tensorial contributions cancel, and an object built from $\partial^{2}g$ alone wouldn’t be a tensor. The two terms separate only in a locally inertial frame at $P$, where $\Gamma(P)=0$ kills the quadratic part at that one event and leaves the curvature there expressed through second derivatives of the metric alone. With that understood, the Riemann tensor is invariant where the raw second derivatives weren’t, and being a tensor, nonzero in one coordinate system means nonzero in all of them.
With it in hand we can say how local a local inertial frame really is. In the locally inertial coordinates at $P$, usually called Riemann normal coordinates, the metric has the expansion $$g_{\mu\nu}(x)=\eta_{\mu\nu}-\tfrac13 R_{\mu\alpha\nu\beta}(P)\,x^{\alpha}x^{\beta}+O(x^{3}).$$ Coordinates of this kind can be built around any event of any spacetime; that’s a standard result which I’ll take on trust rather than prove, since what we need is not the construction but what the expansion says. Read it carefully. The constant term is the first condition free fall achieved and the absence of a linear term is the second, and the content of the expansion is that in such coordinates the first invariant information in the metric appears at quadratic order and is controlled by the curvature. That’s not the same as saying the leftover second derivatives are the curvature components; it’s saying that once the coordinate freedom has been spent, what remains at that order is governed by a tensor. Practically, a laboratory of size $L$ around $P$ departs from flat geometry by a fractional amount of order $\mathcal{R}L^{2}$, where $\mathcal{R}$ stands for the characteristic size of the relevant components of the Riemann tensor rather than the Ricci scalar or any other particular contraction, and that’s the number to hold against your experimental precision. It’s worth keeping this apart from the estimate made by hand much earlier, because the two measure different things. That one was mechanical, a tidal relative acceleration compared with $g$, and gave a few parts in $10^{7}$ across a two-metre laboratory. This one is a dimensionless departure of the metric from Minkowski form, and near the Earth the characteristic curvature is $\mathcal{R}\sim GM_\oplus/c^{2}R_\oplus^{3}\sim10^{-23}\ \mathrm{m^{-2}}$, so $\mathcal{R}L^{2}\sim10^{-22}$ for the same laboratory. The gap between the two numbers is a factor of order $gL/c^{2}$, which is what separates a statement about accelerations from a statement about geometry.
The polar example settles immediately, which is a useful check on the machinery. Substituting $\Gamma^{r}_{\ \theta\theta}=-r$ and $\Gamma^{\theta}_{\ r\theta}=1/r$ into the definition, the derivative terms cancel against the quadratic terms and every component vanishes, so $R^{\rho}_{\ \sigma\mu\nu}=0$ for a plane in any coordinates whatsoever. The connection was a fact about the grid; the curvature is a fact about the surface.
Quantity
What it tells us
What it means
$\Gamma^{\rho}{}_{\mu\nu}\neq 0$
The coordinate basis changes from point to point.
Not necessarily curvature. Christoffel symbols can be nonzero even in flat spacetime, as in polar coordinates.
$R^{\rho}{}_{\sigma\mu\nu}\neq 0$
Parallel transport depends on the path.
Intrinsic curvature. No coordinate transformation can make the Riemann tensor vanish at that event.
Nearby freely falling worldlines accelerate relative to one another.
An observable tidal effect. For a chosen four-velocity and separation, curvature produces geodesic deviation.
These statements are related but not equivalent. Connection coefficients describe how the chosen coordinate basis changes; the Riemann tensor describes intrinsic curvature; geodesic deviation is one operational way that curvature is measured using neighbouring freely falling bodies.
What remains is to connect this tensor to the thing we started from, the relative acceleration of two nearby falling bodies, and the derivation is short enough to do properly. Take a smooth one-parameter family of geodesics and label them by $s$, so that $x^{\mu}(\tau,s)$ is the point reached at proper time $\tau$ along the geodesic numbered $s$. The family sweeps out a two-dimensional surface in spacetime on which $\tau$ and $s$ serve as coordinates, and it carries two natural vector fields, $$u^{\mu}=\frac{\partial x^{\mu}}{\partial\tau},\qquad \xi^{\mu}=\frac{\partial x^{\mu}}{\partial s}.$$ The first is the four-velocity, tangent to each geodesic. The second points from a given geodesic to its neighbour at the same value of $\tau$, and it’s the mathematical version of the vector between the two balls in the elevator. Because $\tau$ and $s$ are coordinates on that surface, $u$ and $\xi$ are the basis fields belonging to those two directions, and mixed partial derivatives commute. Said without notation: move a little way along one geodesic and then step across to its neighbour, or step across first and move along afterwards, and you arrive at the same event either way. The object measuring the failure of two such routes to agree is called the Lie bracket $[u,\xi]$, and what we’ve just observed is that here it vanishes, $[u,\xi]=0$. For a torsion-free connection the bracket can be written with covariant derivatives, $[u,\xi]^{\mu}=u^{\nu}\nabla_{\nu}\xi^{\mu}-\xi^{\nu}\nabla_{\nu}u^{\mu}$, since the $\Gamma$ terms cancel by symmetry in the lower indices, and therefore $$\nabla_{u}\xi=\nabla_{\xi}u.$$ That’s the one nontrivial input, and it’s again the torsion-free choice doing work.
Now differentiate the separation twice along the worldline. Using the relation just obtained, $$\frac{D^{2}\xi^{\mu}}{D\tau^{2}}=\nabla_{u}\nabla_{u}\xi^{\mu}=\nabla_{u}\nabla_{\xi}u^{\mu}.$$ To evaluate the right side we swap the order of the two derivatives, which is precisely what the curvature tensor was built to account for, $$\nabla_{u}\nabla_{\xi}u^{\mu}=\nabla_{\xi}\nabla_{u}u^{\mu}+u^{\alpha}\xi^{\beta}\left(\nabla_{\alpha}\nabla_{\beta}-\nabla_{\beta}\nabla_{\alpha}\right)u^{\mu},$$ where the bracket term is the commutator acting on $u^{\mu}$ and there’s no extra contribution because $[u,\xi]=0$. The first term on the right vanishes, since each curve of the family is a geodesic and $\nabla_{u}u=0$. The second is the commutator identity, which gives $u^{\alpha}\xi^{\beta}R^{\mu}_{\ \nu\alpha\beta}u^{\nu}$. So $$\frac{D^{2}\xi^{\mu}}{D\tau^{2}}=R^{\mu}_{\ \nu\alpha\beta}\,u^{\nu}u^{\alpha}\xi^{\beta}.$$ One step is left, and it’s the one that’s usually skipped. The Riemann tensor is antisymmetric in its last two indices, which we know because the commutator changes sign when the two derivatives are exchanged, so $R^{\mu}_{\ \nu\alpha\beta}u^{\alpha}\xi^{\beta}=-R^{\mu}_{\ \nu\beta\alpha}u^{\alpha}\xi^{\beta}$. Relabelling the dummy indices, $\beta\to\rho$ and $\alpha\to\sigma$, this is $-R^{\mu}_{\ \nu\rho\sigma}\xi^{\rho}u^{\sigma}$, and the equation takes its standard form $$\frac{D^{2}\xi^{\mu}}{D\tau^{2}}=-R^{\mu}_{\ \nu\rho\sigma}\,u^{\nu}\xi^{\rho}u^{\sigma},$$ the equation of geodesic deviation, with signs fixed by the curvature convention chosen above.8Sign conventions for the Riemann tensor differ between textbooks, and the overall sign of this equation differs with them. If you compare against another source and the signs disagree, check its definition of $R^{\rho}_{\ \sigma\mu\nu}$ and its metric signature before concluding that one of the two is wrong. Every object in it is a tensor, so the statement holds in every coordinate system: the components of $\xi^{\mu}$ and of the relative acceleration are as observer-dependent as the components of any vector, and yet if the right side is nonzero for one observer it’s nonzero for all of them. That is the invariant we went looking for, and it took this long to state because nothing smaller would do.
Now compare it with the Newtonian tidal equation, and do it carefully, because the correspondence is more exact than a slogan. Fix an observer with four-velocity $u^{\mu}$. Then the geodesic deviation equation defines a map that takes the separation and returns the relative acceleration, $$\xi^{\rho}\;\longmapsto\;-R^{\mu}_{\ \nu\rho\sigma}u^{\nu}u^{\sigma}\xi^{\rho},$$ which is linear in $\xi$ and depends on $u$ but not on any property of the two falling bodies. This is the relativistic tidal operator measured along that worldline. Newtonian gravity has the same structure, $$\xi^{j}\;\longmapsto\;-\partial_{i}\partial_{j}\Phi\,\xi^{j},$$ a linear map from separation to relative acceleration built from second derivatives of the potential. The input and output roles match directly: $\rho$ and $j$ label the direction of separation, while $\mu$ and $i$ label the direction of the resulting relative acceleration. The difference is the pair of velocity contractions, which the relativistic version needs because the answer depends on how the pair is moving through spacetime, and which the Newtonian version could do without because it assumed a universal time in which everyone advances into the future at the same rate. So the right statement isn’t that the Riemann tensor is the tidal matrix. It’s that contracting the Riemann tensor twice with a given four-velocity produces the tidal operator for observers with that velocity, and the Newtonian Hessian is what that operator reduces to when the velocity is small and the field is weak. The dictionary, which we’ll now verify, is $R^{i}_{\ 0j0}\leftrightarrow\partial_{i}\partial_{j}\Phi/c^{2}$. Notice also that a single worldline appears nowhere on the right and couldn’t, since the left side needs two of them to be defined at all, which is the earlier remark about pointlike frames, now as an equation.
A small cloud of neighbouring freely falling particles is initially spherical but becomes stretched in the radial direction and compressed transversely. This relative deformation is the observable content of tidal gravity and, in general relativity, is described by geodesic deviation.
The correspondence can be made quantitative, and it’s worth stating the assumptions before starting, since all four are doing work. The field is weak, so $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ with $|h_{\mu\nu}|\ll1$ and we keep only first order in $h$. The field is static, so $\partial_{0}g_{\mu\nu}=0$. The test particle is slow, so $|dx^{i}/dt|\ll c$ and $d\tau\simeq dt$. And the signature is $(-,+,+,+)$, which fixes the signs throughout. In the geodesic equation the dominant term in $\Gamma^{\mu}_{\ \alpha\beta}\dot{x}^{\alpha}\dot{x}^{\beta}$ is the one with both indices zero, since $\dot{x}^{0}=c\,dt/d\tau\simeq c$ swamps the spatial components, leaving $d^{2}x^{i}/d\tau^{2}\simeq-c^{2}\Gamma^{i}_{\ 00}$. From the Levi-Civita formula, $\Gamma^{i}_{\ 00}=\tfrac12 g^{i\sigma}(2\partial_{0}g_{\sigma0}-\partial_{\sigma}g_{00})$, and staticity kills the first part while $g^{ij}\simeq\delta^{ij}$ to the order we’re working, so $$\Gamma^{i}_{\ 00}\simeq-\tfrac12\delta^{ij}\partial_{j}g_{00},\qquad \frac{d^{2}x^{i}}{dt^{2}}\simeq\frac{c^{2}}{2}\partial_{i}g_{00}.$$ Newton says $d^{2}x^{i}/dt^{2}=-\partial_{i}\Phi$, so the two agree if $\tfrac12 c^{2}\partial_{i}g_{00}=-\partial_{i}\Phi$, which integrates to $$g_{00}\simeq-\left(1+\frac{2\Phi}{c^{2}}\right)$$ with the constant of integration fixed by demanding $g_{00}\to-1$ far from any mass, where the geometry must reduce to Minkowski. Substituting back, $\Gamma^{i}_{\ 00}\simeq+\partial_{i}\Phi/c^{2}$, and $d^{2}x^{i}/dt^{2}\simeq-c^{2}\Gamma^{i}_{\ 00}=-\partial_{i}\Phi$, which is the check that the signs are consistent. The Newtonian potential isn’t a separate field sitting alongside the geometry. In this limit it’s the leading correction to the time-time component of the metric, which is also the honest replacement for the rubber sheet, since what makes an apple fall is a gradient in $g_{00}$ and no picture of a curved spatial surface can display that. And because $g_{00}$ is the component that fixes how much proper time a stationary clock accumulates, the potential turns up in exactly the part of the weak-field metric that controls gravitational time dilation: relative to infinity, the fractional clock-rate shift is $\Phi/c^{2}$. Since $\Phi<0$, a clock deeper in the well runs slower, by about $7\times10^{-10}$ at the Earth's surface, which the GPS system has to correct for.
A stationary clock deeper in a gravitational potential accumulates less proper time than one farther from the source. In the weak-field limit this follows from $g_{00}\simeq-(1+2\Phi/c^2)$.
The tidal correspondence follows from the same expansion. For a static weak field we can drop time derivatives and terms quadratic in $\Gamma$, since each $\Gamma$ is already first order in $h$, so the curvature reduces to $R^{i}_{\ 0j0}\simeq\partial_{j}\Gamma^{i}_{\ 00}$, and with $\Gamma^{i}_{\ 00}\simeq\partial_{i}\Phi/c^{2}$ from above this gives $$R^{i}_{\ 0j0}\simeq\frac{1}{c^{2}}\partial_{i}\partial_{j}\Phi,$$ which is the dictionary entry promised a moment ago. Feed it into the geodesic deviation equation with $u^{\mu}\simeq(c,0,0,0)$, so that only the $\nu=\sigma=0$ terms survive, and $$\frac{d^{2}\xi^{i}}{dt^{2}}\simeq-R^{i}_{\ 0j0}\,c^{2}\,\xi^{j}=-\partial_{i}\partial_{j}\Phi\,\xi^{j},$$ the Newtonian tidal equation we wrote down before any geometry had appeared. The two factors of $c$ from the four-velocity cancel the $1/c^{2}$ in the curvature, reproducing a finite Newtonian tidal acceleration in the nonrelativistic limit. For the Earth, with $\Phi=-GM/r$, this reproduces stretching along the radius at rate $2GM/r^{3}$ and squeezing in each transverse direction at $GM/r^{3}$, and the trace vanishes in vacuum, so a small falling cloud holds its volume to leading order while its shape distorts. The loop is closed: the quantity the Newtonian argument identified as the irremovable part of gravity is a set of components of the Riemann tensor, and the relativistic tidal operator reduces to the Newtonian Hessian in the limit where Newton’s assumptions hold.
So it’s now possible to say what gravity is in this framework, in terms that refer only to things one could measure. Spacetime carries a metric, which fixes proper times, distances, orthogonality and causal structure, and which determines a unique connection that’s compatible with it and free of torsion. That connection defines parallel transport, and free particles are those whose velocity is parallel-transported along their own worldline. At any event a freely falling observer can choose coordinates in which the metric components are Minkowskian and the connection coefficients vanish, so the non-gravitational laws take their special-relativistic form and a free particle shows no coordinate acceleration, and such a frame can be maintained all along that observer’s worldline. What can’t be done is extending it to a finite region, and the obstruction is the Riemann tensor, the particular combination of first and second derivatives of the metric that no coordinate change removes. Its observable signature is the relative acceleration of neighbouring freely falling bodies through the geodesic deviation equation, and in the weak-field limit the relevant components reduce to the Hessian of the Newtonian potential. It would be too strong to say that only the curvature is really there. The metric is the fundamental field of the theory, and it’s what everything else is built from; what the Riemann tensor provides is the coordinate-independent local gravitational content, the part of the description that survives every relabelling and every change of observer.
That’s as far as the equivalence principle takes us, and it’s worth being exact about where it stops. It told us to look for local inertial frames, and the failure of those frames to fit together forced a metric, then a connection, then covariant differentiation, then curvature, each one because the mathematics we had couldn’t express what the physics required. What we have at the end is a complete account of how free motion works in a given geometry, and a way of detecting the invariant curvature of that geometry by watching neighbouring free-fall trajectories. What we don’t have is any reason why spacetime should carry one geometry rather than another. Nothing in the argument so far mentions the Earth. The metric is still an unknown function, and we’ve never explained why a particular distribution of matter produces a particular metric, which means we can recognise gravity but can’t yet predict it.
That’s as far as the equivalence principle takes us, and it’s worth being exact about where it stops. It told us to look for local inertial frames, and the failure of those frames to fit together forced a metric, then a connection, then covariant differentiation, then curvature, each one because the mathematics we had couldn’t express what the physics required. What we have at the end is a complete account of how free motion works in a given geometry, and a way of detecting the invariant curvature of that geometry by watching neighbouring free-fall trajectories. What we don’t have is any reason why spacetime should carry one geometry rather than another. Nothing in the argument so far mentions the Earth. The metric is still an unknown function, and we’ve never explained why a particular distribution of matter produces a particular metric, which means we can recognise gravity but can’t yet predict it. Answering that needs an equation with matter on one side and geometry on the other: the stress-energy tensor $T_{\mu\nu}$ to describe the matter, and on the other side something built from the curvature, which turns out to be constrained rather tightly. Contracting the Riemann tensor gives the Ricci tensor $R_{\mu\nu}=R^{\alpha}{}_{\mu\alpha\nu}$ and contracting again gives the Ricci scalar $R=g^{\mu\nu}R_{\mu\nu}$, and the combination that matters,
is singled out because the Bianchi identity forces $\nabla^{\mu}G_{\mu\nu}=0$ identically, matching the conservation law $\nabla^{\mu}T_{\mu\nu}=0$ that any sensible matter has to satisfy. Setting the two proportional,
$$
G_{\mu\nu}=\frac{8\pi G}{c^{4}}T_{\mu\nu},
$$
is the Einstein field equation, and that’s where this goes next.
References and Footnotes
1
This is a reconstruction of the logic rather than the history. Einstein worked at the problem for most of a decade with several wrong turns, and the question here is whether we, knowing modern physics, can reconstruct why this structure makes sense. ↩︎
2
Coordinates of this kind are called Fermi normal coordinates. They put the metric in Minkowski form and set the connection coefficients to zero at every point of the central worldline, with corrections appearing only at second order in the distance away from it. Both of those objects are defined later in this post. ↩︎
3
Two smaller conventions, used without further comment below. An index appearing once raised and once lowered is summed over, so $A_{\mu}B^{\mu}$ means $\sum_{\mu=0}^{3}A_{\mu}B^{\mu}$. And I keep $c$ explicit throughout rather than setting it to one, so that every formula can be checked against a Newtonian limit. ↩︎
4
Einstein-Cartan theory is the standard example, where torsion is sourced by intrinsic spin and vanishes outside matter, so the vacuum predictions coincide with general relativity. ↩︎
5
The steps, for anyone checking: with $L=\tfrac12 g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$ and a dot meaning $d/d\lambda$, we have $\partial L/\partial\dot{x}^{\rho}=g_{\rho\nu}\dot{x}^{\nu}$, whose derivative along the curve is $\partial_{\mu}g_{\rho\nu}\dot{x}^{\mu}\dot{x}^{\nu}+g_{\rho\nu}\ddot{x}^{\nu}$, set against $\partial L/\partial x^{\rho}=\tfrac12\partial_{\rho}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$. The first term is contracted with something symmetric in $\mu\nu$ and may be symmetrised, and multiplying through by $g^{\sigma\rho}$ produces the Christoffel combination. ↩︎
6
Timelike geodesics extremise proper time locally and are genuinely maximal for sufficiently nearby events, but over long paths families of geodesics can refocus at conjugate points, beyond which a geodesic is no longer the longest path between its endpoints. I say extremise rather than maximise throughout for that reason. ↩︎
7
With torsion the antisymmetric part of $\Gamma^{\lambda}_{\ \mu\nu}$ survives here and contributes an extra term $-T^{\lambda}{}_{\mu\nu}\nabla_{\lambda}V^{\rho}$ to the commutator. Dropping it is a consequence of the choice made earlier, not an algebraic accident. ↩︎
8
Sign conventions for the Riemann tensor differ between textbooks, and the overall sign of this equation differs with them. If you compare against another source and the signs disagree, check its definition of $R^{\rho}_{\ \sigma\mu\nu}$ and its metric signature before concluding that one of the two is wrong. ↩︎
The Eddington luminosity is the point at which the radiation coming out of an object pushes on the surrounding gas as hard as the object’s gravity pulls it in. Above that luminosity, material near the source is pushed away instead of falling in. It sets a rough ceiling on how bright an accreting object can be, a rough ceiling on how massive a star can get, and a rough floor on how long a black hole needs to grow. The formula is short enough to memorise, which is a shame, because the derivation takes about ten lines and every one of those lines contains an assumption you can later break on purpose.
Set up the simplest possible situation. A compact object of mass $M$ sits at the origin and radiates isotropically with luminosity $L$. Around it is fully ionised hydrogen, which means free protons and free electrons and nothing else. The gas is optically thin, so a photon travelling outward has a small chance of interacting and mostly just leaves. Ask what happens to one electron-proton pair sitting at radius $r$.
Start with the radiation. Luminosity is energy per unit time, and if the source is isotropic that energy is spread evenly over a sphere of area $4\pi r^2$, so the energy flux at radius $r$ is $$F=\frac{L}{4\pi r^{2}}.$$ That’s energy per unit area per unit time. What we actually need is momentum, because forces come from momentum transfer.
A photon of energy $E$ carries momentum $p=E/c$. That follows from the relativistic energy-momentum relation $E^2=p^2c^2+m^2c^4$ with $m=0$, and it’s also what classical electromagnetism gives you for the ratio of momentum density to energy density in a plane wave. So a stream carrying energy flux $F$ also carries momentum flux $F/c$, and momentum flux has units of pressure. This quantity, $F/c$, is the radiation pressure available to push on anything that gets in the way.1It helps to see $p=E/c$ from more than one direction. In classical electromagnetism it comes from the Poynting vector and the Maxwell stress tensor: a plane wave has momentum density $S/c^2$ and energy density $S/c$, so their ratio is $1/c$. In quantum terms each photon carries $E=h\nu$ and $p=h\nu/c=h/\lambda$. Both give the same momentum flux, which is why the classical derivation of the Eddington limit needs no quantum mechanics at all.
Now we need to know how much of that momentum a single electron actually receives, which means we need a cross section. Calling a cross section an area is useful, but it shouldn’t be taken literally. The electron isn’t a small disk of area $6.65\times10^{-25}\ \mathrm{cm^2}$. A cross section packages an interaction probability into the dimensions of an area: if a beam delivers a certain number of photons per square centimetre per second, multiplying that flux by $\sigma$ gives the rate at which one electron scatters them. The area belongs to the interaction, and it tells you nothing about how big the electron is. With that understood, a particle with momentum-transfer cross section $\sigma$ sitting in a stream of momentum flux $F/c$ receives net momentum at a rate $\sigma F/c$, and a rate of momentum gain is a force.
The relevant process here is Thomson scattering, the elastic scattering of low-energy photons off free electrons, and it’s worth deriving rather than looking up because the derivation tells you exactly when it stops being valid. Take an electromagnetic wave with electric field amplitude $E_0$ passing a free electron. The field accelerates the electron with $a=eE_0\cos\omega t/m_e$. An accelerating non-relativistic charge radiates according to the Larmor formula, $P=2e^2a^2/3c^3$, and time-averaging $\cos^2$ gives a factor of $1/2$, so $$\langle P\rangle=\frac{2e^{2}}{3c^{3}}\cdot\frac{e^{2}E_0^{2}}{2m_e^{2}}=\frac{e^{4}E_0^{2}}{3m_e^{2}c^{3}}.$$ The incident flux is the time-averaged Poynting vector, $\langle S\rangle=cE_0^2/8\pi$. The cross section is the ratio of power scattered to flux incident, $$\sigma_T=\frac{\langle P\rangle}{\langle S\rangle}=\frac{8\pi}{3}\left(\frac{e^{2}}{m_ec^{2}}\right)^{2}=\frac{8\pi}{3}r_e^{2}=6.65\times10^{-25}\ \mathrm{cm^{2}},$$ where $r_e=e^2/m_ec^2=2.82\times10^{-13}$ cm is the classical electron radius.
That derivation hands us something useful for free. The cross section scales as $1/m^2$, so the same calculation applied to a proton gives a cross section smaller by $(m_e/m_p)^2=(1/1836)^2\approx3\times10^{-7}$. Protons scatter light about three million times less effectively than electrons. So when we say the radiation force acts on the electrons and ignore the protons, that’s a conclusion with a number behind it.
So the outward radiation force on one free electron is $$f_{\rm rad}=\sigma_T\,\frac{F}{c}=\frac{\sigma_T L}{4\pi r^{2}c}.$$
The geometry of the calculation. Luminosity $L$ spreads over a sphere of area $4\pi r^2$ to give an energy flux $F$, and since each photon carries momentum $E/c$ the same stream carries a momentum flux $F/c$. A free electron intercepts that stream over an effective area $\sigma_T$ and receives net momentum at a rate that constitutes a force. Nothing in this picture is quantum mechanical.
There’s something slightly suspicious about what comes next. The radiation force we just wrote contains an electron cross section. The gravitational force is about to contain a proton mass. Why are we allowed to balance two forces that act on two different particles?
The answer is electrostatics. Gravity pulls on mass, and essentially all the mass is in the protons since $m_p/m_e\approx1836$, while radiation pushes on the electrons. If the electrons started to drift outward relative to the protons, that would immediately create a charge separation, which creates an electric field, which pulls the electrons back and drags the protons out. The question is how much separation is needed, and the answer is almost none, because the Coulomb interaction is absurdly stronger than gravity. For a single proton-electron pair, $$\frac{e^{2}}{G m_p m_e}=2.3\times10^{39}.$$ A fractional charge imbalance far too small to measure supplies an electric field far larger than anything gravity or radiation is doing. The plasma therefore moves as a neutral fluid, and the correct bookkeeping is to treat a proton-electron pair as the unit: it feels the radiation force of one electron and the gravitational pull of one proton.2You can also make this concrete. To transmit a gravitational force $GMm_p/r^2$ to an electron by electrostatic means, you need a field $E$ with $eE=GMm_p/r^2$, which for a spherical geometry corresponds to a net enclosed charge of $Q=GMm_p/e$. For a solar mass that’s about $10^{21}$ elementary charges, which sounds like a lot until you notice that a single gram of ionised hydrogen contains $6\times10^{23}$ electrons. The required fractional imbalance is tiny, and the plasma establishes it essentially instantly.
Why the mass in the formula is the proton mass while the cross section is the electron’s. Radiation pushes electrons, gravity holds protons, and the Coulomb attraction between them is roughly $10^{39}$ times stronger than their mutual gravity. The pair is welded together, so the balance is written per pair rather than per particle.
The inward gravitational force per pair is then $$f_{\rm grav}=\frac{GM(m_p+m_e)}{r^{2}}\simeq\frac{GMm_p}{r^{2}}.$$ Setting the two forces equal gives $$\frac{\sigma_T L}{4\pi r^{2}c}=\frac{GMm_p}{r^{2}},$$ and the radius cancels from both sides, leaving $$L_{\rm Edd}=\frac{4\pi GMm_pc}{\sigma_T}.$$
That cancellation is easy to read past, and it’s the reason the Eddington limit comes out as a statement about luminosity instead of a statement about some particular radius. Both forces fall off as $1/r^2$: gravity because it’s an inverse square law, and radiation because the flux from a point source dilutes as $1/r^2$. Moving farther away weakens gravity and radiation pressure in exactly the same proportion. Plotted against radius on log axes, both are straight lines of slope $-2$, which means they’re parallel and never cross. The balance is either satisfied at every radius or at no radius. There’s no critical distance inside which radiation wins. There’s a critical luminosity, and above it radiation wins at every radius at once.
Why the result is a critical luminosity rather than a critical radius. Gravity and radiation pressure both dilute as $1/r^2$, so on logarithmic axes both are straight lines of slope $-2$. Changing the luminosity slides the radiation line up and down without ever tilting it, so it either sits entirely above the gravity line or entirely below it.
Putting numbers in, $$L_{\rm Edd}=1.26\times10^{38}\left(\frac{M}{M_\odot}\right)\ \mathrm{erg\,s^{-1}}=3.3\times10^{4}\left(\frac{M}{M_\odot}\right)L_\odot.$$ The Sun’s actual luminosity is about $3\times10^{-5}$ of this value, so it sits very far from the balance point and nothing in its outer layers is close to being pushed off. Massive stars are a different matter, and we’ll come back to that.
Before breaking anything, rewrite the result in a form that makes the assumptions easier to isolate. Define the opacity $\kappa$ as cross section per unit mass, so that the force per gram of material is $\kappa F/c$. For pure ionised hydrogen there’s one electron per proton, so $\kappa_{\rm es}=\sigma_T/m_p=0.40\ \mathrm{cm^2\,g^{-1}}$, and the general statement becomes $$L_{\rm Edd}=\frac{4\pi GMc}{\kappa}.$$ The composition and the physics of the interaction are now quarantined inside a single symbol. Many of the corrections below can be understood either as changing $\kappa$, or as breaking one of the assumptions that let us write this simple force balance in the first place.
Check the units of that expression once, because it shows why $\kappa$ has to be an area per unit mass and not just an area: $$[G][M][c]\,[\kappa]^{-1}=\left(\frac{\mathrm{cm^{3}}}{\mathrm{g\,s^{2}}}\right)(\mathrm{g})\left(\frac{\mathrm{cm}}{\mathrm{s}}\right)\left(\frac{\mathrm{g}}{\mathrm{cm^{2}}}\right)=\frac{\mathrm{g\,cm^{2}}}{\mathrm{s^{3}}}=\mathrm{erg\,s^{-1}}.$$ Notice that the mass in $GM$ and the mass hidden inside $\kappa$ don’t cancel. That’s gravity pulling on the gas and radiation pushing on its cross section, showing up in the units.
Composition is the easiest thing to fix. What matters for scattering is the number of free electrons per gram, and hydrogen is unusual because it supplies one electron per nucleon while everything else supplies roughly one per two nucleons. For a fully ionised mix with hydrogen mass fraction $X$ and helium mass fraction $Y=1-X$, the electron density is $n_e=(\rho/m_p)(X+Y/2)=(\rho/m_p)(1+X)/2$, so $$\kappa_{\rm es}=\frac{\sigma_T}{m_p}\frac{1+X}{2}\simeq0.20\,(1+X)\ \mathrm{cm^{2}\,g^{-1}}.$$ For solar composition, $X=0.7$, this gives $0.34\ \mathrm{cm^2\,g^{-1}}$ and $L_{\rm Edd}=1.5\times10^{38}(M/M_\odot)\ \mathrm{erg\,s^{-1}}$. For hydrogen-free material the opacity halves and the Eddington limit doubles, purely because there are half as many electrons per gram to push on.3The commonly quoted $1.26\times10^{38}(M/M_\odot)$ assumes pure hydrogen. Papers on accreting neutron stars and X-ray bursts often use the solar or the helium value instead, so a factor of two discrepancy between two sources usually means they made different composition assumptions rather than that one of them is wrong. So check which $\kappa$ a quoted Eddington luminosity is built on before comparing two of them.
There’s a distinction buried in the derivation that causes a lot of confusion later, so let’s name it now. The force balance we wrote is local. It compares two forces on one parcel of gas at one place. It only turned into a statement about the total luminosity of the source because spherical symmetry let us write the local flux as $L/4\pi r^2$ everywhere at once. So $L_{\rm Edd}$ is a global luminosity scale derived from a local condition, and the translation between the two depends entirely on the geometry. When a source is described as super-Eddington in the literature, it’s often worth asking which of the two is meant, because a source can exceed the spherical Eddington luminosity as measured by a distant observer while the radiation force on every individual parcel of gas remains sub-Eddington.
That distinction is the first structural assumption to break. The derivation put a point source at the centre of a spherical shell of gas, so radiation and gravity were both radial and could be compared directly. Real accretion usually happens through a disc, and radiation escaping from an accretion flow doesn’t come out isotropically. If the gas arrives along the equatorial plane while the radiation leaves through low-density funnels along the rotation axis, the two never have to fight over the same material, and the luminosity inferred by a distant observer can exceed $L_{\rm Edd}$ while the accretion flow itself stays locally sub-Eddington.
How anisotropy breaks the translation from local to global. The derivation compared an outward radiation flux with an inward gravitational pull along the same radial direction. In a disc geometry the inflow and the outflow occupy different solid angles, so the spherical formula no longer converts the local force balance into a single luminosity. Geometric beaming and anisotropic accretion are important ingredients in models of ultraluminous X-ray sources.
The second structural assumption is that the situation can be treated through an instantaneous force balance. The derivation compares two forces acting on the gas at a given moment, and by itself says nothing about how the gas responds with time. A transient can therefore exceed the limit for as long as it takes the surrounding material to accelerate, expand, or otherwise adjust. Classical novae and X-ray bursts on neutron stars both show episodes consistent with locally super-Eddington fluxes, and the associated observational signature is photospheric radius expansion, in which the envelope is lifted, expands, and settles back as the luminosity drops.
The third assumption, and the one that fails most often in practice, is that electron scattering is the only relevant opacity. For soft photons in a hot, fully ionised plasma, electron scattering provides a useful baseline: it’s what remains when the gas is too hot to have any bound electrons left to absorb anything. Bound-free, free-free, line and dust interactions can all add substantially to it, and since $L_{\rm Edd}\propto1/\kappa$, every addition lowers the true limit below the value the standard formula gives. But even this baseline isn’t universal, and we’ll see shortly that once photon energies approach $m_ec^2$ the Thomson approximation itself fails and the scattering opacity falls below $\kappa_{\rm es}$.
Spectral lines are the clearest example of an addition. A bound electron in an ion can absorb a photon whose energy matches a transition, and the cross section at line centre can exceed the Thomson cross section by many orders of magnitude. In the ultraviolet spectrum of a hot star there are very large numbers of such transitions from iron-group ions, and although each covers a narrow frequency range, together they intercept a substantial fraction of the flux. The effective opacity for a hot star’s wind can be tens to hundreds of times $\kappa_{\rm es}$, which is how O stars can drive strong winds while their continuum Eddington ratios remain well below unity. The star is sub-Eddington for electron scattering and super-Eddington for lines.4This is the basis of the Castor, Abbott and Klein theory of line-driven winds from 1975, which parametrises the extra push through a force multiplier applied to the electron-scattering radiation force. The multiplier depends on the ionisation state and on the velocity gradient in the wind, because a wind that’s accelerating Doppler-shifts each line into fresh continuum photons and keeps the lines from saturating. The velocity gradient is therefore part of the opacity, which is an unusual feature of the problem.
Dust is the more extreme addition. A grain absorbs across a broad continuum rather than in narrow lines, and for gas with an ordinary interstellar dust-to-gas ratio the ultraviolet opacity per gram of gas is of order a few hundred $\mathrm{cm^2\,g^{-1}}$, which is roughly a thousand times $\kappa_{\rm es}$. In dusty gas the relevant Eddington limit therefore drops by about three orders of magnitude, and objects that look comfortably sub-Eddington by the standard formula can be strongly super-Eddington for the material around them. This is why radiation pressure on dust matters for massive star formation and for feedback in galactic nuclei, and why the electron-scattering formula is a poor guide in those settings.5The estimate comes from the standard interstellar extinction per hydrogen column, $A_V/N_H\approx5\times10^{-22}$ mag cm$^2$, converted to an opacity per gram using a mean mass per hydrogen nucleus of about $1.4m_p$, which gives roughly $200\ \mathrm{cm^2\,g^{-1}}$ in the visual and two to three times that in the ultraviolet. It scales linearly with the dust-to-gas ratio, so it drops in low-metallicity environments and vanishes wherever grains have been destroyed.
Opacity for a gram of material, on a logarithmic scale. Electron scattering sits at the left because it’s what remains when the gas is hot enough to have no bound electrons. Lines and dust sit far to the right, and since the Eddington luminosity is inversely proportional to opacity, they lower the relevant limit by factors of tens to thousands. The electron-scattering value is a convenient reference scale rather than a universal ceiling.
The fourth assumption is that the photons are soft. The Thomson derivation treated the electron as a non-relativistic charge oscillating in a classical wave, which requires the photon energy to be small compared to $m_ec^2=511$ keV. Above that, the photon recoils the electron appreciably, the scattering becomes inelastic, and the Klein-Nishina formula replaces the Thomson result. The cross section falls, roughly as $\ln(2x)/x$ for $x=h\nu/m_ec^2\gg1$. Hard photons are therefore less effective at pushing gas, and this is the one correction that pushes the limit up rather than down: a source radiating mainly above a few hundred keV can exceed the naive Eddington luminosity simply because its photons couple weakly to electrons.
The scattering cross section as a function of photon energy, in units of the Thomson value. The Thomson result is the flat portion on the left, valid while $h\nu\ll m_ec^2$. By $100$ keV the cross section is down by about a quarter, at $511$ keV it’s below half, and at several MeV it’s an order of magnitude smaller. This is the reason electron scattering cannot be treated as a hard floor on the opacity.
The fifth assumption is Newtonian gravity, which is uncomfortable given that the objects people apply this to are usually black holes and neutron stars. In the Schwarzschild metric two things change for a static observer at radius $r$. The local gravitational acceleration picks up a factor $(1-r_s/r)^{-1/2}$, and the flux from a source of luminosity $L_\infty$ as measured at infinity picks up a factor $(1-r_s/r)^{-1}$ from gravitational redshift and time dilation together. Balancing them gives $$L_{\rm Edd}^{\infty}=\frac{4\pi GMc}{\kappa}\sqrt{1-\frac{r_s}{r}}.$$ The two corrections partly cancel, which is why the Newtonian answer survives as well as it does, but not completely. At $3r_s$ the limit is $82$ percent of the flat-space value and at $1.5r_s$ it’s $58$ percent. Note also that the answer now depends on $r$, so the radius no longer cancels and the earlier statement about parallel lines fails.
The sixth assumption is that the gas is ordinary matter with protons in it. Within the same simple force-balance picture, an electron-positron pair plasma changes the bookkeeping completely: there are no protons carrying most of the mass, so gravity acts on a mass of order $m_e$ per scattering particle rather than $m_p$, while the scattering cross section remains of order $\sigma_T$. The corresponding Eddington scale therefore drops by roughly $m_e/m_p\approx1/1836$. Real pair plasmas involve additional physics such as pair creation, annihilation and more complicated radiative transfer, so this factor is what you get by pushing the original argument into a new regime, and the real problem carries more than that. Pair plasmas are thought to be relevant in the inner regions of some accretion flows and in gamma-ray burst fireballs, where corrections of this size are not small.
The seventh assumption is that the gas is smooth. The derivation implicitly treats the medium as uniform, so every gram sees the same flux. If the gas is clumpy, radiation escapes preferentially through the low-density channels between the clumps and deposits less momentum than a uniform calculation predicts. The effective opacity of a porous medium is lower than that of a smooth one with the same mean density, which raises the effective limit. This idea, usually called porosity, is one of the mechanisms proposed for the eruptions of luminous blue variables, which appear to radiate above their nominal limits for extended periods.
With the limit itself understood, two more quantities follow, and these are what the Eddington luminosity actually gets used for in practice. If an accreting object converts a fraction $\eta$ of the rest mass energy of the infalling material into radiation, so that $L=\eta\dot Mc^2$, then the accretion rate corresponding to the Eddington luminosity is $$\dot M_{\rm Edd}=\frac{L_{\rm Edd}}{\eta c^{2}}=\frac{4\pi GM}{\eta\kappa c}\approx2.2\times10^{-8}\left(\frac{0.1}{\eta}\right)\left(\frac{M}{M_\odot}\right)M_\odot\,\mathrm{yr^{-1}}.$$ The efficiency $\eta$ is about $0.057$ for a non-rotating black hole and up to about $0.42$ for a maximally rotating one, so this number carries a factor of several of uncertainty depending on spin.
The second derived quantity needs one piece of care that’s easy to miss. If $\dot M$ is the rest mass supplied to the hole and a fraction $\eta$ of it is radiated away, then the mass actually retained is $\dot M_{\rm BH}=(1-\eta)\dot M$. Carrying that through, the e-folding time for the black hole mass under Eddington-limited accretion is $$t_S=\frac{\eta}{1-\eta}\,\frac{\kappa c}{4\pi G}\approx50\ \mathrm{Myr}\quad(\eta=0.1),$$ rather than the $45$ Myr you get by dropping the $(1-\eta)$. The distinction is small at low efficiency and large at high efficiency: for a maximally spinning hole with $\eta=0.42$ it changes the answer from $189$ to $326$ Myr.6The $45$ Myr form is the one usually quoted, and for order-of-magnitude arguments the difference rarely matters. Check which version a paper is using when the argument turns on a factor of two, which it sometimes does in discussions of early black hole growth.
That number is behind a well-known problem. Growing a $10\ M_\odot$ seed into a $10^9\ M_\odot$ quasar takes about $18$ e-foldings, which at $50$ Myr each is roughly $900$ Myr of continuous Eddington-limited accretion with no idle periods at all. Quasars of that mass are observed less than a billion years after the Big Bang, so the budget is uncomfortably tight, which is why heavier seeds and episodes of super-Eddington accretion are both actively discussed.
The limit bites hardest in massive stars. Main sequence luminosity rises steeply with mass, faster than linearly, while the Eddington luminosity rises exactly linearly, so the ratio $\Gamma=L/L_{\rm Edd}$ climbs with mass. Using standard zero-age main sequence luminosities, a $10\ M_\odot$ star has $\Gamma\approx0.03$ and a $100\ M_\odot$ star roughly $0.3$ on the continuum opacity alone. Add line opacity and the effective ratio approaches unity, which is connected both to the large mass loss rates of the most massive stars and to the existence of an upper limit on stellar masses. For these objects the Eddington limit is doing real work on the structure of the star.
References and Footnotes
1
It helps to see $p=E/c$ from more than one direction. In classical electromagnetism it comes from the Poynting vector and the Maxwell stress tensor: a plane wave has momentum density $S/c^2$ and energy density $S/c$, so their ratio is $1/c$. In quantum terms each photon carries $E=h\nu$ and $p=h\nu/c=h/\lambda$. Both give the same momentum flux, which is why the classical derivation of the Eddington limit needs no quantum mechanics at all. ↩︎
2
You can also make this concrete. To transmit a gravitational force $GMm_p/r^2$ to an electron by electrostatic means, you need a field $E$ with $eE=GMm_p/r^2$, which for a spherical geometry corresponds to a net enclosed charge of $Q=GMm_p/e$. For a solar mass that’s about $10^{21}$ elementary charges, which sounds like a lot until you notice that a single gram of ionised hydrogen contains $6\times10^{23}$ electrons. The required fractional imbalance is tiny, and the plasma establishes it essentially instantly. ↩︎
3
The commonly quoted $1.26\times10^{38}(M/M_\odot)$ assumes pure hydrogen. Papers on accreting neutron stars and X-ray bursts often use the solar or the helium value instead, so a factor of two discrepancy between two sources usually means they made different composition assumptions rather than that one of them is wrong. So check which $\kappa$ a quoted Eddington luminosity is built on before comparing two of them. ↩︎
4
This is the basis of the Castor, Abbott and Klein theory of line-driven winds from 1975, which parametrises the extra push through a force multiplier applied to the electron-scattering radiation force. The multiplier depends on the ionisation state and on the velocity gradient in the wind, because a wind that’s accelerating Doppler-shifts each line into fresh continuum photons and keeps the lines from saturating. The velocity gradient is therefore part of the opacity, which is an unusual feature of the problem. ↩︎
5
The estimate comes from the standard interstellar extinction per hydrogen column, $A_V/N_H\approx5\times10^{-22}$ mag cm$^2$, converted to an opacity per gram using a mean mass per hydrogen nucleus of about $1.4m_p$, which gives roughly $200\ \mathrm{cm^2\,g^{-1}}$ in the visual and two to three times that in the ultraviolet. It scales linearly with the dust-to-gas ratio, so it drops in low-metallicity environments and vanishes wherever grains have been destroyed. ↩︎
6
The $45$ Myr form is the one usually quoted, and for order-of-magnitude arguments the difference rarely matters. Check which version a paper is using when the argument turns on a factor of two, which it sometimes does in discussions of early black hole growth. ↩︎
Making a star means compressing gas by something like twenty-four orders of magnitude. The warm neutral interstellar medium, the ordinary atomic gas between the stars, holds roughly a third of a hydrogen atom per cubic centimetre at a temperature of several thousand kelvin, which puts its density near $10^{-24}$ grams per cubic centimetre. The mean density of the Sun is $1.41$ grams per cubic centimetre. The Galaxy takes material from one end of that range and pushes it to the other, currently at a rate of one or two solar masses per year, and it has been doing so for most of its life. The question is what makes a patch of thin interstellar gas stop being gas and become a star.
The problem as a density ladder. Each rung is a physically different state of the gas. Gravity matters throughout, but what resists it changes as the collapse proceeds. Thermal pressure, turbulence and magnetic fields dominate at low density. Cooling and radiative transfer become increasingly important as the gas gets denser. Eventually the gas becomes opaque, heats up, dissociates molecular hydrogen and forms a protostar.
Gravity is the only force available that can gather matter over large distances, and it has to compete with pressure. Pressure doesn’t know in advance that a cloud is trying to collapse. It reacts to compression by launching sound waves. So the first question is a race: can gravity make a density enhancement grow faster than pressure can communicate across it and smooth it away? Jeans turned that into an equation in 1902.
Take an isothermal self-gravitating fluid, $$\frac{\partial\rho}{\partial t}+\nabla\!\cdot(\rho\mathbf{v})=0, \qquad \frac{\partial\mathbf{v}}{\partial t}+(\mathbf{v}\cdot\nabla)\mathbf{v}=-\frac{\nabla P}{\rho}-\nabla\Phi, \qquad \nabla^{2}\Phi=4\pi G\rho,$$ with $P=c_s^{2}\rho$. Take a uniform static background of density $\rho_0$, perturb it as $\rho=\rho_0+\rho_1$, $\mathbf v=\mathbf v_1$, $\Phi=\Phi_0+\Phi_1$, and keep only terms first order in the perturbations.1There is a genuine inconsistency hidden here, traditionally called the Jeans swindle. An infinite medium with constant density cannot simultaneously satisfy $\nabla^2\Phi_0=4\pi G\rho_0$ and have no background gravitational acceleration. The usual calculation ignores the gravity of the uniform background and keeps the gravity of the perturbation. More careful treatments using finite systems or expanding backgrounds recover essentially the same instability criterion, so the answer survives even though the original setup is not mathematically self-consistent.
Differentiate the continuity equation with respect to time, take the divergence of the momentum equation, and use Poisson’s equation to eliminate $\Phi_1$. The three equations become one, $$\frac{\partial^{2}\rho_1}{\partial t^{2}}-c_s^{2}\nabla^{2}\rho_1-4\pi G\rho_0\rho_1=0,$$ and looking for a plane wave $\rho_1\propto\exp[i(\mathbf k\cdot\mathbf x-\omega t)]$ collapses the whole problem to $$\omega^{2}=c_s^{2}k^{2}-4\pi G\rho_0.$$
Take the terms one at a time. Without gravity this is $\omega=c_sk$, an ordinary sound wave. Gravity subtracts the same amount from $\omega^2$ at every wavelength, while pressure’s contribution grows as $k^2$ and so becomes increasingly effective on small scales. At large $k$, meaning short wavelengths, pressure wins, $\omega^2>0$, and the perturbation oscillates. At small $k$ gravity wins, $\omega^2<0$, $\omega$ becomes imaginary, and one of the two solutions grows exponentially.
The Jeans dispersion relation plotted against $k^2$. Gravity contributes the constant term $-4\pi G\rho_0$, while pressure contributes $c_s^2k^2$. At $k=k_J$ the two exactly cancel. Longer wavelength perturbations have $\omega^2<0$ and grow exponentially, while shorter wavelength perturbations remain ordinary sound waves.
The boundary sits at $k_J=\sqrt{4\pi G\rho_0}/c_s$, or in wavelength $$\lambda_J=\frac{2\pi}{k_J}=c_s\sqrt{\frac{\pi}{G\rho_0}}.$$ In this idealised medium, perturbations larger than $\lambda_J$ are gravitationally unstable and smaller ones behave like sound waves. Real molecular clouds are neither uniform nor static nor unmagnetised, so $\lambda_J$ shouldn’t be treated as a sharp boundary in nature. What survives is the competition it exposes: pressure grows stronger on small scales, while self-gravity doesn’t weaken when you look at a larger piece of gas.
Turning the Jeans length into a mass gives the more useful quantity. Taking a sphere whose diameter is $\lambda_J$, $$M_J=\frac{4\pi}{3}\rho_0\left(\frac{\lambda_J}{2}\right)^{3}=\frac{\pi^{5/2}}{6}\frac{c_s^{3}}{G^{3/2}\rho_0^{1/2}}\approx2.9\,\frac{c_s^{3}}{G^{3/2}\rho_0^{1/2}}.$$ The prefactor matters less than the scaling, $M_J\propto c_s^3\rho^{-1/2}\propto T^{3/2}\rho^{-1/2}$, and you can recover that scaling without any perturbation theory at all. The only combination of $G$, $c_s$ and $\rho$ with units of mass is $c_s^3G^{-3/2}\rho^{-1/2}$. The full calculation tells you why that combination matters and supplies the numerical factor.
Put in numbers for a cold dense core. Molecular gas with helium has a mean molecular weight $\mu\simeq2.33$, so at $T=10$ K the sound speed is $c_s=\sqrt{k_BT/\mu m_H}\simeq0.19\ \mathrm{km\,s^{-1}}$. At $n=10^4\ \mathrm{cm^{-3}}$ this gives $\lambda_J\simeq0.21$ pc and $M_J\simeq2.9\,M_\odot$. A few solar masses spread over a few tenths of a parsec isn’t an absurd theoretical scale. It’s close to the scale of the dense structures actually found inside nearby molecular clouds.
The temperature dependence matters even more. Since $M_J\propto T^{3/2}$, going from $8000$ K to $10$ K lowers the characteristic unstable mass by $(800)^{3/2}\approx2\times10^4$. Warm diffuse gas is very hard to make gravitationally unstable on stellar mass scales, and cold gas is much easier, which ties star formation directly to cooling.
The details are a little more interesting than saying gas must become molecular in order to cool. Atomic species such as ionised carbon cool neutral gas efficiently, and molecular hydrogen is not always the dominant coolant in ordinary Galactic molecular clouds. What really matters is shielding from the interstellar ultraviolet field. Once enough column density builds up, photoelectric heating weakens, molecules survive, dust becomes cold, and temperatures near $10$ K become possible. Becoming cold and becoming molecular happen together, though one is not simply the cause of the other.
The Taurus molecular cloud in far infrared dust emission. Real molecular clouds look nothing like the uniform medium used in the Jeans calculation. Much of their dense gas is arranged into filaments, and dense cores are often found along those filaments. Herschel observations led to the suggestion that many nearby filaments have an inner width around $0.1$ pc. Whether that width is truly universal is still debated, because resolution, distance, tracer choice and the way a filament is defined can all change the measured width.2Image credit: ESA/Herschel/PACS, SPIRE/Gould Belt survey Key Programme. The characteristic $0.1$ pc width was emphasised in the Herschel Gould Belt work, but later studies have questioned how universal it really is. The safe statement is that filaments are common and closely connected to dense core formation. The exact distribution of their widths remains an active observational problem.
There’s another way to phrase the same stability problem that sits closer to what an observed core actually looks like. A real core is finite and can be confined by the pressure of the gas around it, and an isothermal pressure-confined sphere has a maximum stable mass, the Bonnor-Ebert mass, $$M_{\rm BE}\simeq1.18\,\frac{c_s^{4}}{G^{3/2}P_{\rm ext}^{1/2}}.$$ Above that, no hydrostatic configuration exists for the assumed temperature and external pressure. The details differ from the infinite-medium Jeans calculation, but the message is nearly the same. A cold enough and massive enough concentration of gas cannot find a pressure-supported equilibrium.
Suppose it crosses that line. Strip away pressure completely, start with a uniform sphere, and a shell at initial radius $r_0$ feels only the mass inside it, so $\ddot r=-GM/r^2$. Integrating once gives $\tfrac12\dot r^2=GM(1/r-1/r_0)$, and a further integration using $r=r_0\cos^2\xi$ gives the free-fall time $$t_{\rm ff}=\sqrt{\frac{3\pi}{32G\rho_0}}.$$ The radius has disappeared. For a uniform sphere the collapse time depends only on density, so every shell arrives at the centre together. At $n=10^4\ \mathrm{cm^{-3}}$ this is $0.34$ Myr, and at the more typical mean density of a molecular cloud, $n\simeq100\ \mathrm{cm^{-3}}$, about $3.4$ Myr.
Those numbers fail badly against observation. The Milky Way contains of order $10^9\,M_\odot$ of molecular gas, and if all of it turned into stars on a cloud free-fall time the star formation rate would be roughly $$\dot M_\star\sim\frac{M_{\rm mol}}{t_{\rm ff}}\sim\frac{10^9\,M_\odot}{3.4\ \mathrm{Myr}}\sim300\ M_\odot\,\mathrm{yr^{-1}}.$$ The real Galactic rate is around $1$ to $2\,M_\odot\,\mathrm{yr^{-1}}$, so the naive estimate misses by roughly two orders of magnitude.
A convenient way to write the discrepancy is the efficiency per free-fall time, $\epsilon_{\rm ff}\equiv\dot M_\star t_{\rm ff}/M_{\rm gas}$. On molecular-cloud scales the inferred values are of order one percent, often somewhere between $0.3$ and $3$ percent, with substantial scatter from cloud to cloud and through a cloud’s lifetime. The point isn’t that nature picked exactly $0.0100$. It’s that the answer is nowhere near unity, and something is preventing most molecular gas from falling straight in.3The low value of $\epsilon_{\rm ff}$ is one of the central empirical constraints on star formation theory. Modern cloud-lifecycle work also makes clear that individual clouds need not form stars at a perfectly constant rate. A cloud can evolve through quiet and active phases, so an average value near one percent does not mean every cloud converts exactly one percent of its mass every free-fall time.
The usual bookkeeping device is the virial theorem. The full scalar version contains kinetic energy, gravitational energy, magnetic terms, surface pressure and time-dependent terms, but ignore most of those and approximate the cloud as a uniform sphere. Then $W=-\tfrac35 GM^2/R$ and $T=\tfrac32 M\sigma^2$, with $\sigma$ the one-dimensional velocity dispersion, and their ratio gives the virial parameter $$\alpha_{\rm vir}\equiv\frac{5\sigma^{2}R}{GM}.$$ For this simple model $\alpha_{\rm vir}=1$ corresponds to $2T=|W|$.
Values of order unity tell us kinetic and gravitational energies are comparable. They don’t prove a cloud is sitting in static equilibrium. A cloud can be accreting, contracting, losing material, confined by external pressure or threaded by a dynamically important magnetic field and still have an order-unity virial parameter. That matters because observed molecular clouds commonly have $\alpha_{\rm vir}$ of order one or a few. They aren’t behaving like cold pressureless spheres in free fall, but the virial parameter alone doesn’t tell us what they are doing.
The most obvious extra ingredient is turbulence. At $10$ K the sound speed is only $0.19\ \mathrm{km\,s^{-1}}$, yet molecular lines are routinely much broader than that, and the non-thermal dispersion increases with the size of the region measured. Larson identified this in 1981, and later CO surveys found relations roughly of the form $$\sigma_v\sim0.7\left(\frac{R}{1\ \mathrm{pc}}\right)^{1/2}\ \mathrm{km\,s^{-1}},$$ although neither the slope nor the normalisation is universal and surface density and environment both matter. Taken literally at $R=10$ pc it gives a few kilometres per second, more than ten times the sound speed of cold molecular gas.
Supersonic turbulence then does two things at once. It resists global collapse by supplying kinetic energy, and it also creates the structures that collapse first, because supersonic flows form shocks and shocks compress gas. In idealised isothermal turbulence without self-gravity the density distribution is approximately lognormal, and once self-gravity becomes important the densest collapsing regions develop a power-law tail. Those are the regions where the local Jeans mass becomes small and collapse runs ahead. Turbulence moves mass around, builds filaments and sheets, produces rare high-density regions and changes where gravity first gets to win, which is why many theories of the star formation rate begin by asking what fraction of a turbulent density field becomes gravitationally unstable.
Magnetic fields complicate the picture again, because a field doesn’t resist all compression equally. Gas moves much more easily along a field than across it, and under ideal flux freezing matter and flux move together, so the useful quantity is the mass-to-flux ratio $M/\Phi$. For a simple sheet-like geometry the critical value is of order $$\left(\frac{M}{\Phi}\right)_{\rm crit}\simeq\frac{1}{2\pi\sqrt G}, \qquad M_\Phi\simeq\frac{\Phi}{2\pi\sqrt G}\sim\frac{R^{2}B}{2\sqrt G}.$$ At $R=1$ pc and $B=10\,\mu\mathrm{G}$ that is about $93\,M_\odot$.
If a cloud is magnetically subcritical, the field can prevent collapse across the field lines in the idealised flux-frozen limit. If it’s supercritical, gravity is strong enough that the field can’t hold the cloud up forever on its own. That doesn’t make the field irrelevant, since a supercritical field can still change the geometry, the collapse rate, the fragmentation and the angular momentum of the gas.
Why the mass-to-flux ratio matters. Gas can move relatively easily along magnetic field lines, so even a subcritical cloud can become flattened. What the field strongly resists is collapse across the field. In a supercritical region gravity is strong enough to pull material across the field lines, dragging them inward as the collapse proceeds.
Observationally the situation is messy. Zeeman splitting gives the line-of-sight field strength, dust polarisation gives the field orientation projected on the sky, and neither gives the full three-dimensional field without assumptions. The broad picture is that magnetic fields are dynamically important and that dense star-forming gas is generally not enormously subcritical. Exactly how close clouds and cores sit to the critical value is harder to establish than a single number makes it sound.
Flux freezing isn’t exact either. Molecular gas is only weakly ionised, so charged particles couple directly to the field while neutrals drift relative to them. This is ambipolar diffusion, and it lets matter move inward without dragging the full flux along, so the mass-to-flux ratio of a dense region can increase. For a long time slow ambipolar diffusion through a magnetically supported cloud was one of the standard explanations for inefficient star formation. It still matters, but that clean picture has been replaced by a mixed one in which turbulence, gravity, ambipolar diffusion, Ohmic dissipation, the Hall effect, accretion from the surrounding cloud and stellar feedback all operate somewhere in the problem, with the dominant one depending on density and environment.
Once a core does collapse, another idealised calculation becomes useful. A uniform sphere collapses everywhere at once because every shell shares a free-fall time, but real dense cores are centrally concentrated. The classic analytic model is Shu’s inside-out collapse of a singular isothermal sphere, which starts from $\rho(r)=c_s^2/2\pi Gr^2$. Once collapse begins at the centre an expansion wave moves outward at the sound speed, gas inside the wave falls inward, gas outside hasn’t yet responded, and the accretion rate is constant: $$\dot M=0.975\,\frac{c_s^{3}}{G}\simeq1.5\times10^{-6}\,M_\odot\,\mathrm{yr^{-1}}\quad(T=10\ \mathrm{K}).$$ A solar mass at that rate takes several hundred thousand years to assemble.
The scaling follows from dimensions again, since the only variables in the isothermal problem are $c_s$ and $G$ and the combination with units of mass per time is $c_s^3/G$. Warmer gas therefore has a much larger characteristic accretion rate. The number $0.975$ belongs to one very specific initial condition, though. Real cores don’t begin as exact singular isothermal spheres and real accretion isn’t constant. Simulations and observations both find time dependence, bursts and strong effects from rotation, magnetic fields and the surrounding environment. Shu’s solution is useful because it sets a clean scale, not because protostars follow it like a timetable.
The infalling material also doesn’t jump straight from molecular-cloud density to stellar density. At first the gas radiates away most of the energy released by compression, so the collapse stays close to isothermal. Eventually the central region becomes optically thick to its own infrared radiation, cooling can no longer keep up, the temperature rises, pressure stiffens, and collapse temporarily stops. A small hydrostatic object a few astronomical units across appears: the first hydrostatic core, predicted by Larson in 1969. Finding one has been much harder, since the phase is short-lived and deeply buried, and although several objects have been proposed, an unambiguous observational identification remains difficult.4Modern radiation-hydrodynamic calculations give a detailed picture of the first and second core stages, but the observational signature of a genuine first hydrostatic core is not unique. Several good candidates exist. It is safer to call them candidates than to say the phase has been securely observed as a class of objects.
The first core doesn’t survive. As its centre approaches roughly $2000$ K, molecular hydrogen begins to dissociate, and breaking an $\mathrm{H_2}$ molecule costs $4.48$ eV. Energy from compression now goes into dissociation instead of efficiently raising the temperature, the effective adiabatic index drops, pressure support weakens, and the gas enters a second collapse. That continues until most of the hydrogen in the centre has become atomic and the equation of state stiffens again. What’s left behind is the second hydrostatic core, which is what we normally call the protostar.
This change in thermodynamics also matters for fragmentation. While gas cools efficiently a collapsing region can keep breaking into smaller unstable pieces, and once it becomes opaque and starts heating rapidly, further fragmentation becomes much harder. That produces an opacity-limited minimum fragment scale. The exact mass depends on the thermodynamics, the rotation and the later accretion, so it shouldn’t be confused with a hard lower limit on the final mass of a brown dwarf. It is the first hint that microscopic physics can leave a preferred mass scale inside a problem that looked scale-free when we started with Jeans.
The Pillars of Creation in M16 seen by JWST’s NIRCam. The pillars are dense parts of a molecular cloud being eroded by ultraviolet radiation from nearby massive stars. Young stars are still forming inside them, and some of the small knots and wave-like structures are associated with outflows from embedded young objects. Star formation and destruction are happening in the same piece of cloud at the same time.5Image credit: NASA, ESA, CSA, STScI; image processing by Joseph DePasquale, Anton M. Koekemoer and Alyssa Pagan. JWST NIRCam image released in October 2022.
There’s another problem waiting for the collapsing gas, and it has nothing to do with whether gravity is strong enough. The gas rotates. A core with radius $0.1$ pc and a characteristic rotational speed of $0.1\ \mathrm{km\,s^{-1}}$ has a specific angular momentum of order $j_{\rm core}\sim Rv\sim3\times10^{21}\ \mathrm{cm^2\,s^{-1}}$, while the specific angular momentum actually stored in the Sun’s spin is only of order $10^{15}\ \mathrm{cm^2\,s^{-1}}$. The precise comparison depends on the internal rotation profile and on how you define the core’s rotation, but the mismatch is around six orders of magnitude, and angular momentum can’t simply be destroyed.
The angular momentum problem on a logarithmic scale. A dense molecular core contains millions of times more specific angular momentum than the spin of the star it can eventually produce. That angular momentum cannot disappear, so it has to be transferred elsewhere. The Solar System gives a familiar example: the Sun contains almost all the mass, while the planetary orbits contain most of the angular momentum. During star formation, discs and outflows provide the main routes for that redistribution.
The first part of the solution is a disc. Gas with too much angular momentum to fall directly onto the centre reaches a centrifugal barrier and settles into rotation, with $\Omega(R)=\sqrt{GM_\star/R^3}$ for a Keplerian disc. Making a disc doesn’t solve anything by itself, because if every ring keeps its angular momentum forever the disc just sits there. For matter to move inward, angular momentum has to move outward.
Shakura and Sunyaev gave a way to describe that transport without knowing its cause, writing the effective viscosity as $\nu=\alpha c_sH$ with $H$ the disc scale height and $\alpha$ packaging the unknown physics into one dimensionless number. For years the most famous candidate for supplying it was the magnetorotational instability. Take two pieces of a differentially rotating disc and connect them weakly with a field line. The inner piece orbits faster, magnetic tension transfers angular momentum from the inner piece to the outer one, the inner piece falls further inward while the outer one moves further out, and the increased differential motion lets the instability grow.
The difficulty is that a protoplanetary disc isn’t an ideal conducting fluid. Much of it is cold and only weakly ionised, and Ohmic resistivity, ambipolar diffusion and the Hall effect all modify the magnetic coupling, suppressing vigorous MRI turbulence over large parts of the disc. The modern picture puts much more weight on magnetised disc winds, where large-scale fields threading the disc carry angular momentum vertically away from it and let gas accrete without strong turbulence throughout. MRI probably still matters in some regions and gravitational torques matter in young massive discs, but there’s no single universal $\alpha$ mechanism.6The history is worth noting. Velikhov and Chandrasekhar studied the instability decades before Balbus and Hawley recognised its importance for accretion discs in 1991, after which MRI became almost synonymous with angular momentum transport in disc theory. Non-ideal MHD changed that. Current work increasingly treats magnetic winds as a major and possibly dominant route for angular momentum loss in protoplanetary discs, with the balance between winds, MRI, hydrodynamic instabilities and gravitational torques depending on radius and evolutionary stage.
The disc around HL Tauri in millimetre continuum emission from ALMA. The bright rings and dark gaps showed that strong disc substructure can appear astonishingly early. Planets are one possible explanation for some of these gaps, though not the only one, since dust growth, pressure structures and other disc processes can also produce rings. The important point here is that the smooth-disc approximation breaks down very quickly in real systems.7Image credit: ALMA (ESO/NAOJ/NRAO), released in 2014 from ALMA’s long-baseline campaign.
A young disc can also become gravitationally unstable on its own, and the disc version of the Jeans argument is Toomre’s criterion. For a thin gas disc with surface density $\Sigma$, sound speed $c_s$ and epicyclic frequency $\kappa$, define $$Q=\frac{c_s\kappa}{\pi G\Sigma},$$ with $\kappa=\Omega$ for a Keplerian disc. Axisymmetric disturbances become gravitationally unstable when $Q$ drops below order unity.
Instability doesn’t automatically mean fragmentation. A disc with $Q\sim1$ can develop spiral structure and use gravitational torques to transport angular momentum without breaking into separate objects. To fragment it also has to shed the heat generated by compression and shocks quickly enough, which is why disc thermodynamics matters as much as the value of $Q$. Young massive discs are therefore the most vulnerable, and if fragmentation does occur it can make stellar or brown-dwarf companions. Forming giant planets this way is possible in some parts of parameter space, especially at large radii, but it isn’t the inevitable outcome of $Q<1$.[mfn]A common way to express the cooling requirement is through $t_{\rm cool}\Omega$. Fragmentation is favoured when the cooling time is only a few orbital times or less, although the exact critical value is not universal and depends on the equation of state, irradiation and numerical setup.[/mfn]
The other half of the angular momentum story sits above and below the disc, where young stars launch jets and winds. The fastest collimated jets reach hundreds of kilometres per second and extend for parsecs, and their connection to accretion isn't an accident, because a magnetic field can extract angular momentum from rotating disc material and put it into an outflow. The cleanest toy model is Blandford and Payne's bead on a wire: take a field line anchored to the disc at radius $r_0$, pretend it's rigid, and force a parcel of gas to slide along it while rotating with the footpoint. In the rotating frame the effective potential is $$\Phi_{\rm eff}(r,z)=-\frac{GM}{\sqrt{r^{2}+z^{2}}}-\frac12\Omega_0^{2}r^{2},$$ and expanding it along the field line near the disc shows that if the line leans more than $30^\circ$ from the rotation axis the potential falls as you move outward along it.
Gas can then accelerate away from the surface without thermal pressure pushing it over a barrier. The wind doesn't get its energy for free: the energy and angular momentum come from the rotating accretion flow, and the magnetic field is the lever that transfers them. Further out the gas passes the Alfvén point and the field can no longer hold it in exact corotation, but by then the outflow has taken angular momentum with it, which is what the accreting material needed.
The magnetocentrifugal launching criterion in the idealised bead-on-a-wire picture. The angle is measured from the direction parallel to the rotation axis. A field line inclined by more than $30^\circ$ provides a downhill path in the effective potential, allowing gas to accelerate away from the disc. The outflow draws its energy and angular momentum from the rotating accretion flow.
For most of the history of astronomy almost all of this happened where we couldn’t see it, because a young protostar sits inside an optically thick dusty envelope and visible light tells you very little. Infrared and submillimetre radiation do much better, since dust absorbs short-wavelength radiation and reradiates the energy at longer wavelengths. That led to classifying young stellar objects by the shape of their spectral energy distributions, using the infrared spectral index $$\alpha=\frac{d\log(\lambda F_\lambda)}{d\log\lambda}$$ measured across the near and mid infrared.
A Class I source has a rising infrared SED, roughly $\alpha>0.3$. Flat-spectrum sources sit in the range $-0.3\lesssim\alpha\lesssim0.3$, Class II sources have $-1.6\lesssim\alpha\lesssim-0.3$, and Class III sources have still more negative slopes. Class 0 is different. It was added later and isn’t properly defined by the same infrared slope, because a genuine Class 0 source can be almost invisible at the wavelengths used to measure $\alpha$. It’s identified instead from its very cold SED and from the fact that a large fraction of its luminosity emerges in the submillimetre, with most of the system’s mass still in the envelope.
The sequence roughly tracks evolution. Class 0 and I objects are embedded protostars, by Class II the envelope is mostly gone and a circumstellar disc produces the infrared excess, and by Class III the excess is weak and most of the primordial disc has disappeared. Class and age aren’t the same thing, though. Turn a disc edge-on and it can look far more embedded than it is, and extra foreground extinction changes the slope again. The classes are useful observational labels, and geometry can make one evolutionary stage masquerade as another.
The observational classes of young stellar objects. The sequence broadly follows the removal of the envelope and then the disc, but it isn’t a perfect clock. Class 0 is defined using the cold, envelope-dominated SED rather than the same near-to-mid-infrared slope used for Classes I, Flat, II and III. Inclination and extinction can also move a source between apparent classes.
Once the main accretion phase ends, a low-mass object still isn’t on the main sequence. It’s larger than an ordinary star of the same mass and is still contracting, with gravitational energy of order $E_{\rm grav}\sim GM^2/R$ available, so dividing by the luminosity gives the Kelvin-Helmholtz time $$t_{\rm KH}\sim\frac{GM^{2}}{RL}\sim3\times10^{7}\ \mathrm{yr}$$ for present-day solar values.
A young solar-mass star first moves down a Hayashi track in the Hertzsprung-Russell diagram, staying at roughly constant effective temperature while its luminosity falls, and it’s mostly convective. As a radiative core develops it moves onto a Henyey track and heads toward the main sequence. Deuterium burning starts earlier, at central temperatures of order $10^6$ K, and acts as a temporary thermostat, while sustained hydrogen burning through the proton-proton chain needs of order $10^7$ K. Once nuclear burning replaces the energy radiated from the surface, contraction largely stops. Below roughly $0.075$ to $0.08\,M_\odot$, depending slightly on composition, the centre never gets hot enough, electron degeneracy pressure becomes important first, and the result is a brown dwarf that spends the rest of its life cooling.
So far this has been about how one collapsing region makes one central object. The harder question is why it makes the masses it does. Count newly formed stars by mass and you get the initial mass function, which above roughly a solar mass approaches the Salpeter form $dN/dM\propto M^{-2.35}$, or equivalently $dN/d\log M\propto M^{-1.35}$. Below a solar mass the distribution flattens and turns over, with a characteristic mass of a few tenths of a solar mass.
The broad shape is remarkably similar across many nearby star-forming regions, which is one of the strangest facts in the subject, since clouds vary in density, temperature, turbulence, radiation field and chemical composition and the stellar mass distribution doesn’t respond as violently as you’d expect. That doesn’t mean universality has been proved everywhere. There are continuing claims of variation, especially in extreme environments and in unresolved old stellar populations. The safer statement is that large variations are surprisingly hard to produce in the environments where we can measure individual young stars well.
A schematic initial mass function plotted as number per logarithmic mass interval. The distribution turns over at a characteristic mass of a few tenths of a solar mass. Below about $0.075\,M_\odot$ lie the brown dwarfs, while at high masses the IMF approaches the Salpeter form, $dN/d\log M\propto M^{-1.35}$. Low-mass stars dominate the number count, while the rare massive stars dominate the luminosity, ionising radiation and eventual supernova feedback.
Why the IMF has this shape is still not settled. For the turnover, one line of argument has become more specific in recent years. The Jeans mass by itself isn’t enough, because it changes continuously as density and temperature change. But the thermodynamics of collapse contains a point where the gas stops behaving nearly isothermally and begins to heat strongly as it becomes opaque, and that change produces the first hydrostatic core and introduces a characteristic mass tied to dust opacity and molecular hydrogen physics. Since those microphysical scales don’t change wildly across ordinary Galactic star-forming environments, they may help explain why the peak is fairly stable.8This is not the only proposed explanation for the IMF peak. Earlier arguments often tied the characteristic mass to the density at which gas and dust thermally couple, while modern calculations increasingly emphasise the thermodynamic transition associated with opacity and the first hydrostatic core. The broader point is the same: some piece of non-scale-free microphysics has to enter if the IMF is going to acquire a preferred mass.
The high-mass power law probably has a different origin. Gravity and turbulence are nearly scale-free over a large range and both naturally generate power-law structures, while turbulent fragmentation, continuing accretion, interactions between neighbouring collapsing objects and stellar feedback can all change the final mass spectrum. Current simulations produce IMFs that look impressively realistic, but there’s still no short derivation that starts from cloud properties and predicts the whole function.
That becomes clearer once massive stars enter. A massive protostar doesn’t finish accreting and then switch on. Its Kelvin-Helmholtz time is so short that it can begin sustained hydrogen burning while material is still falling onto it, so strong radiation, ionisation and outflows appear while the star is still being assembled. That used to look like a serious obstacle to making massive stars at all, since radiation pressure might simply blow away the infalling gas. The answer is that the accretion flow isn’t spherical: dense disc-fed accretion can continue while radiation and ionised gas escape preferentially through lower-density directions. Observations and simulations now suggest more continuity between low-mass and high-mass star formation than older pictures implied, though massive stars live in denser, more strongly accreting and more strongly clustered environments.
The newly formed stars then start destroying the cloud that made them. Take an O star producing roughly $Q\sim10^{49}$ ionising photons per second. In uniform hydrogen gas, ionisations balance recombinations at the Strömgren radius $$R_S=\left(\frac{3Q}{4\pi n^{2}\alpha_B}\right)^{1/3}\sim3\ \mathrm{pc}$$ for $n=100\ \mathrm{cm^{-3}}$ and a standard case-B recombination coefficient. The gas inside is heated to roughly $10^4$ K, far above the temperature of the molecular material around it, so the H II region expands and drives a shock into the surrounding cloud. Radiation pressure, stellar winds and protostellar outflows add more momentum, and a few million years later the most massive stars begin to explode.
Feedback doesn’t explain the entire low value of $\epsilon_{\rm ff}$ by itself. Turbulence, magnetic fields, cloud geometry and the continual assembly and dispersal of gas matter before the first O star appears. What feedback does effectively is put a clock on the process, because once enough stars have formed, especially massive ones, the cloud begins losing the reservoir from which further stars could have formed. That helps resolve something which otherwise sounds contradictory. Star formation is slow measured against the molecular gas available, while individual dense structures can still collapse quickly. The Galaxy doesn’t need every dense core to collapse at one percent of free fall. It needs only a small fraction of the molecular gas to be in a rapidly collapsing state at any one time, with feedback and dynamics recycling or dispersing the rest.
The timescales then form a rough hierarchy. Giant molecular clouds live for something like several to a few tens of millions of years, depending on environment and on how a cloud is defined. Dense prestellar cores evolve over hundreds of thousands of years. The deeply embedded Class 0 and I phases together last of order half a million years for nearby low-mass populations. Protoplanetary discs commonly survive a few million years with a broad spread. A solar-mass star then takes tens of millions of years to settle onto the main sequence, while a $20\,M_\odot$ object reaches it while still accreting and starts ionising and disturbing its birth environment before the formation process around it has finished.
Put the chain back together and the original twenty-four orders of magnitude stop looking like one collapse. The gas first has to become cold and concentrated enough for self-gravity to beat pressure, which is the Jeans condition. Then the obvious prediction fails, because clouds don’t turn themselves into stars in one free-fall time, and turbulence, magnetic fields, cloud structure and feedback enter. Once a dense core collapses, radiative transfer changes the equation of state and produces the first and second hydrostatic cores. Rotation creates a disc, the disc loses angular momentum through gravitational torques, magnetic stresses and winds, and the forming star launches outflows. Massive stars eventually ionise and disrupt the larger cloud. Somewhere inside all of that, the process also produces nearly the same broad distribution of stellar masses again and again.
Some parts of that chain are on firm ground. The Jeans dispersion relation follows directly from the fluid equations, the free-fall time follows from Newtonian gravity, the need to transport angular momentum is unavoidable, and the first and second collapse follow from the thermodynamics of molecular hydrogen and radiative transfer. Other parts are much less settled. Star formation is inefficient on molecular-cloud scales, and there’s no agreed partition of responsibility between turbulence, magnetic fields, cloud evolution and feedback. Discs move angular momentum, and the relative roles of magnetised winds, MRI, hydrodynamic instabilities and gravitational torques vary from one part of a disc to another. The IMF is measured very well, and its full shape still can’t be derived from first principles.
References and Footnotes
1
There is a genuine inconsistency hidden here, traditionally called the Jeans swindle. An infinite medium with constant density cannot simultaneously satisfy $\nabla^2\Phi_0=4\pi G\rho_0$ and have no background gravitational acceleration. The usual calculation ignores the gravity of the uniform background and keeps the gravity of the perturbation. More careful treatments using finite systems or expanding backgrounds recover essentially the same instability criterion, so the answer survives even though the original setup is not mathematically self-consistent. ↩︎
2
Image credit: ESA/Herschel/PACS, SPIRE/Gould Belt survey Key Programme. The characteristic $0.1$ pc width was emphasised in the Herschel Gould Belt work, but later studies have questioned how universal it really is. The safe statement is that filaments are common and closely connected to dense core formation. The exact distribution of their widths remains an active observational problem. ↩︎
3
The low value of $\epsilon_{\rm ff}$ is one of the central empirical constraints on star formation theory. Modern cloud-lifecycle work also makes clear that individual clouds need not form stars at a perfectly constant rate. A cloud can evolve through quiet and active phases, so an average value near one percent does not mean every cloud converts exactly one percent of its mass every free-fall time. ↩︎
4
Modern radiation-hydrodynamic calculations give a detailed picture of the first and second core stages, but the observational signature of a genuine first hydrostatic core is not unique. Several good candidates exist. It is safer to call them candidates than to say the phase has been securely observed as a class of objects. ↩︎
5
Image credit: NASA, ESA, CSA, STScI; image processing by Joseph DePasquale, Anton M. Koekemoer and Alyssa Pagan. JWST NIRCam image released in October 2022. ↩︎
6
The history is worth noting. Velikhov and Chandrasekhar studied the instability decades before Balbus and Hawley recognised its importance for accretion discs in 1991, after which MRI became almost synonymous with angular momentum transport in disc theory. Non-ideal MHD changed that. Current work increasingly treats magnetic winds as a major and possibly dominant route for angular momentum loss in protoplanetary discs, with the balance between winds, MRI, hydrodynamic instabilities and gravitational torques depending on radius and evolutionary stage. ↩︎
7
Image credit: ALMA (ESO/NAOJ/NRAO), released in 2014 from ALMA’s long-baseline campaign. ↩︎
8
This is not the only proposed explanation for the IMF peak. Earlier arguments often tied the characteristic mass to the density at which gas and dust thermally couple, while modern calculations increasingly emphasise the thermodynamic transition associated with opacity and the first hydrostatic core. The broader point is the same: some piece of non-scale-free microphysics has to enter if the IMF is going to acquire a preferred mass. ↩︎
Almost everyone makes this mistake once, and it’s a good one to make. You learn the Einstein field equations. You learn that in vacuum, with no matter and no energy and no cosmological constant, they collapse down to something very simple: $$R_{\mu\nu} = 0.$$ So you think: $R_{\mu\nu}$ is the curvature, it’s zero, so there’s no curvature out there, spacetime is flat, nothing is happening. Then somebody points out that the entire outside of a black hole is vacuum. Every single point out there satisfies $R_{\mu\nu}=0$ exactly. And things fall into it.
So the reasoning broke somewhere, and I want to find out where, because once you see which step failed you’ll understand what the Ricci tensor actually measures, what the rest of the curvature is doing, and why gravity has waves at all. Short version: $R_{\mu\nu}$ isn’t the curvature. It’s a contraction of the curvature, and contracting throws information away.
Before any tensors show up, though, I want to convince you that you already believe this and have believed it since your first mechanics course. Newtonian gravity obeys Poisson’s equation, $$\nabla^2 \Phi = 4\pi G \rho.$$ Where there’s no matter, $\rho = 0$, and you get Laplace’s equation, $\nabla^2\Phi = 0$. Now, does anybody think Laplace’s equation means the gravitational field is zero? No. The potential outside a spherical mass is $\Phi = -GM/r$, you can check for yourself that $\nabla^2(-GM/r) = 0$ everywhere except the origin, and $\Phi$ is obviously not zero, its gradient isn’t zero, and things fall. The air above your head is vacuum. Laplace’s equation holds up there. You’re still stuck to the floor.
The same puzzle in Newtonian language. Outside the source, Poisson’s equation becomes Laplace’s equation. That sets only the trace of the Hessian $\partial_i\partial_j\Phi$ to zero; it does not set the potential, the field, or the trace-free tidal Hessian to zero.
That last bit is the whole answer already, so let me slow down on it. The second derivatives of the potential form a symmetric matrix, the Hessian $\partial_i \partial_j \Phi$, and in three dimensions that has six independent entries. Laplace’s equation sets exactly one number to zero: the trace, $\sum_i \partial_i\partial_i \Phi$. Five entries are left completely alone. Those five are the tidal field, the thing that stretches you along the line toward the Earth and squeezes you sideways, and they’re what gravity does in empty space.
General relativity does the same thing with better bookkeeping. The full curvature of spacetime is the Riemann tensor $R^\rho{}{\sigma \mu \nu}$, and the Ricci tensor isn’t some separate object, it’s what you get by contracting Riemann on a pair of indices: $$R{\mu\nu} = R^\alpha{}{\mu\alpha\nu}.$$ You’re summing over one direction. That’s an averaging operation, and averaging loses information. If I tell you a list of numbers averages to zero, you don’t conclude the numbers were all zero, because ${+1, -1}$ averages to zero too. So $R{\mu\nu} = 0$ says that certain averages of the curvature vanish. It doesn’t say the curvature vanishes.
You can make that exact by counting. Riemann has a lot of symmetry: antisymmetric in the first two indices, antisymmetric in the last two, and symmetric when you swap the two pairs. In four dimensions an antisymmetric pair takes $\binom{4}{2}=6$ values, so the pair symmetry makes Riemann a symmetric $6\times 6$ matrix with $\tfrac{6\cdot 7}{2}=21$ entries, and the first Bianchi identity $R_{\mu[\nu\rho\sigma]}=0$ knocks off one more. So the Riemann tensor of a four dimensional spacetime has $$20 \ \text{independent components.}$$ The Ricci tensor is symmetric, so it has $\tfrac{4\cdot 5}{2}=10$. Setting $R_{\mu\nu}=0$ puts ten conditions on twenty numbers.1Worth doing this count yourself once. The first Bianchi identity $R_{\mu[\nu\rho\sigma]}=0$ looks like it should impose a lot of conditions, but nearly all of them already follow from the pair symmetries. The only new content is the totally antisymmetric part $R_{[\mu\nu\rho\sigma]}=0$, and in four dimensions there’s exactly $\binom{4}{4}=1$ way to pick four antisymmetrised indices, so exactly one new constraint. In general the answer is $\frac{n^2(n^2-1)}{12}$, which gives 1 for $n=2$, 6 for $n=3$, 20 for $n=4$, and 50 for $n=5$. Ten are left over. Whatever those ten are, they’re the curvature that empty space is allowed to have.
At a point, the vacuum condition removes the ten Ricci-determined components of the four-dimensional Riemann tensor while permitting ten independent Weyl components to remain. “Permitting” does not mean unconstrained dynamics: the Bianchi identities govern how the Weyl curvature propagates through vacuum.
The ten leftovers have a name. There’s exactly one way to split the Riemann tensor into a piece built out of Ricci and a remainder that’s completely traceless, and the remainder is the Weyl tensor $C_{\mu\nu\rho\sigma}$. In $n$ dimensions the split reads
$$
R_{\mu\nu\rho\sigma} = C_{\mu\nu\rho\sigma} + \frac{2}{n-2}\Bigl(g_{\mu[\rho}R_{\sigma]\nu} – g_{\nu[\rho}R_{\sigma]\mu}\Bigr) – \frac{2R}{(n-1)(n-2)}\, g_{\mu[\rho}g_{\sigma]\nu},
$$
where the square brackets mean antisymmetrise. Every trace of $C$ vanishes by construction, $C^\alpha{}{\mu\alpha\nu}=0$, and that’s exactly why the field equations can’t see it. So in vacuum, where $R{\mu\nu}=0$ and therefore $R=0$ as well, the whole thing collapses to
$$
R_{\mu\nu\rho\sigma} = C_{\mu\nu\rho\sigma}.
$$
The curvature of empty space is pure Weyl. Not mostly Weyl, not approximately Weyl. Exactly and only Weyl.2The Weyl tensor has $\frac{n(n+1)(n+2)(n-3)}{12}$ independent components in $n$ dimensions. It’s also conformally invariant: rescale the metric by any positive function, $g_{\mu\nu} \to \Omega^2 g_{\mu\nu}$, and the Weyl tensor with one index raised doesn’t change. That’s why it sometimes gets called the conformal curvature tensor, and why a spacetime is conformally flat exactly when its Weyl tensor vanishes.
Now stare at that component formula for a second, because there’s something odd hiding in it. The factor $(n-3)$ means the Weyl tensor has zero components whenever $n \le 3$. Not zero for some particular metric, zero identically, always, for every three dimensional geometry there is. Which means the Riemann tensor in three dimensions is completely fixed by the Ricci tensor. Which means that in three dimensions, $R_{\mu\nu}=0$ really does force spacetime to be flat.
Independent local curvature data by dimension. In two and three dimensions the Riemann tensor contains no independent Weyl part, so Ricci-flat implies flat. Four dimensions are the first in which curvature can remain after all Ricci curvature has vanished. The 2D Ricci entry counts curvature information, not the components of an arbitrary symmetric rank-two tensor.
This still surprises me. The wrong intuition isn’t wrong because anyone reasoned badly. It’s wrong because it’s reasoning about a universe with one dimension too few. In three dimensional gravity there’s no curvature outside a source, no light bending around a star, no structure in a black hole exterior, and no gravitational radiation, because there’s simply nowhere for that information to sit. The theory is topological. Our fourth dimension is what opens up ten free slots in the curvature and turns gravity into a field with a life of its own.3Three dimensional gravity isn’t boring, and the qualification matters: with a negative cosmological constant you get the BTZ black hole, with a horizon and a temperature and an entropy. But that isn’t a counterexample. BTZ is locally identical to anti de Sitter space everywhere, with constant curvature set entirely by $\Lambda$; what makes it a black hole is a global identification of points, which is topology rather than local geometry. There are still no local gravitational degrees of freedom in three dimensions and no gravitational waves.
So Weyl is what survives. The question that matters is what it does, and here the cleanest statement in the whole subject is waiting. Take a small ball of test particles, freely falling, starting at rest relative to each other. Coffee grounds released in a spacecraft will do fine. Watch the ball. Its volume is controlled by the expansion $\theta = \nabla_\mu u^\mu$ of the congruence of worldlines, which obeys the Raychaudhuri equation $$\frac{d\theta}{d\tau} = -\frac{1}{3}\theta^2 – \sigma_{\mu\nu}\sigma^{\mu\nu} + \omega_{\mu\nu}\omega^{\mu\nu} – R_{\mu\nu}u^\mu u^\nu,$$ with $\sigma_{\mu\nu}$ the shear and $\omega_{\mu\nu}$ the vorticity. Start the ball at rest with no shear and no rotation, so at the first instant $\theta$, $\sigma$ and $\omega$ are all zero, and every term on the right dies except the last one. You’re left with $$\left.\frac{\ddot V}{V}\right|{\tau=0} = -R{\mu\nu}u^\mu u^\nu.$$ The Ricci tensor is the thing, and the only thing, that changes the volume of a ball of freely falling matter.4Getting from $\theta$ to volume: $\theta = \dot V/V$, so $\dot\theta = \ddot V/V – (\dot V/V)^2$, and at the initial instant where $\dot V=0$ that’s just $\ddot V / V$. Feed in Einstein’s equation and the geometric statement becomes a statement about matter: for a perfect fluid you get $\ddot V/V = -4\pi G(\rho + 3p/c^2)$, which is a very direct way to see why pressure gravitates in general relativity, and why no equation of state stiff enough can save a star from collapsing.
Now put the ball in vacuum. There $R_{\mu\nu}u^\mu u^\nu = 0$ identically, so the volume doesn’t change at all. But something obviously does happen to coffee grounds falling toward the Earth, because the ones nearer the Earth fall faster and the ones off to the sides fall along slightly converging lines. The ball deforms. It just deforms at fixed volume. In vacuum the geodesic deviation equation reads $$\frac{D^2 \xi^\mu}{d\tau^2} = -R^\mu{}{\alpha\nu\beta},u^\alpha \xi^\nu u^\beta = -E^\mu{}\nu, \xi^\nu, \qquad E_{\mu\nu} \equiv C_{\mu\alpha\nu\beta}u^\alpha u^\beta,$$ where $E_{\mu\nu}$ is called the electric part of the Weyl tensor. It’s symmetric, and since every trace of Weyl vanishes it’s traceless too: $E^\mu{}\mu = 0$. A traceless matrix of relative accelerations is exactly the statement that volume is preserved while shape isn’t. That’s the division of labour, and it’s a theorem, not a slogan. Ricci changes volume. Weyl changes shape.5Weyl also has a magnetic part, $B{\mu\nu} = \tfrac{1}{2}\epsilon_{\mu\alpha\beta\gamma}C^{\alpha\beta}{}_{\nu\delta}u^\gamma u^\delta$, and the two together account for all ten components, five each, both symmetric and traceless. The electric and magnetic language isn’t a loose analogy either: $E$ and $B$ obey equations closely parallel to Maxwell’s, and it’s the two of them oscillating into each other that makes a gravitational wave. A purely electric Weyl field is just the static tidal field of a mass. You need the magnetic part to get radiation.
For an initially round, comoving ball of freely falling particles, Ricci curvature controls the trace of the relative-acceleration tensor and therefore the initial volume acceleration. In vacuum that trace vanishes, while Weyl curvature can immediately produce shear. Thus $\ddot V(0)=0$; the figure should not be read as claiming exact volume conservation at all later times, because Weyl-generated shear can subsequently enter the Raychaudhuri equation.
You can check the whole story with nothing but calculus, and I’d recommend actually doing it, because it makes the abstraction concrete. Take $\Phi = -GM/r$ and compute the Hessian: $$\partial_i \partial_j \Phi = GM\left(\frac{\delta_{ij}}{r^3} – \frac{3 x_i x_j}{r^5}\right).$$ The trace is $GM(3/r^3 – 3r^2/r^5) = 0$, which is Laplace’s equation turning up as promised. The eigenvalues are easy to read off: along the radial direction you get $-2GM/r^3$, and in each of the two transverse directions you get $+GM/r^3$. Since the relative acceleration of nearby particles is $\ddot\xi^i = -\partial_i\partial_j\Phi\,\xi^j$, the ball gets pulled apart radially at a rate $2GM/r^3$ and squeezed inward in each transverse direction at $GM/r^3$. One stretch, two squeezes, and because $2 = 1+1$ the volume stays put. That’s the traceless condition, sitting there in numbers you can work out in a minute.
The Newtonian tidal field of a point mass. The relative-acceleration eigenvalues are $+2GM/r^3$ radially and $-GM/r^3$ in each of the two transverse directions. Their trace is zero, so an initially comoving small ball has zero initial volume acceleration even though it immediately begins to distort.
Some numbers help here. Near the Earth’s surface the tidal rate is $2GM_\oplus/R_\oplus^3 \approx 3.1\times 10^{-6}\ \mathrm{s^{-2}}$, so across the two metres of a standing person the difference in acceleration between head and feet is about $6\times 10^{-6}\ \mathrm{m\,s^{-2}}$, which is why you’ve never noticed it. Go to a neutron star of $1.4$ solar masses and stand $100$ kilometres away and the same expression gives $3.7\times 10^{5}\ \mathrm{s^{-2}}$, so the head to foot difference is around seventy five thousand times Earth gravity and you get pulled into a thread. Both of those places are perfect vacuum. Both have $R_{\mu\nu}=0$ exactly. The difference between an effect you can’t feel and one that kills you is entirely a difference in the Weyl tensor, which the field equations never mention.
If you want a proof rather than an argument that Schwarzschild isn’t flat, compute a curvature invariant, because scalars can’t be changed by any coordinate transformation. The Kretschmann scalar is $$K = R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma} = \frac{48G^2M^2}{c^4 r^6}.$$ That’s nonzero at every finite radius. If spacetime were flat then $R_{\mu\nu\rho\sigma}=0$ and $K$ would have to vanish, and it doesn’t. There’s no clever choice of coordinates that turns a nonzero scalar into a zero one. So the outside of a star is curved, curved despite being empty, and since Riemann and Weyl coincide there we can also write $K = C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}$.6The same scalar settles the other classic Schwarzschild confusion. At the horizon $r = r_S$ the metric components go haywire, but $K = 48G^2M^2/(c^4 r_S^6)$ is perfectly finite there, which tells you the horizon is a problem with the coordinates and not with the geometry. As $r \to 0$ the scalar blows up, and that singularity is real and no choice of coordinates removes it.
The most extreme version of all this is gravitational radiation. A gravitational wave crossing interstellar space solves $R_{\mu\nu}=0$ at every point along its path. No matter anywhere in it, nothing on the right hand side of the field equations, nothing at all in the naive sense. It’s pure Weyl curvature moving at the speed of light, made entirely out of the components the vacuum equations decline to constrain. And when one arrives it does exactly what we derived above: changes the shape of things while leaving their volume alone, stretching along one axis and squeezing along the perpendicular one. That’s precisely what LIGO measures, a difference in the lengths of two perpendicular arms, and in September 2015 that difference was about $4\times 10^{-18}$ metres over a four kilometre baseline, roughly a thousandth of the width of a proton. Empty space, curved, and curved strongly enough to detect.
So the corrected statement is this. $R_{\mu\nu}=0$ doesn’t say spacetime is flat. It says spacetime is Ricci flat, which is a real and restrictive condition, but one that leaves half the curvature untouched in four dimensions. The half it kills is the half that responds directly to local matter and changes the volumes of falling objects. The half it leaves is the tidal, shape distorting, freely propagating half, and that half isn’t sourced locally at all. It’s set by matter somewhere else and carried through the vacuum by the Bianchi identities, in exactly the way that $\nabla^2\Phi=0$ in the air above your head is perfectly compatible with the Earth underneath pulling on you.7The mechanism by which distant matter fixes Weyl in a vacuum region is the second Bianchi identity, $\nabla_{[\lambda}R_{\mu\nu]\rho\sigma}=0$. Contract it in vacuum and you get $\nabla^\mu C_{\mu\nu\rho\sigma}=0$, a wave equation for the Weyl tensor where the matter enters only through boundary conditions. That’s the precise sense in which Weyl curvature is the free gravitational field, and the parallel with the source free Maxwell equations $\nabla^\mu F_{\mu\nu}=0$ is very close.
One last thing, for anyone drifting toward pure mathematics. Ricci flat manifolds aren’t something physicists invented to describe black holes, they’re a central object in Riemannian geometry, and $R_{\mu\nu}=0$ is a hard nonlinear equation with a big space of solutions. Calabi conjectured and Yau proved that every compact Kähler manifold with vanishing first Chern class carries a Ricci flat metric, and the resulting Calabi-Yau manifolds are highly curved things that happen to have no Ricci curvature anywhere. The K3 surface is a four dimensional example: compact, Ricci flat, and very much not flat. When a geometer says Ricci flat they never mean flat, and it wouldn’t occur to them that anyone might. It only occurs to us because we met the Ricci tensor for the first time in an equation with the word vacuum sitting next to it.
There’s a general lesson in here somewhere, and it’s the thing I’d want you to keep. When a theory hands you an equation of the form “some contraction of the interesting object is zero,” the first question isn’t what the equation forbids rather it’s what it leaves alone. Count the components. Find the part the equation never touches. In general relativity that part is the Weyl tensor, it has ten components, it exists only because we live in four dimensions instead of three, and it’s where black holes and tides and gravitational waves all live.
References and Footnotes
1
Worth doing this count yourself once. The first Bianchi identity $R_{\mu[\nu\rho\sigma]}=0$ looks like it should impose a lot of conditions, but nearly all of them already follow from the pair symmetries. The only new content is the totally antisymmetric part $R_{[\mu\nu\rho\sigma]}=0$, and in four dimensions there’s exactly $\binom{4}{4}=1$ way to pick four antisymmetrised indices, so exactly one new constraint. In general the answer is $\frac{n^2(n^2-1)}{12}$, which gives 1 for $n=2$, 6 for $n=3$, 20 for $n=4$, and 50 for $n=5$. ↩︎
2
The Weyl tensor has $\frac{n(n+1)(n+2)(n-3)}{12}$ independent components in $n$ dimensions. It’s also conformally invariant: rescale the metric by any positive function, $g_{\mu\nu} \to \Omega^2 g_{\mu\nu}$, and the Weyl tensor with one index raised doesn’t change. That’s why it sometimes gets called the conformal curvature tensor, and why a spacetime is conformally flat exactly when its Weyl tensor vanishes. ↩︎
3
Three dimensional gravity isn’t boring, and the qualification matters: with a negative cosmological constant you get the BTZ black hole, with a horizon and a temperature and an entropy. But that isn’t a counterexample. BTZ is locally identical to anti de Sitter space everywhere, with constant curvature set entirely by $\Lambda$; what makes it a black hole is a global identification of points, which is topology rather than local geometry. There are still no local gravitational degrees of freedom in three dimensions and no gravitational waves. ↩︎
4
Getting from $\theta$ to volume: $\theta = \dot V/V$, so $\dot\theta = \ddot V/V – (\dot V/V)^2$, and at the initial instant where $\dot V=0$ that’s just $\ddot V / V$. Feed in Einstein’s equation and the geometric statement becomes a statement about matter: for a perfect fluid you get $\ddot V/V = -4\pi G(\rho + 3p/c^2)$, which is a very direct way to see why pressure gravitates in general relativity, and why no equation of state stiff enough can save a star from collapsing. ↩︎
5
Weyl also has a magnetic part, $B{\mu\nu} = \tfrac{1}{2}\epsilon_{\mu\alpha\beta\gamma}C^{\alpha\beta}{}_{\nu\delta}u^\gamma u^\delta$, and the two together account for all ten components, five each, both symmetric and traceless. The electric and magnetic language isn’t a loose analogy either: $E$ and $B$ obey equations closely parallel to Maxwell’s, and it’s the two of them oscillating into each other that makes a gravitational wave. A purely electric Weyl field is just the static tidal field of a mass. You need the magnetic part to get radiation. ↩︎
6
The same scalar settles the other classic Schwarzschild confusion. At the horizon $r = r_S$ the metric components go haywire, but $K = 48G^2M^2/(c^4 r_S^6)$ is perfectly finite there, which tells you the horizon is a problem with the coordinates and not with the geometry. As $r \to 0$ the scalar blows up, and that singularity is real and no choice of coordinates removes it. ↩︎
7
The mechanism by which distant matter fixes Weyl in a vacuum region is the second Bianchi identity, $\nabla_{[\lambda}R_{\mu\nu]\rho\sigma}=0$. Contract it in vacuum and you get $\nabla^\mu C_{\mu\nu\rho\sigma}=0$, a wave equation for the Weyl tensor where the matter enters only through boundary conditions. That’s the precise sense in which Weyl curvature is the free gravitational field, and the parallel with the source free Maxwell equations $\nabla^\mu F_{\mu\nu}=0$ is very close. ↩︎
Cosmic rays arrive at Earth from every direction almost perfectly evenly, to about one part in a thousand. That’s a disaster if you wanted to know where they came from, because they’re charged and the Galaxy is magnetised, so whatever direction they started in was erased long ago. The usual response is to say fine, we’ll describe them statistically instead, with a diffusion coefficient. Then somebody hands you a formula, $D \approx 10^{28}\ \mathrm{cm^2\,s^{-1}}$ scaling as $E^{1/3}$, and you’re expected to nod.
I don’t want to do that. A diffusion coefficient isn’t a number you look up, it’s something you derive, and the derivation is genuinely beautiful. It starts from the Lorentz force with nothing else assumed, and along the way the resonance condition that everyone quotes falls out as a delta function rather than being asserted, and the mysterious factor of one third in $D = \lambda v / 3$ turns out to be $\int_{-1}^{1}(1-\mu^2)\,d\mu = 4/3$ divided by 8. I want to walk the whole chain, because once you’ve seen it you never have to memorise any of it again. You can reconstruct it.
So here’s the plan. Write the exact equation of motion for a particle in a mean field plus a fluctuation. Read off how the pitch angle changes. Turn that into a diffusion coefficient in pitch angle using a correlation integral. Discover the resonance. Then turn pitch angle diffusion into spatial diffusion by solving a transport equation. Then, and only then, put in numbers.
Step one. A charged particle in a magnetic field obeys
and because the magnetic force is always perpendicular to $\mathbf{v}$, the speed $v$ never changes. Only the direction does. So the natural variable isn’t the velocity, it’s the pitch angle cosine
$$
\mu\equiv\frac{v_z}{v},
$$
measuring how much of the motion is along the mean field, which I’ll put along $\hat{z}$.
Write the field as $\mathbf{B}=B_0\hat{z}+\delta\mathbf{B}$. Since $v$ is fixed,
Now look hard at that last expression, because it already tells you something important. It only involves the transverse components of $\mathbf{B}$. The mean field is purely along $\hat{z}$, so $B_{0x}=B_{0y}=0$ and the mean field contributes nothing at all. Only $\delta\mathbf{B}$ can change the pitch angle.
1This is worth sitting with, because it’s the reason the whole subject exists. A uniform field cannot scatter anything. It rotates the velocity around $\hat{z}$ at a constant rate and leaves the angle to $\hat{z}$ untouched forever, so a particle in a uniform field slides along its field line for eternity with perfect memory of where it was going. Every bit of transport physics that follows comes from the difference between the real field and a uniform one.The unperturbed motion carries the velocity vector around a cone of fixed opening angle, so $\mu$ is a constant of the motion. Nothing in a uniform field can change it. Only the transverse part of $\delta\mathbf B$ appears in $d\mu/dt$, which is why the entire transport problem is a problem about the fluctuating field and not about the mean one.
Let’s make that explicit. Write the gyrofrequency $\Omega = qB_0/\gamma m c$, and use the unperturbed orbit for the velocity, which is legitimate as long as the fluctuation is a small perturbation over one turn. On that orbit $v_x = v_\perp\cos\psi$ and $v_y = -v_\perp\sin\psi$ with phase $\psi = \Omega t + \psi_0$ and $v_\perp = v\sqrt{1-\mu^2}$. Substituting, $$\boxed{\ \frac{d\mu}{dt} = \frac{\Omega\sqrt{1-\mu^2}}{B_0}\Big[\cos(\Omega t + \psi_0)\,\delta B_y + \sin(\Omega t + \psi_0)\,\delta B_x\Big].\ }$$ That single line contains everything. Three features of it will each become a piece of physics later, so let me point them out now. The factor $\sqrt{1-\mu^2}$ vanishes when $\mu = \pm 1$, meaning a particle streaming exactly along the field has no transverse velocity for the fluctuation to push against, so it cannot be scattered. The fluctuation enters only through its transverse components. And the whole kick oscillates at the gyrofrequency, which is going to select which fluctuations matter.
Step two is to turn a fluctuating kick into a diffusion coefficient, and the recipe is completely general. Any quantity being pushed around randomly performs a random walk, and the accumulated change over time $t$ is $\Delta\mu = \int_0^t \dot\mu\,dt’$. That has zero mean, but its square doesn’t: $$\left\langle (\Delta\mu)^2 \right\rangle = \int_0^t!!\int_0^t \left\langle \dot\mu(t’)\,\dot\mu(t”) \right\rangle dt’\,dt”.$$ Now suppose the correlation function $C(\tau) = \langle\dot\mu(0)\dot\mu(\tau)\rangle$ depends only on the time difference and dies away after some correlation time $\tau_c$. Then for $t \gg \tau_c$ the double integral collapses: for each $t’$, the inner integral over $t”$ contributes the same finite amount $\int C(\tau)d\tau$, and there are $t$ worth of choices of $t’$. So $\langle(\Delta\mu)^2\rangle$ grows linearly in $t$, which is exactly what diffusion means, and defining $D_{\mu\mu} = \langle(\Delta\mu)^2\rangle/2t$ gives $$D_{\mu\mu} = \frac{1}{2}\int_{-\infty}^{\infty} \left\langle \dot\mu(0)\,\dot\mu(\tau) \right\rangle d\tau.$$ A transport coefficient is the area under a correlation function. This is a Green-Kubo relation, and the same structure gives you electrical conductivity, viscosity and thermal conductivity in other contexts. The physical content is that transport is set by how long a fluctuation remembers itself.2The condition $t \gg \tau_c$ is not a technicality, it’s the entire definition of diffusive behaviour. On timescales shorter than the correlation time the motion is ballistic and $\langle(\Delta\mu)^2\rangle \propto t^2$, not $t$. Everything in this post applies only after the particle has been scattered many times. That’s also why diffusion fails at the highest cosmic ray energies, where a particle crosses the Galaxy before it gets scattered even once.
The Green-Kubo picture. Take the quantity that’s being kicked around, correlate it with itself at a later time, and integrate. The area is the diffusion coefficient. The same relation gives conductivity, viscosity and every other transport coefficient in physics, which is a good hint that the structure is more fundamental than the application.
Step three is where the resonance appears, and I think it’s the prettiest moment in the whole derivation. Put the expression for $\dot\mu$ into the correlation function. Two things get correlated. The gyration phase gives $\langle\cos\psi(0)\cos\psi(\tau)\rangle$, which after averaging over the initial phase is $\tfrac12\cos\Omega\tau$. The field gives $\langle \delta B(z_0)\,\delta B(z_0 + v\mu\tau)\rangle$, because in time $\tau$ the particle has moved a distance $v\mu\tau$ along the field and is sampling the fluctuation at a new place. Writing that correlation in Fourier space with power spectrum $P(k)$ turns it into $\int dk\, P(k)\,e^{ikv\mu\tau}$. So the thing we have to integrate over all $\tau$ is $$\int_{-\infty}^{\infty} d\tau\ \cos(\Omega\tau)\, e^{i k v\mu \tau} = \pi\Big[\delta(kv\mu – \Omega) + \delta(kv\mu + \Omega)\Big].$$ There it is. A delta function. Not a rule of thumb, not a scaling argument, an exact statement that out of the entire turbulent spectrum the particle feels only the wavenumber satisfying $kv\mu = \Omega$, which rearranges to $$k_{\rm res} = \frac{\Omega}{v\mu} = \frac{1}{r_g \mu}, \qquad r_g \equiv \frac{v}{\Omega} = \frac{pc}{qB_0}.$$ The physical reading is exactly what you’d hope. The particle turns at rate $\Omega$ and drifts along the field at speed $v\mu$. If the wave it’s passing through has a wavelength such that it advances by one wavelength in exactly one gyroperiod, then it meets the same phase of the wave at the same phase of its own orbit every single turn, and the kicks add coherently. Any other wavelength and the phase slips, and over many turns the kicks average to nothing.
Where the delta function comes from. The top curve is the particle’s transverse velocity oscillating at $\Omega$, the middle is the fluctuating field it samples as it drifts along the line, and the bottom is their product, which is what actually pushes the pitch angle. When the two are locked in phase the product never changes sign and the kicks accumulate. When they aren’t, the product oscillates about zero and the time integral kills it. That is the entire content of the resonance condition.
Carrying the delta function through the $k$ integral is now just bookkeeping, and it gives $$D_{\mu\mu} = \frac{\pi}{4}\,\Omega\,(1-\mu^2)\ \mathcal{P}(k_{\rm res}), \qquad \mathcal{P}(k) \equiv \frac{k\,P(k)}{B_0^{2}},$$ where $\mathcal P$ is the fraction of the mean field energy sitting in one logarithmic band of wavenumber around $k$, a dimensionless number.3Prefactors of order unity in this formula are convention dependent, because different books normalise the spectrum differently, use one sided or two sided $k$, and include one or both transverse components in $\delta B^2$. I’m using $\int P(k)\,dk = \delta B^2$ with $k$ running over positive values only. If you compare with a textbook and find a factor of two, that’s why, and it doesn’t change any scaling. It’s worth checking this against intuition before going on. It says the pitch angle randomises at a rate equal to the gyrofrequency multiplied by the fractional turbulent power at resonance. If the turbulence at the resonant scale were as strong as the mean field, $\mathcal P \sim 1$, the pitch angle would be scrambled in a single orbit. That’s as fast as scattering can possibly go, and we’ll come back to it.
Step four is to get from pitch angle diffusion to actual spatial transport, and this is the step almost every treatment skips. It shouldn’t, because the answer is not what you’d naively guess. Write the phase space density $f(z,\mu,t)$ and let it obey streaming along the field plus diffusion in pitch angle, $$\frac{\partial f}{\partial t} + v\mu \frac{\partial f}{\partial z} = \frac{\partial}{\partial \mu}\left( D_{\mu\mu} \frac{\partial f}{\partial \mu} \right).$$ Split $f = f_0(z,t) + f_1(z,\mu,t)$ into a nearly isotropic part and a small anisotropy. To leading order the anisotropy is set by balancing the streaming of $f_0$ against the scattering of $f_1$, $$v\mu \frac{\partial f_0}{\partial z} = \frac{\partial}{\partial\mu}\left( D_{\mu\mu}\frac{\partial f_1}{\partial\mu} \right).$$ Integrate once in $\mu$, fixing the constant by demanding nothing blows up at $\mu = \pm 1$, and you get $$\frac{\partial f_1}{\partial \mu} = -\frac{v(1-\mu^2)}{2 D_{\mu\mu}}\frac{\partial f_0}{\partial z}.$$ The flux along the field is $S = \tfrac{v}{2}\int_{-1}^{1}\mu f_1 \,d\mu$, and integrating that by parts to bring in $\partial f_1/\partial\mu$ gives $$S = -\frac{v^2}{8}\left[ \int_{-1}^{1}\frac{(1-\mu^2)^2}{D_{\mu\mu}(\mu)}\,d\mu \right] \frac{\partial f_0}{\partial z}.$$ Comparing with Fick’s law $S = -D_\parallel \,\partial f_0/\partial z$ identifies $$D_\parallel = \frac{v^2}{8}\int_{-1}^{1} \frac{(1-\mu^2)^2}{D_{\mu\mu}(\mu)}\, d\mu.$$
Stop and look at that, because it’s counterintuitive in an important way. The spatial diffusion coefficient contains $D_{\mu\mu}$ in the denominator. Fast pitch angle scattering means slow spatial transport, which makes sense, but more than that, the integral is dominated by whichever pitch angles scatter worst. It’s a bottleneck. If there’s a value of $\mu$ where $D_{\mu\mu}$ becomes small, that region controls the whole answer no matter how efficiently everything else is being scattered, in the same way that a single closed lane sets the travel time on an otherwise empty motorway. And we already know $D_{\mu\mu}$ has a factor $(1-\mu^2)$, so it vanishes at $\mu = \pm1$, and from the resonance it also carries a factor $|\mu|^{q-1}$, so it vanishes at $\mu = 0$ too. Particles moving nearly perpendicular to the field are the bottleneck, and they’re the hardest ones to scatter because the wavelength they resonate with runs off to infinity.
The pitch angle diffusion coefficient across the range of $\mu$. It dies at $\mu = \pm1$, where a particle aligned with the field has no transverse velocity to push on, and it dies at $\mu = 0$, where a particle moving perpendicular to the field would need an infinitely long wave to resonate with. Because $D_\parallel$ integrates the reciprocal, these minima dominate, which is why so much effort in this field goes into what happens near ninety degrees.
Now we can settle the factor of one third that gets quoted everywhere without justification. Take the simplest case, $D_{\mu\mu} = D_0 (1-\mu^2)$, which is what you get if the resonant power happens to be flat in $\mu$. Then the integral is elementary: $$D_\parallel = \frac{v^2}{8 D_0}\int_{-1}^{1}\left(1-\mu^2\right)d\mu = \frac{v^2}{8D_0}\cdot\frac{4}{3} = \frac{v^2}{6 D_0}.$$ If we define a scattering rate $\nu = 2D_0$, the rate at which $\mu$ loses memory, and a mean free path $\lambda = v/\nu$, this is $$D_\parallel = \frac{v^2}{3\nu} = \frac{1}{3}\lambda v.$$ So the one third isn’t a three dimensional geometry factor and isn’t a hand wave. It’s $\tfrac{4}{3} \div 8 \times 2$. The number came out of an integral over pitch angle, and if $D_{\mu\mu}$ had a different shape you’d get a different number.
Everything so far is exact given the model, and none of it needed a spectrum. Now we choose one. Take a power law $P(k) \propto k^{-q}$ from an outer scale $L$ downwards, normalised so the total power is $\delta B^2$, which fixes $$\mathcal{P}(k) = (q-1)\frac{\delta B^2}{B_0^2}\left( kL \right)^{1-q} \quad\Longrightarrow\quad \mathcal{P}(k_{\rm res}) = (q-1)\frac{\delta B^2}{B_0^2}\left( \frac{r_g|\mu|}{L} \right)^{q-1}.$$ Feed that into $D_{\mu\mu}$, feed that into the $D_\parallel$ integral, and the $\mu$ dependence collects into a single number: $$D_\parallel = \frac{I(q)}{2\pi (q-1)}\ v\, r_g \left(\frac{B_0}{\delta B}\right)^{2}\left(\frac{L}{r_g}\right)^{q-1}, \qquad I(q) = \int_{-1}^{1}(1-\mu^2)|\mu|^{1-q}d\mu = \frac{2}{2-q}-\frac{2}{4-q}.$$ For a Kolmogorov cascade, $q = 5/3$, that gives $I = 36/7$ and a prefactor of $1.23$. Notice that $I(q)$ blows up as $q \to 2$: the bottleneck at $\mu = 0$ becomes fatal and the theory stops giving a finite answer, which is the ninety degree problem showing up as a divergent integral rather than as a footnote.
The energy dependence now falls out with no extra input. Since $r_g \propto E$ for relativistic particles, $$D_\parallel \propto r_g \cdot r_g^{-(q-1)} = r_g^{2-q} \propto E^{2-q},$$ so Kolmogorov gives $E^{1/3}$ and a Kraichnan cascade with $q = 3/2$ gives $E^{1/2}$. Measurements of the boron to carbon ratio, which tell you how much interstellar gas cosmic rays have ploughed through and therefore how long they were confined, want an exponent between about $0.3$ and $0.6$. Both cascades sit inside that, which is honest to report and mildly annoying, because it means the data cannot yet decide which cascade the interstellar medium actually has.
Now, finally, numbers, and I want to stress that these are a check on the derivation rather than facts to carry around. Take a $10$ GeV proton in a mean field of $3\ \mu$G, turbulence as strong as the mean field so $\delta B = B_0$, an outer scale of $L = 100$ pc which is roughly the size of a supernova driven eddy, and Kolmogorov. The gyroradius from $r_g = pc/(qB_0)$ is $1.11\times10^{13}$ cm, about three quarters of an astronomical unit. Then $$D_\parallel = 1.23\, c\, r_g \left(\frac{L}{r_g}\right)^{2/3} \approx 3.8\times 10^{28}\ \mathrm{cm^2\,s^{-1}},$$ with a corresponding mean free path $\lambda = 3D/c \approx 1.2$ pc. The value inferred from cosmic ray composition is a few times $10^{28}\ \mathrm{cm^2\,s^{-1}}$. Three inputs you could look up in an afternoon, one derivation, and the answer lands on the measurement. Note also that $\lambda/r_g \approx 3\times10^5$, so the particle really does complete hundreds of thousands of clean orbits between scatterings, which is exactly the condition that made the perturbative treatment legitimate in the first place. The theory checks its own assumptions.
From here the consequences follow in a line. Confinement in a magnetised halo of half thickness $H \approx 4$ kpc takes $\tau \sim H^2/2D \approx 64$ Myr, and radioactive beryllium-10, with a $1.4$ Myr half life, independently reads tens of millions of years from completely unrelated nuclear physics. In that time the particle covers a path length $c\tau \approx 20$ Mpc while getting $4$ kpc from where it started, a ratio of about five thousand. It crossed the Local Group many times over and never left the Galaxy.
I’ve been writing $D_\parallel$ throughout, and the subscript matters. Scattering randomises motion along the field efficiently, but to move across field lines a particle has to hop between lines, which is much harder. So diffusion is anisotropic, and the diffusion coefficient is not a number but a rank two tensor $$D_{ij} = D_\perp \delta_{ij} + \left(D_\parallel – D_\perp\right) b_i b_j, \qquad \mathbf{b} = \mathbf{B}/|\mathbf{B}|.$$ Contract it with $\mathbf b$ twice and you recover $D_\parallel$; contract it with anything perpendicular and you get $D_\perp$. Which is to say the field direction picks out the principal axes and the tensor is diagonal in that frame with eigenvalues $(D_\parallel, D_\perp, D_\perp)$. Classical transport gives $D_\perp/D_\parallel \approx [1 + (\lambda/r_g)^2]^{-1}$, and with our $\lambda/r_g \approx 3\times10^5$ that would be $10^{-11}$, absurdly small.4It is absurdly small, and it’s wrong, for an instructive reason: that estimate assumes the field lines themselves are straight. They’re not. Turbulent field lines separate from each other exponentially, so a particle sliding along one is carried sideways for free without ever hopping. Including this field line random walk puts the ratio nearer $10^{-2}$. Perpendicular transport remains the least settled part of the theory and the place where numerical simulations have contributed most.
Because the field breaks the isotropy of space, the diffusion coefficient has to be a tensor, and the field direction is one of its principal axes. Drawing the surface of constant diffusion gives a cigar pointing along the field line. This is the same structure as the polarizability of a crystal: a direction picked out by the medium, and a tensor that becomes diagonal when you use it.
Strong fields are where the derivation earns its keep, because you can see immediately what happens rather than having to look anything up. Every appearance of $B_0$ is inside $r_g = pc/qB_0$, so raising the field shrinks the gyroradius, shrinks the resonant wavelength, and shrinks the mean free path. But there’s a floor. Look back at $D_{\mu\mu} = \tfrac{\pi}{4}\Omega(1-\mu^2)\mathcal P$: the largest $\mathcal{P}$ can be is of order one, which puts the scattering rate at the gyrofrequency, the mean free path at $r_g$, and $$D_{\rm Bohm} = \tfrac{1}{3} r_g v.$$ You cannot scatter a particle in less than one orbit. Bohm diffusion is the tightest confinement physics allows, and it isn’t an extra assumption, it’s the ceiling on $\mathcal P$.
That floor is what makes supernova remnants work as cosmic ray sources, and the calculation is one line. In diffusive shock acceleration a particle gains energy crossing and recrossing the shock, taking a time $t_{\rm acc} \approx 20 D(E)/u_s^2$ for a shock of speed $u_s$. Take $u_s = 5000$ km s$^{-1}$ and ask how long to reach the knee at $10^{15}$ eV. In the ambient $3\ \mu$G field, $r_g = 1.1\times10^{18}$ cm and even Bohm diffusion gives $D = 1.1\times10^{28}\ \mathrm{cm^2\,s^{-1}}$ and $t_{\rm acc} \approx 28{,}000$ years, which is far longer than the remnant stays fast. But cosmic rays streaming ahead of the shock amplify the field they’re moving through, and at $100\ \mu$G the gyroradius drops by a factor of thirty, so $D$ drops by thirty, and $$t_{\rm acc} \approx 840\ \mathrm{years},$$ comfortably inside a young remnant’s life. The entire case for supernovae as the origin of Galactic cosmic rays rests on that factor of thirty, and you can see exactly where it enters: through $r_g$, through $\mathcal P$, through $D$.5The amplification is the non-resonant streaming instability worked out by Bell in 2004. The cosmic ray current ahead of the shock drives a return current in the background plasma, and the resulting force is unstable to transverse modes on scales much shorter than the cosmic ray gyroradius. Thin non-thermal X-ray filaments in remnants like Cassiopeia A independently point to fields of hundreds of microgauss.
Push the field far higher and the derivation tells you it stops applying, which is more useful than a formula that keeps returning numbers. Near a magnetar at $10^{15}$ G a 1 GeV proton has $r_g = 3\times10^{-9}$ cm, smaller than an atom. There is no turbulence at that scale, so $\mathcal{P}(k_{\rm res})$ is essentially zero, so there is no scattering and no diffusion. The particle is welded to one field line and the problem becomes one dimensional motion along a prescribed curve. Diffusion needs the gyroradius to sit comfortably inside the turbulent inertial range, and outside that window the diffusion coefficient stops meaning anything at all.
The same thing happens at the other end, where it’s more famous. As energy climbs, $r_g$ grows until it’s comparable to the system, at which point $t \gg \tau_c$ fails and the particle simply leaves. Requiring $r_g > R$ for a source of size $R$ and field $B$ gives the Hillas criterion, $E_{\max} \approx 0.9\,Z\beta\,(B/\mu\mathrm{G})(R/\mathrm{kpc})$ EeV. The whole Galactic disk manages a few times $10^{18}$ eV. We detect particles at $10^{20}$ eV with gyroradii of about a hundred kiloparsecs, larger than the Galaxy, so they were neither made here nor confined here, and they arrive barely deflected. The highest energy cosmic rays are the only ones that still remember their direction, which is exactly why so much effort goes into catching the rarest particles in the sky.
Let me be honest about what’s shaky, because a derivation you can’t criticise is a derivation you don’t understand. Quasi-linear theory assumes the fluctuation is a small perturbation over one orbit, and we then applied it with $\delta B \approx B_0$, which is outside its formal domain. Real magnetohydrodynamic turbulence is not the isotropic slab we assumed either; in the Goldreich and Sridhar picture the eddies are stretched along the local field, which starves the resonant parallel wavenumber and predicts much weaker scattering than we computed. And cosmic rays are not test particles. They carry an energy density of about $1$ eV cm$^{-3}$, essentially equal to the magnetic field and the turbulent gas, so they generate the very waves that scatter them and the problem is really nonlinear. Every one of those objections attaches to a specific line in the derivation above, which is the point of having done it line by line.
What I’d want you to take from this isn’t the value of $D$. It’s the shape of the argument, because it recurs everywhere. Write the exact microscopic equation of motion. Identify the variable that random walks. Integrate a correlation function to get a transport coefficient. Discover which piece of the environment the system actually couples to, and let a delta function tell you rather than guessing. Then close the loop by asking which part of the phase space is the bottleneck, because that’s what the answer will be controlled by. You can run that program on cosmic rays in the Galaxy, on electrons in a metal, on momentum in a fluid, and on heat in a solid, and it is the same program every time. The particles change and the correlation function changes. Nothing else does.
References and Footnotes
1
This is worth sitting with, because it’s the reason the whole subject exists. A uniform field cannot scatter anything. It rotates the velocity around $\hat{z}$ at a constant rate and leaves the angle to $\hat{z}$ untouched forever, so a particle in a uniform field slides along its field line for eternity with perfect memory of where it was going. Every bit of transport physics that follows comes from the difference between the real field and a uniform one. ↩︎
2
The condition $t \gg \tau_c$ is not a technicality, it’s the entire definition of diffusive behaviour. On timescales shorter than the correlation time the motion is ballistic and $\langle(\Delta\mu)^2\rangle \propto t^2$, not $t$. Everything in this post applies only after the particle has been scattered many times. That’s also why diffusion fails at the highest cosmic ray energies, where a particle crosses the Galaxy before it gets scattered even once. ↩︎
3
Prefactors of order unity in this formula are convention dependent, because different books normalise the spectrum differently, use one sided or two sided $k$, and include one or both transverse components in $\delta B^2$. I’m using $\int P(k)\,dk = \delta B^2$ with $k$ running over positive values only. If you compare with a textbook and find a factor of two, that’s why, and it doesn’t change any scaling. ↩︎
4
It is absurdly small, and it’s wrong, for an instructive reason: that estimate assumes the field lines themselves are straight. They’re not. Turbulent field lines separate from each other exponentially, so a particle sliding along one is carried sideways for free without ever hopping. Including this field line random walk puts the ratio nearer $10^{-2}$. Perpendicular transport remains the least settled part of the theory and the place where numerical simulations have contributed most. ↩︎
5
The amplification is the non-resonant streaming instability worked out by Bell in 2004. The cosmic ray current ahead of the shock drives a return current in the background plasma, and the resulting force is unstable to transverse modes on scales much shorter than the cosmic ray gyroradius. Thin non-thermal X-ray filaments in remnants like Cassiopeia A independently point to fields of hundreds of microgauss. ↩︎
You have seen this demonstration. Somebody stretches a rubber sheet over a frame, drops a bowling ball in the middle, and the sheet sags. Then they roll a marble across it and the marble curves toward the bowling ball, maybe even loops around it once or twice before spiralling in. And the narrator says: that is gravity. Mass curves spacetime, and curved spacetime tells matter how to move. Everyone nods. It is on the cover of textbooks. It is in every documentary. I have used it myself when someone asks me at a party what general relativity is about, because it takes ten seconds and it gets the shape of the idea across.
But it is wrong. Not wrong in the sense of being a simplification that leaves out details, which would be fine, but wrong in a specific and interesting way: it is an accurate picture of something that is not the thing causing gravity. And once you see exactly what it has left out, you understand general relativity considerably better than the demonstration was ever going to teach you.
The objection you usually hear is that the demonstration is circular. The marble rolls toward the bowling ball because real gravity is pulling it down into the dip, so you are using gravity to explain gravity. This is true, and it is the weakest objection available. It is an inelegance in the staging, not an error in the physics. You could imagine running the whole thing in free fall with the sheet under tension and a marble constrained to its surface, and you would still have a geometry, and the marble would still follow the shortest path across it. The circularity is a presentational problem. I want to make a harder objection.
Here it is. Take an apple and hold it at rest above the ground. Let go. It falls. Now go back to the rubber sheet and put a marble on it, at rest, on a flat part of the sheet far from the bowling ball. It does not move. To make the marble move on the sheet, you have to give it a velocity, or you have to put it where the sheet is already sloped, at which point real gravity does the work. But the apple was at rest. It had no velocity. It was not sitting on a slope in space. So what curved geometry acted on an object that was not moving?
The rubber sheet cannot answer this, and the reason is that it has deleted the axis along which the apple was moving. An apple at rest in space is not at rest in spacetime. It is moving through the time direction, and moving through it fast. That motion is what the curvature acts on. The rubber sheet shows you two space directions and no time direction at all, which means it has thrown away the entire mechanism.1Sometimes people defend the demonstration by saying the sheet is meant to represent spacetime, not space, with time somehow implicit. It is not. The surface drawn in every version of this demonstration is a two dimensional surface, and both of its dimensions are spatial. There is no axis on the sheet along which a stationary object advances. If there were, a stationary marble would move along it, which is exactly what does not happen.
The rubber sheet is a two dimensional spatial surface. Every direction it shows is a direction of space. The direction along which a stationary apple is actually travelling, and along which the curvature acts, has been removed from the picture before the demonstration begins.
Let me make that precise, because “the time direction matters more” is the kind of claim that should come with an equation attached. Take a weak, static gravitational field, meaning a Newtonian potential $\Phi$ with $|\Phi|/c^2 \ll 1$, which covers the Earth, the Sun, and essentially every situation in which anyone has ever demonstrated a rubber sheet. To first order in $\Phi/c^2$, the metric of spacetime is $$ds^2 = -\left(1 + \frac{2\Phi}{c^2}\right)c^2 dt^2 + \left(1 – \frac{2\Phi}{c^2}\right)\left(dx^2 + dy^2 + dz^2\right).$$ There are two corrections to flat spacetime here. One sits on the time part, $g_{00}$, and one sits on the space part, $g_{ij}$. They are the same size, both of order $\Phi/c^2$. The rubber sheet is a picture of the second one. Let us find out what the first one does.2This form of the metric is written in what is called the Newtonian gauge or the conformal Newtonian gauge. The two potentials appearing in $g_{00}$ and $g_{ij}$ are in general independent functions, often written $\Phi$ and $\Psi$, and they are only forced to be equal when the stress energy tensor has no anisotropic stress. For ordinary matter, dust and stars and planets, this condition holds, so a single $\Phi$ does the job.
A free particle in general relativity moves along a geodesic, which is the curve satisfying $$\frac{d^2 x^\mu}{d\tau^2} + \Gamma^\mu{}_{\alpha\beta}\, \frac{dx^\alpha}{d\tau}\frac{dx^\beta}{d\tau} = 0,$$ where $\tau$ is proper time along the path. This is just the statement that the particle goes as straight as the geometry allows. Now take our apple, which is moving slowly, so that its spatial velocity is negligible compared to $c$. Then $dx^i/d\tau \approx 0$ while $dx^0/d\tau = c\, dt/d\tau \approx c$, and the double sum over $\alpha$ and $\beta$ collapses to the single term with $\alpha = \beta = 0$: $$\frac{d^2 x^i}{d\tau^2} \approx -\Gamma^i{}_{00}\, c^2.$$ Everything now depends on one Christoffel symbol. For a static metric, $$\Gamma^i{}_{00} = -\tfrac{1}{2} g^{ij} \partial_j g_{00},$$ and substituting $g_{00} = -(1 + 2\Phi/c^2)$ and $g^{ij} \approx \delta^{ij}$ gives $\Gamma^i{}_{00} = \partial_i \Phi / c^2$. Therefore $$\frac{d^2 x^i}{dt^2} = -\partial_i \Phi, \qquad \text{that is,} \qquad \vec{a} = -\nabla\Phi.$$ That is Newton’s law of gravitation, recovered exactly, and I want you to look carefully at where it came from. It came from $g_{00}$. The spatial part of the metric, the part the rubber sheet is drawing, never appeared in the calculation. Not once. It contributed nothing.
You can see why from the structure of the geodesic equation. The spatial components of the metric couple to $dx^i/d\tau$, the rate at which the particle moves through space. For an apple, that rate is tiny. The time component couples to $dx^0/d\tau$, the rate at which the particle moves through time, and that rate is enormous. The two corrections to the metric are the same size, but they get multiplied by wildly different things. The relative contribution of the spatial curvature to the motion of a particle moving at speed $v$ is suppressed by a factor of order $(v/c)^2$. An apple that has fallen a metre is moving at about $4.4$ metres per second, so $(v/c)^2 \approx 2 \times 10^{-16}$. The rubber sheet is showing you the part of the effect that contributes roughly one part in ten thousand million million.
If you want independent confirmation that this accounting is right, look at the bending of starlight, which is the one case where the spatial part genuinely matters. A photon has $v = c$, so the suppression factor is gone and both parts of the metric contribute equally. Compute the deflection using only the time curvature, which is what you get from a naive Newtonian calculation treating light as a fast particle, and you find $\Delta\theta = 2GM/(c^2 b)$ for impact parameter $b$. Compute it in full general relativity and you find $$\Delta\theta = \frac{4GM}{c^2 b},$$ exactly twice as much. That missing factor of two is the spatial curvature, contributing its equal half. For the Sun this is $1.75$ arcseconds against a Newtonian $0.87$, and Eddington’s 1919 expedition measured the larger number, which is what made Einstein famous.3The 1919 measurement was considerably less precise than the confident retellings suggest, with error bars large enough that the two candidate values were not separated as cleanly as one would like. The result has since been confirmed to far better precision by radio interferometry of quasars occulted by the Sun, which now agrees with the general relativistic prediction to better than one part in ten thousand. The factor of two is not in doubt. So the spatial curvature is real, and it is measurable, and for light it accounts for half of everything. For an apple it accounts for essentially nothing, and an apple is what the demonstration claims to be explaining.
None of this means the rubber sheet is a fantasy. It is a picture of a real geometric object, and here it is worth being precise about which one. Take the Schwarzschild solution, freeze the time coordinate at some value $t = \text{const}$, and look at the equatorial plane $\theta = \pi/2$. The geometry of that two dimensional slice is $$dl^2 = \left(1 – \frac{r_S}{r}\right)^{-1} dr^2 + r^2 d\phi^2,$$ and if you ask what surface embedded in ordinary flat three dimensional space has this same intrinsic geometry, you get $z(r) = 2\sqrt{r_S(r – r_S)}$, a paraboloid opening outward. That is the funnel. That is the rubber sheet, and it is not an approximation or an artist’s impression: it is the exact embedding diagram of a spatial slice of Schwarzschild spacetime, computed by Flamm in 1916. The demonstration is drawing a correct picture. It is a correct picture of the wrong thing.
So what is the right thing? The curvature that makes apples fall lives in $g_{00}$, and $g_{00}$ is the thing that tells you how fast clocks run. Write out the proper time for a stationary observer at position $x$: setting $dx = dy = dz = 0$ in the metric gives $$d\tau = \sqrt{1 + \frac{2\Phi}{c^2}}\; dt \approx \left(1 + \frac{\Phi}{c^2}\right) dt.$$ A clock deeper in the potential well, where $\Phi$ is more negative, ticks slower. Near the Earth’s surface, $\Phi = gh$ for height $h$, and the fractional difference in tick rate between two clocks separated vertically by $h$ is $gh/c^2$. This is not a metaphor and it is not a small print correction. Pound and Rebka measured it in 1959 using a $22.5$ metre tower at Harvard, where the predicted fractional shift is $$\frac{gh}{c^2} = \frac{9.81 \times 22.5}{8.99 \times 10^{16}} = 2.46 \times 10^{-15},$$ and they confirmed it to about one percent. Fifty years later, optical lattice clocks at NIST resolved the same effect over a height difference of $33$ centimetres. Your head ages faster than your feet by about half a microsecond over an eighty year life, and that is a real, measurable, in principle observable fact about you.4The half microsecond figure comes from $gh/c^2 \approx 1.9 \times 10^{-16}$ for a height difference of $1.7$ metres, multiplied by roughly $2.5 \times 10^9$ seconds in eighty years. Whether this is worth mentioning at parties depends heavily on the party.
The curvature that makes things fall is curvature of the time direction, and it shows up as clocks running at different rates at different heights. This is measured routinely and to high precision. It is also, unlike the sag of a rubber sheet, the actual cause of gravity.
Now, why should a gradient in clock rate make anything move? Here is the principle, and it is one of the most beautiful statements in physics: a free particle follows the path through spacetime that maximises the proper time elapsed along it. Not minimises. Maximises. Of all the ways of getting from one event to another, an unforced particle takes the one on which its own watch reads the longest. Let me show you that this single statement contains Newtonian mechanics whole. For a slow particle in a weak field, $$d\tau = dt\sqrt{1 + \frac{2\Phi}{c^2} – \frac{v^2}{c^2}} \approx dt\left(1 + \frac{\Phi}{c^2} – \frac{v^2}{2c^2}\right),$$ so the total proper time along a path is $$\int d\tau = \int dt + \frac{1}{c^2}\int\left(\Phi – \frac{v^2}{2}\right) dt.$$ The first term is the same for every path between the same two events, so maximising $\int d\tau$ means maximising the second integral, which means minimising $$\int \left(\frac{v^2}{2} – \Phi\right) dt.$$ That integrand is the kinetic energy minus the potential energy, per unit mass. It is the Lagrangian. Extremising proper time is the principle of least action, and running it through the Euler-Lagrange equation gives back $\ddot{x}^i = -\partial_i \Phi$, which is where we came in.5The sign flip between maximising proper time and minimising the action is not a slip. It comes from the overall minus sign in the timelike part of the metric signature, and it means that free particles in relativity extremise a quantity that reduces to the classical action with the conventional sign.
So here is the picture that actually explains gravity. Everything in the universe is moving through spacetime at speed $c$, always. An object sitting still on your desk is not stationary; it is travelling through the time direction at the full speed, and through the space directions at nothing. Gravity does not push on it. What gravity does is make the time direction run at different rates in different places, and a path that keeps going straight through a region where time runs at a gradient does not stay parallel to where it started. Some of the object’s motion gets rotated out of the time direction and into a space direction. That rotation is what you call falling.
And that gives us the analogy I actually want, which is not a sheet at all. It is a cart with two wheels on a common axle, rolling forward. The cart has no steering. Nobody pushes it sideways. If both wheels turn at the same rate, it goes dead straight. Now let the ground under the right wheel be slightly slower, so that the right wheel covers a little less distance per turn than the left one. The cart veers to the right. No sideways force acted on it. The cart was going straight the entire time. It curved because “straight” in a medium with a speed gradient means “bending toward the slow side.”
A cart with no steering, rolling forward. On uniform ground it goes straight. Where the ground under one wheel is slower, it turns toward the slow side without any sideways force acting on it. Replace “distance rolled” with “proper time elapsed” and this is exactly what gravity does to a worldline.
The translation is direct. The cart’s forward motion is the object’s motion through time, which never stops and never slows. The speed gradient across the axle is the gradient in clock rate, $gh/c^2$. The turning is the worldline bending out of the pure time direction and acquiring a spatial component, which is to say the object accelerating downward. And notice what this analogy explains that the rubber sheet cannot. It explains why a stationary object falls, because the cart was never stationary, it was rolling forward the whole time. It explains why gravity is universal and independent of mass or composition, because the turning depends only on the ground, not on what the cart is made of or how heavy it is. And it needs no external gravity to run, because nothing in the setup falls.
It is worth getting a feel for how gentle this bending is. Take the apple falling under $g = 9.81$ and draw its worldline on a spacetime diagram with the time axis measured in metres, so that one second of elapsed time is $c \times 1\,\text{s} = 3 \times 10^8$ metres. In that second the apple moves $4.9$ metres sideways. Its worldline is a parabola $x = gT^2/2c^2$ where $T = ct$, and the radius of curvature of that parabola at its vertex is $$R = \frac{c^2}{g} = \frac{8.99 \times 10^{16}}{9.81} \approx 9.2 \times 10^{15}\,\text{m},$$ which is $0.97$ light years. The worldline of a falling apple is bent on a scale of about one light year. It is almost perfectly straight. Gravity feels overwhelming to us not because the curvature is large but because we persist through enormous stretches of the time direction, and a tiny angle applied over a light year of lever arm still knocks you off a ladder.
A falling apple’s worldline, drawn honestly on the left and exaggerated on the right. In one second the apple advances $3\times10^8$ metres through time and about five metres through space. The curvature is real, and its radius is roughly one light year.
Now let me be honest about where my cart analogy breaks, because an analogy you cannot break is an analogy you do not understand. A uniform speed gradient across the axle is not curvature. It is the flat spacetime of an accelerating observer, which is what the equivalence principle tells you a uniform gravitational field is indistinguishable from. Genuine curvature is what you get when the gradient itself varies from place to place, so that two carts started off parallel do not merely both turn but turn by different amounts and converge. That is the real signature, and it is called geodesic deviation. If $\xi^\mu$ is the separation between two neighbouring geodesics with four velocity $u^\mu$, then $$\frac{D^2 \xi^\mu}{d\tau^2} = -R^\mu{}_{\alpha\nu\beta}\, u^\alpha \xi^\nu u^\beta,$$ where $R^\mu{}_{\alpha\nu\beta}$ is the Riemann tensor. This equation is the honest definition of gravity. In the weak field limit its leading component reduces to $R^i{}_{0j0} \approx \partial_i \partial_j \Phi / c^2$, the Newtonian tidal tensor, which is the thing that stretches you along the radial direction and squeezes you sideways as you fall toward a planet. A uniform field can be transformed away by jumping. Tides cannot, and that is how you know the curvature is real.
Two objects released from rest side by side, drawn in spacetime with time running upward. Both worldlines are geodesics, and yet they converge. In flat spacetime, straight lines that start parallel stay parallel. The convergence is the curvature, and it is measured by the Riemann tensor. This is the part no uniform gradient can imitate.
So the summary is this. The rubber sheet draws a real object, the spatial geometry of a slice of Schwarzschild spacetime, and it draws it accurately. But that object contributes essentially nothing to why things fall, suppressed relative to the effect that matters by a factor of order $(v/c)^2$, which for anything slower than light is a very small number indeed. The gravity you feel comes from the curvature of the time direction, from the plain fact that clocks run at different rates at different heights, and from the principle that a free object takes the path through spacetime along which its own clock reads the most. Nothing pulls the apple. The apple goes straight. Straight, in a place where time runs at a gradient, means down.
I still use the rubber sheet at parties. It takes ten seconds and it gets across the idea that geometry is doing the work, which is the part most people have never heard. But if the conversation lasts longer than ten seconds, I put the sheet away and start talking about clocks, because that is where the physics actually is. The demonstration answers the question “what does curved space look like.” The question we asked was “why do things fall,” and those turn out to be almost entirely different questions.
References and Footnotes
1
Sometimes people defend the demonstration by saying the sheet is meant to represent spacetime, not space, with time somehow implicit. It is not. The surface drawn in every version of this demonstration is a two dimensional surface, and both of its dimensions are spatial. There is no axis on the sheet along which a stationary object advances. If there were, a stationary marble would move along it, which is exactly what does not happen. ↩︎
2
This form of the metric is written in what is called the Newtonian gauge or the conformal Newtonian gauge. The two potentials appearing in $g_{00}$ and $g_{ij}$ are in general independent functions, often written $\Phi$ and $\Psi$, and they are only forced to be equal when the stress energy tensor has no anisotropic stress. For ordinary matter, dust and stars and planets, this condition holds, so a single $\Phi$ does the job. ↩︎
3
The 1919 measurement was considerably less precise than the confident retellings suggest, with error bars large enough that the two candidate values were not separated as cleanly as one would like. The result has since been confirmed to far better precision by radio interferometry of quasars occulted by the Sun, which now agrees with the general relativistic prediction to better than one part in ten thousand. The factor of two is not in doubt. ↩︎
4
The half microsecond figure comes from $gh/c^2 \approx 1.9 \times 10^{-16}$ for a height difference of $1.7$ metres, multiplied by roughly $2.5 \times 10^9$ seconds in eighty years. Whether this is worth mentioning at parties depends heavily on the party. ↩︎
5
The sign flip between maximising proper time and minimising the action is not a slip. It comes from the overall minus sign in the timelike part of the metric signature, and it means that free particles in relativity extremise a quantity that reduces to the classical action with the conventional sign. ↩︎
Mary Somerville is so difficult to put into words that it’s funny. To be honest, who was she? Although you would be correct in calling her a mathematician, you would then be required to clarify the astronomy. The word “astronomer” sounds a bit tame when you consider all the physics, geography, and geology involved in her work. This is no longer an issue because the term “scientist” adequately describes those who spend their lives seeking to comprehend the natural world. The only catch is that there were no scientists when Mary Somerville was born. Of course, people were conducting scientific experiments; it’s just that the term “science” didn’t exist, which is a little strange.
Although she was born in Scotland in 1780, Somerville’s upbringing was far from ideal for a mathematician. Her family initially saw her growing interest in mathematics as a mix of a bothersome curiosity and a reason for concern, likely stemming from their belief that girls weren’t supposed to excel in math. Nonetheless, she gained knowledge. Much of her mathematical education was something she had to seek out on her own; there’s a slightly embellished tale about how she studied geometry and algebra in the middle of the night, but the essential idea remains the same. This is my favorite part of her story because it describes the math learning process so realistically. Some symbols appear, and you find yourself confused. This annoys you. You come across an explanation for them. There are five new things in that explanation that you don’t get. Somehow, you’ve picked up a completely new subject a few months down the road. It appears pointless that Somerville accomplished this without Google.
Over time, she excelled at it. Competent enough to tackle Laplace’s Mécanique Céleste, conversant enough to correspond with prominent European scientists, and capable of delving deeply into astronomical mathematics. The mathematical machinery of the heavens, including the gravitational motion of the planets, was something Laplace was attempting to describe. Famously, his writing was not something you could read on a leisurely basis before turning in for the night. The Mechanism of the Heavens was written by Somerville. This, which is sometimes called a translation, is disrespectful to her writing. While it’s true that she was translating the French, she was also filling in the math and reasoning so that others could learn from it. Cambridge decided to use the book because of how popular it was. As someone who was advised against studying algebra in the past, that path seems rather peculiar.
However, the truly peculiar historical accident was brought about by a book. She released On the Connexion of the Physical Sciences in 1834. Her priorities are clear from the title. Not physics or math, but rather the fields that form their connections. Physics, chemistry, physics, heat, light, electricity, and magnetism. These were getting into increasingly complex topics, but Somerville hoped to demonstrate that they were not independent worlds. They were just two perspectives on the same cosmos. Remember that the best ideas don’t always sound obvious until someone has spent years pointing them out; this may sound obvious when you say it now, but it’s worth remembering.
Because it is necessary, we continue to divide information into different categories. Someone has to make a call at some point: schools have schedules, universities have departments, to designate which classroom will be used for physics and which for mathematics. Nature was never sluggish with this setup. All these things, chemistry, electromagnetism, gravity, and nuclear physics, are happening simultaneously inside a star. It is not a chemistry problem when the chemist walks in. Geometry, thermodynamics, quantum physics, and information theory can all be occupied for an afternoon by contemplation of a black hole. When you delve deeply enough into either mathematics or physics, the lines between them start to blur in an odd way. At the same time that Somerville was becoming fascinated by these connections, the scientific community was heading in the opposite social direction, with men developing into specialists.
As a result, an unexpectedly mundane issue arose. Please tell me what you intended to refer to each of these individuals as. Physicists, geologists, astronomers, mathematicians, and chemists were all present. “Man of science” and “natural philosopher” were both in use in the past, but neither of these terms adequately described what was happening in English at the time. When discussing this naming issue in 1834 while reviewing Somerville’s On the Connexion of the Physical Sciences, William Whewell proposed a rather strange new term that was derived from the word “artist”: scientist. At first, I thought the word was hilarious. Obviously so. So long as everyone agrees to stop making fun of it, every word is slightly absurd.
Whewell basically looks at Mary Somerville, realizes that she knows far too many subjects to fit into any existing category, and has to invent the word “scientist” for her on the spot, according to a very alluring online version of this story. While I would love it if that were the whole story, it’s that good, history is usually less accommodating. Whewell was addressing a broader issue with scientific terminology; the term was not coined solely due to Somerville’s gender or the fact that her field of study was a matter of debate. The genuine serendipitous event is nearly more delightful. In a dispute sparked by a book written specifically to unite increasingly fragmented scientific disciplines, the word “scientist” was first used in print. A new term for scientists was introduced, and a book was also released, urging us to view science as an integrated whole.
Somerville also has a surprisingly contemporary vibe to it. This is the age of ridiculous specialization. Most of the time, this will happen. But there is now an overwhelming amount of information that no one individual can adequately process. No matter how much time you devote to studying a specific area of theoretical physics, there will always be some literature that your neighbors assume you are familiar with. If you haven’t figured out how to make the day last much longer than twenty-four hours, the ideal of being mathematically and physically literate is meaningless. It seems, though, that the intriguing concepts did not get the message. But they persist in traversing borders. Geometry deviated into weightlessness. When it entered particle physics, group theory made a mistake. A collision with condensed matter occurred in topology. Black holes were discovered by information theory. Occasionally, two fields develop independently for many decades until someone realizes that they have been discussing remarkably similar topics in different languages.
I suppose that explains my fondness for Mary Somerville. She “refused to pick one thing,” but that’s just for the Instagram slide. This is because, thankfully, she had developed the habit of asking for connections between concepts, which proved to be quite useful. Something is different. Possessing a broad range of interests can make you seem diverse, but it takes actively seeking out connections to make that diversity matter. If you want to be a good physicist, you need to know how to ask certain kinds of questions, and math is a necessary part of that. When you read about astronomy, a mechanical problem takes on a new perspective. You learn the language of a thought that was hazy in one field to another. As time goes on, it becomes increasingly difficult to distinguish between one interest and the next.
Even in her latter years, ninety-one-year-old Mary Somerville continued to pursue mathematics. She passed away in 1872. A lot bigger scientific world was emerging than the one she had just entered, and the little word that Whewell had proposed was gradually becoming the norm. The term “scientist” is so commonplace now that a six-year-old can confidently declare their goal of becoming one without realizing that the same statement would have had no meaning two centuries ago. It has an eerie beauty to it. Mary Somerville’s life happened at the exact moment when our conception of science was shifting, not because she was the sole inventor of the word (though she certainly wasn’t). Amidst all this discord, someone was penning a book about the interconnectedness of all things, and new terms were being coined.
I don’t want to put life through that, but maybe there’s a lesson there. However, if you want to become an expert in a certain field, specialization is inevitable. That is an undeniable fact. However, you shouldn’t allow the label of your field to dictate the limits of your interest. Sitting next to you in the other class can be just what you need sometimes. Each time, it’s five topics up there. Starting with a question can lead you astray from your original intent; it is only after you’ve gone on that you come to appreciate the journey as its most captivating aspect.
I used to think complex numbers were the end of the line. You start with counting, patch in zero and the negatives, then the fractions, then the irrationals to fill the gaps between them, and finally you throw in $i$ to answer $x^2 = -1$. At that point the number system feels closed. Complete. Nothing left to want. I believed that for a long time.
Then I was talking with a friend about where complex numbers actually show up in the physical world, and he asked, almost in passing, whether I had ever looked into quaternions. I had heard the word. I could have told you they were “like complex numbers but bigger” and left it there, which is another way of admitting I did not understand them at all. He does. He has a PhD in computer science and the rare gift of making a hard idea feel inevitable instead of arbitrary, and over one long conversation, followed by most of a weekend spent chasing down the details he left me with, quaternions turned from a word I recognized into a structure I could see. This post is my attempt to hand you that same shift.
The place to start is not with quaternions but with a small miracle hiding inside the complex numbers you already know, because quaternions exist only to reproduce that miracle one dimension up. Take a point in the plane, written as $z = x + iy$, and multiply it by a unit complex number $w = \cos\phi + i\sin\phi$. Multiply the two out, keeping $i^2 = -1$:
Look at the two real numbers that came out. The new $x$-coordinate is $x\cos\phi – y\sin\phi$, the new $y$-coordinate is $x\sin\phi + y\cos\phi$, and that is exactly the rotation you met in the vectors post, the same $R_{ij}$ dictionary that turned one observer’s components into another’s. Multiplying a point by a unit complex number rotates that point about the origin by the angle $\phi$. Rotation, one of the most geometric operations there is, turns out to be nothing more than multiplication in disguise.
Figure 1: Multiplying the complex number z by the unit complex number w = cos φ + i sin φ rotates it about the origin through the angle φ while preserving its distance from the origin. The dashed arc shows the circular path of the endpoint.
Two things make this work, and both matter for what comes next. The first is that the length of $z$ is unchanged: a unit $w$ rotates without stretching. That is a statement about a norm, $|wz| = |w|\,|z|$, and squaring it gives the identity $(a^2+b^2)(c^2+d^2) = (ac-bd)^2 + (ad+bc)^2$, which says that the product of two sums of two squares is again a sum of two squares.1This is the Brahmagupta–Fibonacci identity, known in essence to Diophantus. It is exactly the statement that the complex norm is multiplicative, $|w z|^2 = |w|^2 |z|^2$, written out in real coordinates with $w = a+bi$ and $z = c+di$. Keep an eye on the number of squares: two here, and the question of whether the same trick works for three and for four squares is the whole story below. The second is that rotations compose the way multiplication does, and this is worth doing rather than asserting. Rotate by $\phi_1$ and then by $\phi_2$, which means multiply by the two unit numbers in turn:
The angle-addition formulas collapse the right-hand side to $\cos(\phi_1+\phi_2) + i\sin(\phi_1+\phi_2)$, a single rotation through the summed angle. Composing two rotations is nothing but multiplying the two numbers that represent them. Hold onto that sentence. It is the one demand quaternions will have to meet in three dimensions, where composing rotations is a famously trickier business than adding angles. In the plane, rotation is easy and complex numbers are the reason.
Now the obvious hunger sets in. We live in three dimensions, not two, and rotations in space are far richer and far more troublesome than rotations in a plane. If a two-dimensional number system hands us plane rotations for free, surely a three-dimensional one would hand us space rotations the same way. This is precisely the thought that seized William Rowan Hamilton in the 1830s, and he chased it for the better part of a decade.2Hamilton was searching for what he called “triples,” three-dimensional numbers $a + bi + cj$ that could be added and, crucially, multiplied while preserving length, so that unit triples would rotate space the way unit complex numbers rotate the plane. He worked at it for years; there is a much-repeated story that his children would ask him at breakfast whether he could multiply triples yet, and he would have to answer that he could still only add and subtract them.
Set up his problem the natural way. A triple is $a + bi + cj$, with two imaginary units obeying $i^2 = j^2 = -1$, and we want a multiplication that keeps the norm $a^2 + b^2 + c^2$ multiplicative, so that unit triples act as rotations. The entire difficulty collapses onto a single question: what is the product $ij$? For the triples to be closed under multiplication, $ij$ has to be another triple, some combination
$$ ij = \alpha + \beta i + \gamma j $$
with real coefficients $\alpha, \beta, \gamma$. Watch what that assumption forces. Multiply both sides on the left by $i$, and use $i^2 = -1$ together with ordinary associativity:
$$ i(ij) = i^2 j = -j, $$
while the right-hand side becomes
$$ i(\alpha + \beta i + \gamma j) = \alpha i + \beta i^2 + \gamma\, ij = -\beta + \gamma\alpha + (\alpha + \gamma\beta)\,i + \gamma^2 j, $$
where in the last step I substituted the assumed form of $ij$ back into itself. These two expressions must be equal, so match the coefficient of $j$ on each side:
$$ -1 = \gamma^2. $$
There is no real number $\gamma$ whose square is $-1$. The assumption that $ij$ is a triple has destroyed itself. No matter how you try to define the product, $ij$ simply refuses to live in the three-dimensional space spanned by $1$, $i$, and $j$; the multiplication insists on pointing somewhere outside. Hamilton spent years trying to force it back inside, and it could not be done. There is even a clean reason from number theory that it can never be done for any definition at all: the product of two sums of three squares is not always a sum of three squares, so no multiplication on triples can keep the norm multiplicative.3Take $3 = 1^2+1^2+1^2$ and $5 = 0^2+1^2+2^2$, both sums of three squares. Their product is $15$, and $15$ is not a sum of three squares: by Legendre’s three-square theorem a positive integer fails to be a sum of three squares exactly when it has the form $4^a(8b+7)$, and $15 = 8(1)+7$. So the three-square analogue of the identity in the previous footnote is simply false, and a length-preserving multiplication on triples cannot exist. The obstruction is not a failure of imagination; it is a theorem.
The way out is the move that looks like a defeat and turns out to be the discovery. If $ij$ will not stay inside three dimensions, then stop fighting it and give it a home. Call it a new unit, $k := ij$, orthogonal to $1$, $i$, and $j$, and accept that the numbers we are building are four-dimensional, not three. This is what Hamilton finally saw on 16 October 1843, walking with his wife along the Royal Canal in Dublin, and he was so certain it was right and so afraid he would forget it that he carved the defining relation into the stone of Broom Bridge with a knife:
$$ \boxed{\; i^2 = j^2 = k^2 = ijk = -1. \;} $$
A quaternion is a number $q = a + bi + cj + dk$, one real coordinate and three imaginary ones, and everything about how they multiply is squeezed out of that one carved formula. The remarkable thing is how little you have to be told. Those four equalities are the entire definition; every product of $i$, $j$, and $k$ follows from them by algebra alone, and it is worth turning the crank once to see that it really does.
Start from $ijk = -1$ and multiply on the right by $k$. Since $k^2 = -1$,
$$ ijk\cdot k = -k \;\Longrightarrow\; ij\,(k^2) = -k \;\Longrightarrow\; ij\,(-1) = -k \;\Longrightarrow\; ij = k. $$
Now multiply the same $ijk = -1$ on the left by $i$, using $i^2 = -1$:
$$ i\cdot ijk = -i \;\Longrightarrow\; (i^2) jk = -i \;\Longrightarrow\; -jk = -i \;\Longrightarrow\; jk = i. $$
The third product comes out by the same style of manipulation. Multiply $ij = k$ on the left by $k$, so that $k(ij) = k^2 = -1$; reading the left side as $(ki)j$ by associativity gives $(ki)j = -1$, and right-multiplying by $j$ with $j^2 = -1$ turns this into $(ki)(j^2) = -j$, that is $-(ki) = -j$, so
$$ ij = k, \qquad jk = i, \qquad ki = j. $$
The three units cycle, each one times its neighbor returning the third, in the order $i \to j \to k \to i$. Now the order of multiplication, which is where quaternions earn their keep. These products do not commute, and that is the point of the construction rather than a blemish on it. Left-multiply $jk = i$ by $j$:
$$ j(jk) = ji \;\Longrightarrow\; (j^2)k = ji \;\Longrightarrow\; ji = -k, $$
the exact negative of $ij$. The same one-line move on the other two, left-multiplying $ij = k$ by $i$ and right-multiplying $ij = k$ by $j$, gives $ik = -j$ and $kj = -i$. Reversing any product flips its sign, so the forward cycle $i \to j \to k \to i$ carries a plus and the backward cycle a minus,
$$ ij = k,\quad jk = i,\quad ki = j, \qquad ji = -k,\quad kj = -i,\quad ik = -j. $$
Figure 2: The multiplication rule in one picture. Following the arrows, $ij = k$, $jk = i$, $ki = j$. Going against them flips the sign: $ji = -k$, $kj = -i$, $ik = -j$. Order matters, which is exactly the feature that will let quaternions carry rotations that themselves depend on order.
That last sentence deserves weight. In the vectors post we watched two spatial rotations fail to commute: rotate a book about one axis then another, reverse the order, and the book lands somewhere else, while ordinary addition never cares about order. A number system meant to encode rotations therefore cannot be commutative, and the quaternions are not. The noncommutativity that looks at first like a bug is the precise property that fits them to the job. They gave up an old comfort, $xy = yx$, and bought a new capability with it.
Before going further it pays to learn to multiply any two quaternions at once, not just the units, because a single formula will carry the entire rest of the post. Split a quaternion into its real and imaginary parts exactly as you split a complex number, and give the imaginary part the name it deserves: for $q = a + bi + cj + dk$ write $q = a + \mathbf{u}$, where $a$ is the scalar part and $\mathbf{u} = bi + cj + dk$ is a pure quaternion, an honest three-dimensional vector wearing $i, j, k$ as its axes. Multiplying two pure quaternions is now just the multiplication table applied to every cross term. Take $\mathbf{u} = u_1 i + u_2 j + u_3 k$ and $\mathbf{v} = v_1 i + v_2 j + v_3 k$ and expand the nine products. The three like terms collapse by $i^2 = j^2 = k^2 = -1$ into a single real number, $-(u_1 v_1 + u_2 v_2 + u_3 v_3)$. The six unlike terms fold into the cycle: $u_1 v_2\,(ij) + u_2 v_1\,(ji) = (u_1 v_2 – u_2 v_1)\,k$, and its two rotations give the $i$ and $j$ parts. Collecting everything,
Both pieces are old acquaintances in disguise, and I will name them properly at the very end of this post, once they can land with full force. For now, call the scalar part $-\,\mathbf{u}\cdot\mathbf{v}$ and the vector part $\mathbf{u}\times\mathbf{v}$, so that the product of two pure quaternions is
One consequence is immediate and we will lean on it hard. Set $\mathbf{v} = \mathbf{u}$: the cross product of a vector with itself vanishes, so $\mathbf{u}^2 = -\,\mathbf{u}\cdot\mathbf{u} = -|\mathbf{u}|^2$, and a pure unit quaternion squares to $-1$, exactly as $i$ does. There is not a single square root of $-1$ living in the quaternions but a whole sphere of them, one for every direction in space. From here the product of two arbitrary quaternions follows by distributing over the split. With $p = a + \mathbf{u}$ and $q = b + \mathbf{v}$,
$$ pq = (a + \mathbf{u})(b + \mathbf{v}) = ab + a\mathbf{v} + b\mathbf{u} + \mathbf{u}\mathbf{v} = (ab – \mathbf{u}\cdot\mathbf{v}) + (a\mathbf{v} + b\mathbf{u} + \mathbf{u}\times\mathbf{v}). $$
The scalar part of a product is $ab – \mathbf{u}\cdot\mathbf{v}$ and the vector part is $a\mathbf{v} + b\mathbf{u} + \mathbf{u}\times\mathbf{v}$. Every calculation left in this post is a special case of that one line.
The norm we were chasing all along now falls out cleanly, and it is genuinely the multiplicative norm. Define the conjugate of $q = a + bi + cj + dk$ by flipping the sign of the imaginary part, $\bar q = a – bi – cj – dk$, in exact analogy with the complex conjugate. Multiply $q$ by $\bar q$ and every cross term cancels in pairs, because $ij$ and $ji$ carry opposite signs, leaving only the squares:
a real, nonnegative number, the squared length of $q$ as a four-dimensional vector. Multiplicativity now rests on a single fact about conjugation, and with the product formula in hand we can prove it outright instead of borrowing it: conjugating a product reverses the order of the factors, $\overline{pq} = \bar q\,\bar p$.4This order reversal is the same phenomenon as $(AB)^\top = B^\top A^\top$ for matrices, and it is not a coincidence: quaternions can be represented by $2\times 2$ complex matrices, with the conjugate becoming the conjugate transpose. It also rhymes with the covariant bookkeeping from the vectors post, where a transforming quantity reverses order relative to the basis. Order reversal under an involution is a recurring signature of noncommutative structures. Conjugation flips the sign of the vector part, so from the general product $pq = (ab – \mathbf{u}\cdot\mathbf{v}) + (a\mathbf{v} + b\mathbf{u} + \mathbf{u}\times\mathbf{v})$ we read off
Now multiply the conjugates in the opposite order. Using $\mathbf{v}\mathbf{u} = -\mathbf{v}\cdot\mathbf{u} + \mathbf{v}\times\mathbf{u} = -\mathbf{u}\cdot\mathbf{v} – \mathbf{u}\times\mathbf{v}$, since the dot product is symmetric and the cross product changes sign when its arguments swap,
$$ \bar q\,\bar p = (b – \mathbf{v})(a – \mathbf{u}) = ab – b\mathbf{u} – a\mathbf{v} + \mathbf{v}\mathbf{u} = (ab – \mathbf{u}\cdot\mathbf{v}) – (a\mathbf{v} + b\mathbf{u} + \mathbf{u}\times\mathbf{v}), $$
which is exactly $\overline{pq}$. The order reversal is not an analogy imported from matrices; it is a two-line consequence of how the vector part multiplies. Because $N(p) = p\bar p$ is a real scalar and commutes with everything, the multiplicativity is then three short steps:
$$ N(pq) = (pq)\overline{(pq)} = pq\,\bar q\,\bar p = p\,N(q)\,\bar p = N(q)\,p\bar p = N(p)\,N(q). $$
The norm is multiplicative. Written out in the sixteen real products this becomes Euler’s four-square identity, the four-dimensional cousin of the two-square identity we started with, and the reason four is the magic dimension where three failed.5Euler’s four-square identity, that a product of two sums of four squares is again a sum of four squares, predates Hamilton by nearly a century; he found it in 1748 while working on Lagrange’s four-square theorem. Hamilton’s quaternions are, in a sense, the hidden reason the identity is true: it is just $N(pq) = N(p)N(q)$ in coordinates. The pattern of which dimensions admit such an identity is severe. Sums-of-squares multiply cleanly only in dimensions $1$, $2$, $4$, and $8$, a result of Hurwitz, and the four cases are the real numbers, the complex numbers, the quaternions, and the octonions. There is no five-dimensional or seven-dimensional number system of this kind, and there never will be.
One consequence is worth pulling out on its own, because it is the property that makes quaternions a genuine number system and not merely an interesting array. Since $q\bar q = N(q)$ is a positive real number for any $q \ne 0$, we can divide by it, and
$$ q^{-1} = \frac{\bar q}{N(q)} $$
is a true multiplicative inverse: $q\,q^{-1} = q\bar q / N(q) = 1$. Every nonzero quaternion can be divided by. A system where you can add, subtract, multiply, and divide, with division always possible except by zero, is called a division algebra, and division algebras are extraordinarily rare. Over the real numbers, if you insist on the ordinary law $x(yz) = (xy)z$, there are exactly three: the real numbers, the complex numbers, and the quaternions, and nothing else, ever.6This is Frobenius’s theorem, 1877: the only finite-dimensional associative division algebras over the reals are $\mathbb{R}$, $\mathbb{C}$, and $\mathbb{H}$ (the $\mathbb{H}$ is for Hamilton). Each rung of the ladder costs you a property. Passing from $\mathbb{R}$ to $\mathbb{C}$ you lose ordering: there is no consistent way to say one complex number is less than another. Passing from $\mathbb{C}$ to $\mathbb{H}$ you lose commutativity, as we have seen. If you are willing to give up associativity as well you can reach one dimension further, to the eight-dimensional octonions $\mathbb{O}$, and there the ladder ends for good. This is the same “the language runs out” story that ran through the earlier posts, now with a precise list of exactly how far it can possibly go.
So we have our four-dimensional number system, and it does everything complex numbers did. The question that started Hamilton off is still open, though: does it actually rotate space? It does, and the way it does is stranger and more beautiful than the two-dimensional case, so we should build it carefully rather than assert it.
First, identify ordinary three-dimensional vectors with the quaternions that have no real part. A vector $\mathbf{v} = (v_1, v_2, v_3)$ becomes the pure quaternion $v = v_1 i + v_2 j + v_3 k$, and a quaternion is pure exactly when $\bar v = -v$. These pure quaternions are our three-dimensional space, sitting inside the four-dimensional quaternions as the imaginary part. The recipe for rotating them is the sandwich product: pick a unit quaternion $q$, meaning $N(q) = 1$, and send
$$ v \;\longmapsto\; q\,v\,q^{-1} = q\,v\,\bar q, $$
where the second equality holds because $q^{-1} = \bar q$ when $q$ is a unit quaternion. I claim this map is a rotation of three-dimensional space, and there are three things to check, each of which is a short calculation rather than an act of faith.
It preserves length, because the norm is multiplicative: $N(qv\bar q) = N(q)\,N(v)\,N(\bar q) = 1\cdot N(v)\cdot 1 = N(v)$, using $N(q) = N(\bar q) = 1$. It sends pure quaternions to pure quaternions, so a vector maps to a vector and never sprouts a real part, which you can see by conjugating the result and finding it flips sign:
$$ \overline{q v \bar q} = \overline{\bar q}\;\bar v\;\bar q = q\,(-v)\,\bar q = -\,q v \bar q, $$
where $\bar v = -v$ because $v$ is pure. A length-preserving linear map of three-dimensional space that fixes the origin is either a rotation or a reflection, and connectedness rules out the reflection.7Such a map has determinant $+1$ (a rotation) or $-1$ (a reflection), with nothing in between allowed. The map $v \mapsto q v \bar q$ depends continuously on $q$, and at $q = 1$ it is the identity, whose determinant is $+1$. The unit quaternions form the three-dimensional sphere $S^3$, which is connected, so as $q$ ranges over it the determinant cannot jump from $+1$ to $-1$ without passing through some forbidden value in between. It stays $+1$ throughout, and every sandwich by a unit quaternion is a rotation. The only question left is which rotation, and the quaternion announces its axis and angle directly. Write the unit quaternion in a form that should look suspiciously like Euler’s formula,
where $\hat{\mathbf{n}} = n_1 i + n_2 j + n_3 k$ is a pure unit quaternion built from a unit vector, so that $\hat{\mathbf{n}}^2 = -1$ by the sphere-of-square-roots fact from earlier. Two short calculations settle exactly what the sandwich does, and neither asks you to take anything on faith.
First, the axis does not move. Because $\hat{\mathbf{n}}^2 = -1$, the quaternion $q$ commutes with $\hat{\mathbf{n}}$:
so $q\,\hat{\mathbf{n}}\,\bar q = \hat{\mathbf{n}}\,q\bar q = \hat{\mathbf{n}}$. A vector pointing along $\hat{\mathbf{n}}$ is left precisely where it was, which is what it means for $\hat{\mathbf{n}}$ to be the axis.
Second, and here the half-angle finally explains itself, take a vector $\mathbf{v}$ perpendicular to the axis and watch it turn. Perpendicularity means $\hat{\mathbf{n}}\cdot\mathbf{v} = 0$, so by the product formula $\hat{\mathbf{n}}\mathbf{v} = \hat{\mathbf{n}}\times\mathbf{v}$ while $\mathbf{v}\hat{\mathbf{n}} = \mathbf{v}\times\hat{\mathbf{n}} = -\hat{\mathbf{n}}\times\mathbf{v}$; the two anticommute, $\mathbf{v}\hat{\mathbf{n}} = -\hat{\mathbf{n}}\mathbf{v}$, and a short step gives $\hat{\mathbf{n}}\mathbf{v}\hat{\mathbf{n}} = -\hat{\mathbf{n}}^2\mathbf{v} = \mathbf{v}$. Abbreviate $c = \cos\tfrac{\theta}{2}$ and $s = \sin\tfrac{\theta}{2}$ and expand the sandwich term by term:
Feed in the two facts just found, $\mathbf{v}\hat{\mathbf{n}} = -\hat{\mathbf{n}}\mathbf{v}$ and $\hat{\mathbf{n}}\mathbf{v}\hat{\mathbf{n}} = \mathbf{v}$. The two middle terms now add instead of cancelling, and the last becomes $-s^2\mathbf{v}$:
where the double-angle identities $c^2 – s^2 = \cos\theta$ and $2cs = \sin\theta$ did the folding, and $\hat{\mathbf{n}}\mathbf{v} = \hat{\mathbf{n}}\times\mathbf{v}$ because the two are perpendicular. Read the answer. The vector $\mathbf{v}$ has become $\cos\theta\,\mathbf{v} + \sin\theta\,(\hat{\mathbf{n}}\times\mathbf{v})$, which is $\mathbf{v}$ swung through the angle $\theta$ inside the plane perpendicular to the axis, with $\hat{\mathbf{n}}\times\mathbf{v}$ supplying the perpendicular direction it swings toward. That is a rotation by $\theta$, no more and no less. The half-angle is now no mystery at all: the sandwich touches $\mathbf{v}$ with $q$ on the left and $\bar q$ on the right, so the angle stored in $q$ is spent twice, and the double-angle formula is precisely the arithmetic of spending it twice. A general vector splits into a piece along the axis, which we showed is fixed, and a piece perpendicular to it, which turns by $\theta$; together they carry the whole vector through $\theta$ about $\hat{\mathbf{n}}$.8The result $q\mathbf{v}\bar q = \cos\theta\,\mathbf{v} + \sin\theta\,(\hat{\mathbf{n}}\times\mathbf{v})$ for $\mathbf{v}\perp\hat{\mathbf{n}}$, extended to a general vector by adding the fixed axial part, is Rodrigues’ rotation formula, usually derived with far more trigonometry. The quaternion sandwich produces it in four lines. A unit quaternion storing half the angle for a two-sided action is the same bookkeeping as a square root storing half an exponent; both act, in effect, on both sides.
The general formula is easy to trust once you have watched it collapse to numbers. Let me rotate the $x$-axis by ninety degrees about the $z$-axis, which should carry it onto the $y$-axis. Take $v = i$, the $x$-direction, and $\hat{\mathbf{n}} = k$, the $z$-direction, with $\theta = 90^\circ$, so the quaternion uses half of that, forty-five degrees:
$$ q v \bar q = \tfrac{1}{2}(i+j)(1-k) = \tfrac{1}{2}\,(i – ik + j – jk) = \tfrac{1}{2}\,(i + j + j – i) = j. $$
The $x$-axis lands exactly on the $y$-axis. The forty-five degrees inside the quaternion produced a ninety-degree turn in space, and the arithmetic that got us there used nothing but the multiplication table carved on the bridge.
Figure 3: The sandwich product $q\,v\,\bar q$ turns the vector $v$ about the axis $\hat{\mathbf{n}}$ through the angle $\theta$, while the quaternion itself stores only half that angle. The dashed ellipse is the circle $v$ sweeps out as it rotates around the axis.
This construction also settles the demand we carried over from the complex numbers, that composing rotations should amount to multiplying the numbers that represent them. Rotate first by $q_1$ and then by $q_2$, and stack the sandwiches:
$$ v \;\longmapsto\; q_2\,(q_1 v \bar q_1)\,\bar q_2 = (q_2 q_1)\,v\,(\bar q_1 \bar q_2) = (q_2 q_1)\,v\,\overline{(q_2 q_1)}, $$
where the last step is the order-reversal $\bar q_1 \bar q_2 = \overline{q_2 q_1}$ we proved earlier. The composite rotation is a single sandwich by the product $q_2 q_1$. Composition of rotations is multiplication of quaternions, just as it was for complex numbers, with the one difference that is now the entire point: quaternion multiplication does not commute, so $q_2 q_1 \ne q_1 q_2$ in general, and neither do the rotations. The book that landed in two different places depending on the order of two turns, back in the vectors post, was showing you the noncommutativity of this product before we had the product to show. Associativity does survive, since quaternion multiplication is associative, so long chains of rotations compose without any ambiguity about grouping; only their order matters. That is exactly the structure that rotations of space actually have.
The half-angle has a consequence that reaches straight into physics. Because the quaternion holds $\theta/2$, turning all the way around by $\theta = 360^\circ$ makes the quaternion carry $180^\circ$, at which point $\cos 180^\circ + \sin 180^\circ\,\hat{\mathbf{n}} = -1$. The quaternion has become $-1$, not $+1$, even though space itself has returned exactly to where it started. You have to go around twice, a full $720^\circ$, before the quaternion comes back to $+1$. Two different quaternions, $q$ and $-q$, describe the very same rotation of space, so the map from quaternions to rotations is two-to-one.9This is the famous double cover, written $SU(2) \to SO(3)$ in the language of groups. It is not a defect of the description; it is a real feature of rotations that only becomes visible once you have the right language. The plate trick, or Dirac’s belt trick, makes it physical: hold a belt buckle fixed, rotate the free end by $360^\circ$, and the belt is twisted; rotate by another $360^\circ$, the same way, and the twist comes out. Space returns to itself after $360^\circ$, but the connection to its surroundings does not, and needs $720^\circ$. In quantum mechanics this is exactly the behavior of a spin-$\tfrac12$ particle, whose wavefunction changes sign under a single full rotation and only returns after two, and Hamilton’s $i$, $j$, $k$ are, up to a factor of $i$, the Pauli matrices that describe it. Space forgets the difference between a $360^\circ$ turn and no turn at all; the quaternions remember it. That memory is not an artifact. It is the same double-valuedness that a spinning electron carries in quantum mechanics, and quaternions saw it a century before anyone knew to look for it.
Now the promise I made in the middle of the post, to name the two pieces of the pure-quaternion product properly once they could land with full force. Recall what we found the moment we first multiplied two pure quaternions, two ordinary vectors dressed in $i, j, k$:
The scalar part is minus the dot product, the very inner product the vectors post built an entire discussion around. The vector part is the cross product, component for component. Both of the products that structure all of three-dimensional vector algebra are already inside quaternion multiplication, the scalar shadow and the vector shadow of one operation, and we have been using them without comment for the whole second half of this post: they were in the norm, in the order-reversal proof, in the rotation formula. This is not an analogy. It is the historical source. The dot product and the cross product were carved out of exactly this quaternion product by Gibbs and Heaviside in the 1880s, who kept the two halves and discarded the quaternion framework around them, which is why students meet the dot and cross products today as two separate definitions dropped from the sky rather than as the scalar and vector shadows of a single multiplication.10I mentioned in the vectors post that Gibbs and Heaviside carved modern vector analysis out of Hamilton’s quaternions; this identity is precisely what they carved. There was a genuine and rather heated dispute in the 1890s between the quaternion loyalists, led by Peter Tait, and the vector-analysis reformers, over which language physics should speak. The reformers won for practical purposes, and the dot and cross products became standard, but the quaternions never really left. They returned through the back door in quantum mechanics as the algebra of spin, and through the front door in computer graphics, robotics, and spacecraft attitude control, where a rotation stored as a single unit quaternion avoids the gimbal lock and interpolation headaches that plague rotation matrices and Euler angles. Your friend with the computer science PhD almost certainly met them there.
That is where my weekend of reviewing left me, and where I will leave you. Complex numbers were not the end of the number line rather only a rung on it. One dimension up, at the price of commutativity, the quaternions bought us the ability to turn objects in space, a second copy of the rotation miracle that complex numbers perform in the plane, and they folded the dot and cross products into a single act of multiplication as a kind of parting gift. The ladder climbs one rung further, to the octonions, at the price of associativity, and then, by a theorem, it stops. Not because no one has been clever enough to continue it, but because the mathematics has been shown to have no further rung to offer. Some staircases really do have a top, and it is a strange comfort to know exactly where it is.
References and Footnotes
1
This is the Brahmagupta–Fibonacci identity, known in essence to Diophantus. It is exactly the statement that the complex norm is multiplicative, $|w z|^2 = |w|^2 |z|^2$, written out in real coordinates with $w = a+bi$ and $z = c+di$. Keep an eye on the number of squares: two here, and the question of whether the same trick works for three and for four squares is the whole story below. ↩︎
2
Hamilton was searching for what he called “triples,” three-dimensional numbers $a + bi + cj$ that could be added and, crucially, multiplied while preserving length, so that unit triples would rotate space the way unit complex numbers rotate the plane. He worked at it for years; there is a much-repeated story that his children would ask him at breakfast whether he could multiply triples yet, and he would have to answer that he could still only add and subtract them. ↩︎
3
Take $3 = 1^2+1^2+1^2$ and $5 = 0^2+1^2+2^2$, both sums of three squares. Their product is $15$, and $15$ is not a sum of three squares: by Legendre’s three-square theorem a positive integer fails to be a sum of three squares exactly when it has the form $4^a(8b+7)$, and $15 = 8(1)+7$. So the three-square analogue of the identity in the previous footnote is simply false, and a length-preserving multiplication on triples cannot exist. The obstruction is not a failure of imagination; it is a theorem. ↩︎
4
This order reversal is the same phenomenon as $(AB)^\top = B^\top A^\top$ for matrices, and it is not a coincidence: quaternions can be represented by $2\times 2$ complex matrices, with the conjugate becoming the conjugate transpose. It also rhymes with the covariant bookkeeping from the vectors post, where a transforming quantity reverses order relative to the basis. Order reversal under an involution is a recurring signature of noncommutative structures. ↩︎
5
Euler’s four-square identity, that a product of two sums of four squares is again a sum of four squares, predates Hamilton by nearly a century; he found it in 1748 while working on Lagrange’s four-square theorem. Hamilton’s quaternions are, in a sense, the hidden reason the identity is true: it is just $N(pq) = N(p)N(q)$ in coordinates. The pattern of which dimensions admit such an identity is severe. Sums-of-squares multiply cleanly only in dimensions $1$, $2$, $4$, and $8$, a result of Hurwitz, and the four cases are the real numbers, the complex numbers, the quaternions, and the octonions. There is no five-dimensional or seven-dimensional number system of this kind, and there never will be. ↩︎
6
This is Frobenius’s theorem, 1877: the only finite-dimensional associative division algebras over the reals are $\mathbb{R}$, $\mathbb{C}$, and $\mathbb{H}$ (the $\mathbb{H}$ is for Hamilton). Each rung of the ladder costs you a property. Passing from $\mathbb{R}$ to $\mathbb{C}$ you lose ordering: there is no consistent way to say one complex number is less than another. Passing from $\mathbb{C}$ to $\mathbb{H}$ you lose commutativity, as we have seen. If you are willing to give up associativity as well you can reach one dimension further, to the eight-dimensional octonions $\mathbb{O}$, and there the ladder ends for good. This is the same “the language runs out” story that ran through the earlier posts, now with a precise list of exactly how far it can possibly go. ↩︎
7
Such a map has determinant $+1$ (a rotation) or $-1$ (a reflection), with nothing in between allowed. The map $v \mapsto q v \bar q$ depends continuously on $q$, and at $q = 1$ it is the identity, whose determinant is $+1$. The unit quaternions form the three-dimensional sphere $S^3$, which is connected, so as $q$ ranges over it the determinant cannot jump from $+1$ to $-1$ without passing through some forbidden value in between. It stays $+1$ throughout, and every sandwich by a unit quaternion is a rotation. ↩︎
8
The result $q\mathbf{v}\bar q = \cos\theta\,\mathbf{v} + \sin\theta\,(\hat{\mathbf{n}}\times\mathbf{v})$ for $\mathbf{v}\perp\hat{\mathbf{n}}$, extended to a general vector by adding the fixed axial part, is Rodrigues’ rotation formula, usually derived with far more trigonometry. The quaternion sandwich produces it in four lines. A unit quaternion storing half the angle for a two-sided action is the same bookkeeping as a square root storing half an exponent; both act, in effect, on both sides. ↩︎
9
This is the famous double cover, written $SU(2) \to SO(3)$ in the language of groups. It is not a defect of the description; it is a real feature of rotations that only becomes visible once you have the right language. The plate trick, or Dirac’s belt trick, makes it physical: hold a belt buckle fixed, rotate the free end by $360^\circ$, and the belt is twisted; rotate by another $360^\circ$, the same way, and the twist comes out. Space returns to itself after $360^\circ$, but the connection to its surroundings does not, and needs $720^\circ$. In quantum mechanics this is exactly the behavior of a spin-$\tfrac12$ particle, whose wavefunction changes sign under a single full rotation and only returns after two, and Hamilton’s $i$, $j$, $k$ are, up to a factor of $i$, the Pauli matrices that describe it. ↩︎
10
I mentioned in the vectors post that Gibbs and Heaviside carved modern vector analysis out of Hamilton’s quaternions; this identity is precisely what they carved. There was a genuine and rather heated dispute in the 1890s between the quaternion loyalists, led by Peter Tait, and the vector-analysis reformers, over which language physics should speak. The reformers won for practical purposes, and the dot and cross products became standard, but the quaternions never really left. They returned through the back door in quantum mechanics as the algebra of spin, and through the front door in computer graphics, robotics, and spacecraft attitude control, where a rotation stored as a single unit quaternion avoids the gimbal lock and interpolation headaches that plague rotation matrices and Euler angles. Your friend with the computer science PhD almost certainly met them there. ↩︎