Derivation of the Cosmic-Ray Diffusion Coefficient from Pitch-Angle Scattering

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

$$
\frac{d\mathbf{p}}{dt}=\frac{q}{c}\mathbf{v}\times\mathbf{B},
$$

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,

$$\frac{d\mu}{dt}=\frac{1}{v}\frac{dv_z}{dt}=\frac{q}{\gamma mcv}(\mathbf{v}\times\mathbf{B})_z=\frac{q}{\gamma mcv}(v_xB_y-v_yB_x).$$

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.

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.²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.

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.

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.³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. 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.

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.⁴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.

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$.⁵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.

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. ↩︎