23 Cosmic Ray
Cosmic rays (CRs), that is the population of charged, relativistic particles with non-thermal spectra, are ubiquitous in the Universe. They pervade systems of all sizes, from stellar systems to whole galaxies, from galaxy clusters to the intercluster medium. Modelling of cosmic ray transport and interpretation of cosmic ray data ultimately rely on a solid understanding of the interactions of charged particles with turbulent magnetic fields.
Mertsch (2020) presents a thorough guide of performing test particle simulation for cosmic rays.
23.1 Quasilinear Theory
More than 50 years since its inception, quasi-linear theory (QLT) [Jokipii (1966); Kennel and Engelmann 1966; Hall and Sturrock 1967; Hasselmann and Wibberenz 1970] is still very much the paradigm for phenomenological applications to Galactic cosmic rays. In QLT, the Fokker-Planck equation for the temporal evolution of the phase space density of CRs is derived in a perturbative approach where the force on a particle due to a turbulent magnetic field is evaluated along the unperturbed trajectory in a regular background field. The Fokker-Planck coefficients, most prominently the components of the spatial diffusion tensor, can be computed for a given model of turbulence, parametrised by the two-point function of the turbulent magnetic field. Famously, in QLT the interactions between plasma waves and particles are found to be resonant, meaning that particles of a certain gyroradius \(r_L = v/\Omega\), \(\Omega\) denoting the gyrofrequency, are only affected by waves with a wavenumber k that satisfies \(k r_L \mu \approx 1\) (for low-frequency waves like Alfvén waves), where \(\mu\) is the cosine of the pitch-angle, that is the angle between the particle momentum \(\mathbf{p}\) and the regular magnetic field \(\left<\mathbf{B}\right>\), \(\mu \equiv \mathbf{p}\cdot\left< \mathbf{B} \right> / (|\mathbf{p}||\left< \mathbf{B} \right>|)\).
While some of QLT’s predictions are qualitatively confirmed by data, e.g. the rigidity-dependence of the diffusion coefficients, there are a number of concerns. The most famous one is the \(90^\circ\) problem: Due to the resonance condition, particles with pitch-angle close to \(90^\circ\) (\(\mu\approx 0\)) can only be in resonance with very large wavenumbers k which for the usual turbulent spectra contain little energy. In the limit \(\mu\rightarrow 0\), the scattering rate vanishes and particles cannot change direction along the background field resulting in ballistic transport. This is obviously at variance with the diffusive transport inferred from observations.
The root cause of the \(90^\circ\) problem is the assumption of unperturbed trajectories in QLT. This is remedied in non-linear theories, where the decay of correlations leads to a broadening of the resonance condition which allows for efficient enough scattering through \(90^\circ\). However, any such extension of QLT requires additional assumptions, for instance on the form of temporal decorrelation of the particle’s trajectory.
The success and popularity of QLT can be ascribed to its conceptual simplicity and validity in a number of important environments, including the solar wind, the interstellar medium and galaxy clusters. In addition, QLT is simple in principle and thus allows for a straight-forward computation of the transport parameters, albeit it can become arbitrarily complex in practice. Finally, these results can be found to agree with inferences from observations, e.g. the normalisation and power law shape of the Galactic diffusion coefficient.
23.1.1 Derivation of the Fokker-Planck equation
For self-consistency, I will repeat all the equations here instead of linking to other parts of the note. Philipp Mertsch uses CGS unit in the following equations.
Charged particles in electric and magnetic fields \(\mathbf{E}\) and \(\mathbf{B}\) are subject to the Lorentz force, \[ \mathbf{F}_L = e\left( \mathbf{E} + \mathbf{v}\times\mathbf{B}/c \right) \tag{23.1}\]
It is customary to decompose the magnetic field into a large-scale, homogeneous, regular background field, \(\langle \mathbf{B} \rangle\) and a small-scale, turbulent, random field \(\delta\mathbf{B}\), that is \(\mathbf{B} = \langle\mathbf{B}\rangle\) with \(\langle\delta\mathbf{B}\rangle = 0\). The angled brackets denote averages over an ensemble of turbulent magnetic fields.
Without loss of generality, we assume in the following that the regular field is oriented along the z-direction, \(\langle\mathbf{B}\rangle = B_z\hat{z}\), unless stated otherwise. Large-scale electric fields are usually ignored, \(\langle\mathbf{E}\rangle = 0\), as the large mobility of charges in astrophysical plasmas is efficiently shielding against regular electric fields (that is on scales much larger than the Debye length). Small-scale electric fields \(\delta\mathbf{E}\) are necessarily present, but from Faraday’s induction law, their magnitude can be estimated to be \(|\delta\mathbf{E}| \sim ( v_{\text{A}} / c) |\delta\mathbf{B}|\) with \(v_\text{A}\) the Alfvén velocity and \(v_A/c \ll 1\) in most astrophysical environments.1 Thus, to lowest order in \((v_\text{A}/c)\), there is no electric field and as the magnetic force is not performing any work on the particle, particle energy is consequently conserved. Note that at higher orders in \((v_\text{A}/c)\), the particle energy can change in a second-order Fermi type process. For simplicity, we constrain ourselves here to considering the lowest order case which results in pitch-angle scattering.
A charged particle in a magnetic field forms a Hamiltonian system as long as dissipative processes (or any form of energy losses) can be ignored. A consequence of this is Liouville’s theorem, that is the conservation of phase space volume under canonical transformations. As time evolution is a canonical transformation, phases space volume is conserved in time. Together with particle number conservation this implies the conservation of phase space density \(f = f(\mathbf{r},\mathbf{p},t)\). This is conveniently captured by what we will call Liouville’s equation, \[ \frac{\mathrm{d}f}{\mathrm{d}t} = \frac{\partial f}{\partial t} + \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \cdot \mathbf{\nabla}_{\mathbf{r}} f + \frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t} \cdot \mathbf{\nabla}_{\mathbf{p}} f = 0 \tag{23.2}\] encoding the incompressibility of the phase space flow. Here, \[ \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} = \mathbf{v},\quad \frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t} = e\left( \mathbf{E} + \frac{\mathbf{v}\times\mathbf{B}}{c} \right) \tag{23.3}\] are the equations of motion. Note that a necessary condition for a Hamiltonian system is that the forces are conservative and differentiable (“p-divergence-free”, see more in (Bellan 2008)).
A collisionless plasma under the influence of external \(\mathbf{E}\) and \(\mathbf{B}\) fields is an example of a Hamiltonian system. Its Hamiltonian is (Jackson 1998) \[ H = \sqrt{(c\mathbf{P} - e\mathbf{A})^2 + m^2c^4} + e\phi \tag{23.4}\] where \(\mathbf{P} = \mathbf{p} + (e/c)\mathbf{A}\) is the canonical momentum, \(\mathbf{A}\) the potential vector, m the particle mass, e its charge and \(\phi\) the electric potential. Therefore, the phase space density of this collisionless plasma satisfies Equation 23.2 and substituting the Lorentz force, Equation 23.1, in Equation 23.2 gives the Vlasov equation, \[ \frac{\partial f}{\partial t} + \frac{\mathrm{d} \mathbf{r}}{\mathrm{d} t} \cdot \mathbf{\nabla}_{\mathbf{r}} f + e \left ( \mathbf{E} + \frac{\mathbf{v} \times \mathbf{B}}{c} \right ) \cdot \mathbf{\nabla}_{\mathbf{p}} f = 0 \tag{23.5}\] which together with Maxwell’s equations forms the basis of plasma kinetic theory. For a collisional plasma, a term needs to be added to the right-hand side, the famous collision operator. For a collisionless plasma (as appropriate for CRs) the right-hand side remains zero.
Considering turbulent fields, the phase space density also becomes a random field, \(f = \langle f\rangle + \delta f\), with an expectation value, \(\langle f \rangle\), and fluctuations around it, \(\delta f\), that satisfy \(\langle \delta f \rangle=0\).
In any realistic astrophysical situation, it is of course impossible to know the small-scale turbulent field at all positions in order to exactly solve Equation 23.5. Instead, one can only hope to predict statistical moments of the phase space density for a statistical ensemble of turbulent magnetic fields. Traditionally, one is mostly interested in the first moment, the ensemble average, while higher order moments is also possible to obtain.
For the reason explained above, we will ignore \(\mathbf{E}\) in the following. Averaging Equation 23.5 we find (Jokipii 1972) \[ \begin{aligned} \frac{\mathrm{d} \langle f \rangle}{\mathrm{d} t} &= \frac{\partial \langle f \rangle }{\partial t} + \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \cdot \mathbf{\nabla}_{\mathbf{r}} \langle f \rangle + e \frac{\mathbf{v} \times \langle \mathbf {B} \rangle}{c} \cdot \mathbf{\nabla}_{\mathbf{p}} \langle f \rangle \\ &= - \left \langle e \frac{\mathbf{v} \times \mathbf{\delta B}}{c} \cdot \nabla _{\mathbf{p}} \delta f \right \rangle \neq 0 \end{aligned} \tag{23.6}\]
Note that unlike the phase space density \(f\), the ensemble averaged phase space density \(\langle f \rangle\) is not conserved, \(\mathrm{d}\langle f\rangle /\mathrm{d}t \neq 0\). (More on this later!)
One way to glean some physical insight from Equation 23.6 is to identify its right-hand-side with a damping term, (Earl+ 1988; Webb 1989), \[ \left \langle e \frac{\mathbf{v} \times \mathbf{\delta B}}{c} \cdot \nabla_{\mathbf{p}} \delta f \right \rangle \to \nu \left( \langle f \rangle - \frac{1}{4\pi}\int\mathrm{d}\hat{\mathbf{p}} \, \langle f \rangle \right) \] (where \(\hat{\mathbf{p}}=\mathbf{p}/|\mathbf{p}|\)), that is driving the phase space density towards isotropy at a rate \(\nu\), an approach that can also be motivated by gas kinetic theory (Bhatnagar+ 1954). This way, Equation 23.6 can be solved and shown to lead to a spatial diffusion equation. The parallel diffusion coefficient can be identified as \(\kappa_\parallel = v^2 / (3\nu)\) whereas the perpendicular diffusion coefficient satisfies \(\kappa_\perp / \kappa_\parallel = (1 + \Omega^2 / \nu^2)^{-1}\), a result referred to as the classical scattering limit (Gleeson 1969). Here, \(\Omega\) is the particle’s gyrofrequency.
In QLT, however, a more systematic solution for \(f\) is sought through an equation for the temporal evolution of the fluctuations \(\delta f\). Such an equation can be obtained by subtracting the ensemble-averaged Vlasov Equation 23.6 from the original Vlasov Equation 23.5, \[ \begin{aligned} & \frac{\partial \delta f}{\partial t} + \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \cdot \mathbf{\nabla}_{\mathbf{r}} \delta f + e \left(\frac{\mathbf{v} \times \langle \mathbf{B}\rangle}{c}\right) \cdot \mathbf{\nabla}_{\mathbf{p}} \delta f \\ & \simeq -e \left(\frac{\mathbf{v} \times \mathbf{\delta B}}{c} \right) \cdot \mathbf{\nabla}_{\mathbf{p}} \langle f \rangle \end{aligned} \tag{23.7}\] Here, we have chosen to ignore the difference \[ e \frac{\mathbf{v} \times \mathbf{\delta B}}{c} \cdot \nabla_{\mathbf{p}} \delta f - \left \langle e \frac{\mathbf{v} \times \mathbf{\delta B}}{c} \cdot \nabla_{\mathbf{p}} \delta f \right \rangle \] which is second order in perturbed quantities, \(\delta\mathbf{B}\) and \(\delta f\). This assumes, of course, that \(|\delta\mathbf{B}|\ll |\langle \mathbf{B} \rangle|\) and therefore \(\delta f \ll \langle f \rangle\). Equation 23.7 can now be integrated with the method of characteristics, the formal solution being \[ \delta f = \delta f_{0} - \int _{t_{0}}^{t} \mathrm {d} t' \left[ e \left( \frac{\mathbf{v} \times \mathbf{\delta B}}{c} \right) \cdot \mathbf{\nabla}_{\mathbf{p}} \langle f \rangle \right]_{P(t')} \tag{23.8}\] Here, \(\delta f_{0} \equiv \delta f(\mathbf{r}, \mathbf{p}, t_{0})\) denotes the phase space density at time \(t_0\) and the subscript \(P(t')\) indicates that positions and momenta in the square brackets are to be evaluated along the characteristics of Equation 23.7, that is the solutions of the equations of motions, Equation 23.3 with \(\mathbf{B}\) replaced by the regular field \(\langle\mathbf{B}\rangle\) only (and again no electric field). These solutions \(P\) are commonly referred to as unperturbed orbits or unperturbed trajectories. For the homogeneous regular magnetic field \(\langle\mathbf{B}\rangle = B_z\hat{z}\) assumed here they are of course helices along the z-direction.
We can now substitue Equation 23.8 into Equation 23.6, \[ \begin{aligned} & \frac{\partial \langle f \rangle }{\partial t} + \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \cdot \mathbf{\nabla}_{\mathbf{r}} \langle f \rangle + e \frac{\mathbf{v} \times \langle \mathbf {B} \rangle}{c} \cdot \mathbf{\nabla}_{\mathbf{p}} \langle f \rangle \\ &\simeq \int _{t_{0}}^{t} \!\! \mathrm {d} t' \! \left \langle e \frac{\mathbf{v} \! \times \! \mathbf{\delta B}}{c} \! \cdot \! \nabla _{\mathbf{p}} \left [ e \frac{\mathbf{v} \! \times \! \mathbf{\delta B}}{c} \! \cdot \! \mathbf{\nabla}_{\mathbf{p}} \langle f \rangle \right ]_{P(t')} \right \rangle \end{aligned} \tag{23.9}\] where we have dropped the term \(\propto \delta f_0\). At this stage, we can already see that the right-hand side will lead to diffusion terms (courtesy of the two momentum derivatives) and that it depends on the turbulent magnetic field’s two-point function, integrated along the unperturbed trajectory \(P(t')\). To make further progress, we consider the momentum \(\mathbf{p}\) in spherical coordinates, that is \(\mathbf{p} = p\left( \sqrt{1-\mu^2}\cos\phi, \sqrt{1-\mu^2}\sin\phi, \mu \right)^T\) and introduce the correlation lengths \(l_c\) and correlation times \(\tau_c\) of the turbulent magnetic field, defined through2 \[ \langle \delta B^{2} \rangle l_{\text{c}} \equiv \int _{0}^{\infty} \mathrm {d} \Delta r \, \langle \delta B(\mathbf{r}) \delta B(\mathbf{r}+\mathbf{\Delta r}) \rangle \tag{23.10}\] \[ \langle \delta B^{2} \rangle \tau_{\text{c}} \equiv \int _{0}^{\infty} \mathrm {d} \Delta t \, \langle \delta B(t) \delta B(t+\Delta t) \rangle \tag{23.11}\] Here, \(\delta B(t)\) is short-hand for \(\left[\delta B(\mathbf{r}) \right]_{P(t)}\), that is \(\delta B(\mathbf{r})\) evaluated along the unperturbed trajectory at time t.
The right-hand side of Equation 23.9 is still rather unwieldy and further progress requires a number of assumptions. In addition to
- Smallness of perturbations, \(|\mathbf{\delta B}| \ll |\langle \mathbf{B} \rangle|\)
these are:
Gyrotropy: The ensemble-averaged phase space density \(\langle f \rangle\) does not depend on the azimuthal angle \(\phi\), so \(\langle f \rangle (\mathbf{r}, p, \mu, \phi, t) \to \langle f \rangle (\mathbf{r}, p, \mu, t)\).
Adiabatic approximation: The phase space density only varies on time-scales much larger than the correlation time of the turbulent magnetic field, \(\tau_c\), \[ \langle f \rangle \left / \frac{\partial \langle f \rangle }{\partial t} \right . \gg \tau _{\text{c}} \tag{23.12}\]
Finite correlation times: The correlation times of the turbulent magnetic field are much larger than the Larmor time, \(\tau_c \gg \Omega^{-1}\).
Homogeneous and stationary turbulence: the statistical moments are invariant under translations in space and time. In natural language, it means that the statistics of turbulence do not depend on when and where to take measurements.
Under these conditions, the ensemble averaged Vlasov equation ultimately results in a Fokker-Planck type equation (Fokker 1914; Planck 1917), also known as the Kolmogorov forward (Kolmogorov 1931) or as the Smoluchowski equation (Bogoliubov and Krylov 1939), describing diffusion in pitch-angle, \[ \frac{\partial \langle f \rangle }{\partial t} + v \mu \frac{\partial \langle f \rangle }{\partial z} = \frac{\partial }{\partial \mu } \left ( D_{\mu \mu } \frac{\partial \langle f \rangle }{\partial \mu } \right ) \, . \tag{23.13}\]
Following the approach sketched above, the pitch-angle diffusion coefficient, \[ D_{\mu \mu} \equiv \frac{\langle (\Delta \mu )^{2} \rangle }{2 \Delta t} \tag{23.14}\] can be expressed in terms of the correlation function of the magnetic field. (HOW???)
If we had not decided to ignore any electric field, additional terms would have appeared in the Fokker-Planck Equation 23.13, relating to changes in momentum \(\mathbf{p}\) and pitch-angle, with diffusion coefficients \(D_{\mu p} = D_{p\mu}\) and \(D_{pp}\) defined analogously to Equation 23.14.
We have furthermore assumed that \(v_A / v \ll 1\) in order for \(D_{xx}, D_{yy}, D_{xy}\) and \(D_{yx}\) to be negligible. Not doing so, would have resulted in the additional terms \[ \frac{\partial}{\partial x} \left (D_{xx} \frac{\partial \langle f \rangle }{\partial x} + D_{xy} \frac{\partial \langle f \rangle }{\partial y} \right ) + \frac{\partial}{\partial y} \left (D_{yx} \frac{\partial \langle f \rangle }{\partial x} + D_{yy} \frac{\partial \langle f \rangle }{\partial y} \right ) \] to be added to the right-hand side of Equation 23.13.
In summary, under the influence of a turbulent magnetic field, charged particles are performing a random walk in pitch-angle which in the ensemble average results in diffusion in pitch-angle (cosine).
23.1.2 The diffusion approximation
Particle transport can be conveniently categorised if the mean-square displacement in direction i, \(\langle \Delta x_i^2 \rangle\), has a power law dependence,3 \[ \langle \Delta {r}_{i}^{2} \rangle \propto (\Delta t)^{\alpha} \tag{23.15}\] as \[ \begin{array}{r@{\quad } l} \alpha < 1: & \text{sub-diffusive,} \\ \alpha = 1: & \text{diffusive,} \\ \alpha > 1:& \text{super-diffusive, in particular} \\ \alpha = 2: & \text{ballistic.} \end{array} \]
It seems clear that transport in any perturbative theory with \(|\mathbf{\delta B}| \ll |\langle \mathbf{B} \rangle|\) must be ballistic at early enough times: Particles just gyrate around \(\langle\mathbf{B}\rangle\) and \(\langle \Delta z^{2} \rangle = (v \mu \Delta t)^{2}\) while \(\langle \Delta x^{2} \rangle = \langle \Delta y^{2} \rangle = 0\) when integrated over full gyroperiods. At late times, that is for \(t \gg D_{\mu\mu}^{-1}\), we would expect diffusive behaviour for the transport along the field.
In order to formalise this picture, we derive a spatial diffusion equation from the Fokker-Planck Equation 23.13. To this end, we decompose \(\langle f\rangle\) into an isotropic part, \(g\), and an anisotropic part, \(h\), \[ \langle f \rangle (p, \mu, t) = g(p, t) + h(p, \mu, t) \] where \[ g(p, t) = \frac{1}{2}\int_{-1}^{1} \mathrm{d}\mu \langle f \rangle (p, \mu, t) \] and \[ \int_{-1}^{1} \mathrm{d}\mu\, h(p, \mu, t) = 0 \]
If \(g\) varies only slowly with time and position, \[ g \left / \frac{\partial g}{\partial t} \right . \gg \tau_{\text{sc}} \quad \text{and} \quad g \left / \frac{\partial g}{\partial z} \right . \gg \lambda_{\text{sc}} \] where \(\tau_{\text{sc}} \sim D_{\mu \mu }^{-1}\) and \(\lambda_{\text{sc}} \sim v \tau_{\text{sc}}\) are the scattering time and mean-free path, respectively, the phase space density will be very isotropic, \(h\ll g\). In this case, we can derive a spatial diffusion equation for the isotropic part \(g\) (e.g. Hasselmann and Wibberenz 1970), \[ \frac{\partial g}{\partial t} - \frac{\partial }{\partial z} \left( \kappa_{\parallel} \frac{\partial g}{\partial z} \right) = 0 \tag{23.16}\] with the parallel diffusion coefficient \[ \kappa_{\parallel} = \frac{v^{2}}{8} \int_{-1}^{1} \mathrm{d} \mu \frac{(1- \ \mu^{2})^{2}}{D_{\mu\mu}} \tag{23.17}\]
Furthermore, we would expect the anisotropic part \(h\) to be dominated by the dipole anisotropy (???), that is \(h\approx h_1\mu\) with \[ h_{1} = \frac{3}{2} \int_{-1}^{1} \mathrm{d} \mu \, \mu \, h(\mu) = - \frac{2}{v} \kappa_{\parallel} \frac{\partial g}{\partial z} \tag{23.18}\]
23.1.3 Computation of transport coefficients
So far, we have not specified the functional form of the Fokker-Planck coefficients, e.g. the pitch-angle diffusion coefficient \(D_{\mu\mu}\), and its dependence on the two-point correlation function of turbulence \(P_{ij}(\mathbf{k})\) that emerges in the derivation of the Fokker-Planck Equation 23.13. An alternative to the derivation of Section 23.1.1 is to directly compute the Fokker-Planck coefficients from solutions of the equations of motion. In fact, an arbitrary Fokker-Planck coefficient \(D_{PQ}\) can be defined in terms of the mean displacements of the variables in question, P and Q. For instance, the pitch-angle diffusion coefficient can be derived as the \(t\to\infty\) limit of the running diffusion coefficient, \[ d_{\mu\mu}(t) = \frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t}\left \langle (\Delta \mu )^{2} \right \rangle \tag{23.19}\]
This is a consequence of the Taylor-Green-Kubo formula (Taylor 1922; Green 1951; Kubo 1957), \[ D_{\mu\mu} = \int_{0}^{\infty} \mathrm{d} t \langle \dot{\mu}(0) \dot{\mu}(t) \rangle \tag{23.20}\] where the dots denote derivatives with respect to time. For diffusive transport, Equation 23.14 and Equation 23.19 coincide, of course. Moreover, this allows computing the parallel diffusion coefficient \(\kappa_\parallel\) without the detour of computing \(D_{\mu\mu}\) first and then applying the diffusion approximation, Equation 23.17.
From the equations of motion Equation 23.1, we find (HOW???) \[ \dot{\mu} = \frac{e}{c p} \left( \mathbf{v} \times \mathbf{B} \right)_{z} = \frac{1}{r_{\text{g}} B_z} \left(v_{x} \delta B_{x}(\mathbf{r}) - v_{y} \delta B_{y}(\mathbf{r}) \right) \tag{23.21}\] and thus \[ D_{\mu\mu} = \frac{1}{ B_{z} ^{2} r_{\text{g}}^{2}} \int_{0}^{\infty }\mathrm{d} t \left[ v_{x}(t) v_{x}(0) \mathcal{P}_{yy}(t) + v_{y}(t) v_{y}(0) \mathcal{P}_{xx}(t) \right] \tag{23.22}\]
Here, we have defined \[ \mathcal{P}_{ij}(t) \equiv \langle \delta B_{i}(0) \delta B_{j}(t) \rangle \tag{23.23}\] and both the velocities and the magnetic fields are to be evaluated along unperturbed trajectories. Note that the fact that the Fokker-Planck coefficients only depend on the two-point function means that we can constrain ourselves to the Gaussian part of the turbulent magnetic field.
23.1.4 Turbulence geometries and spectra
To make further progress, we need to specify the turbulence correlation tensor \(\mathcal{P}_{ij}\). In the derivation of the Fokker-Planck equation we had to assume that turbulence is homogeneous and stationary (assumption 5). In this case, the field can be represented very economically in Fourier space. To this end, we introduce the Fourier transform pair \[ \begin{aligned} \delta \tilde{B}_{j}(\mathbf{k}, t) &= (2\pi)^{-3/2} \int_{-\infty }^{\infty} \mathrm{d}^{3} r \, \delta B_{j}(\mathbf{r}, t) \mathrm{e} ^{\imath \mathbf{k} \cdot \mathbf{r}} \\ \delta B_{j}(\mathbf{x}, t) &= (2\pi)^{-3/2} \int_{-\infty}^{\infty} \mathrm{d}^{3} k \, \delta \tilde{B}_{j}(\mathbf{k}, t) \mathrm{e} ^{-\imath \mathbf{k} \cdot \mathbf{r}} \end{aligned} \tag{23.24}\]
Note that for the magnetic field to have real values, \(\delta B_{j}(\mathbf{r}) = \delta B^{*}_{j}(\mathbf{r})\), requires a relation between the Fourier components and their complex conjugates, \[ \delta \tilde{B}^{*}_{j}(\mathbf{k}) = \delta \tilde{B}_{j}(-\mathbf{k}) \]
The homogeneity and stationarity now guarantee that the two-point functions \(\langle \delta B_{i}(\mathbf{r}, t) \delta B_{j}(\mathbf{r}', t') \rangle\) depend on the positions \(\mathbf{r}\) and \(\mathbf{r}'\) and times \(t\) and \(t'\) only through the differences \(\Delta \mathbf{r} \equiv (\mathbf{r} - \mathbf{r}')\) and \((t-t')\). It is then easy to see that the two-point function in Fourier space is diagonal, \[ \begin{aligned} & \langle \delta \tilde{B}_{i}(\mathbf{k}, t) \delta \tilde{B}^{*}_{j}(\mathbf{k}', t') \rangle \\ &= \! (2 \pi )^{-3} \!\!\! \int \!\! \mathrm {d} ^{3} r \, \mathrm {e} ^{ \imath \mathbf{k} \cdot \mathbf{r}} \!\!\! \int \!\! \mathrm {d} ^{3} r' \mathrm {e} ^{-\imath \mathbf{k}' \cdot \mathbf{r}'} \!\! \langle \delta B_{i}(\mathbf{r}, t) \delta B_{j}(\mathbf{r}'\!, t') \rangle \\ &= \delta ^{(3)} (\mathbf{k} - \mathbf{k}') P_{ij}(\mathbf{k}, t-t') \end{aligned} \] where the turbulence correlation tensor \(P_{ij}(\mathbf{k}, \Delta t)\) is the Fourier transform of the coordinate space two-point function \[ \begin{aligned} P_{ij}(\mathbf{k}, \Delta t) \equiv & (2 \pi )^{-3/2} \int _{-\infty}^{\infty} \mathrm {d} ^{3} (\Delta r) \, \mathrm {e} ^{\imath \mathbf{k} \cdot \Delta \mathbf{r}} \\ & \times \langle \delta B_{i}(\mathbf{r}, t) \delta B_{j}(\mathbf{r} - \Delta \mathbf{r}\!, t') \rangle \end{aligned} \tag{23.25}\]
It contains all the (statistical) information on the magnetic turbulence that enters into the computation of the Fokker-Planck coefficients:
- information on the turbulence geometry, for instance whether there is a preferred direction for the propagation of waves;
- information on the turbulence spectrum, that is the distribution of energy among different turbulent scales;
- information on the time-dependence of the correlations.
Oftentimes, it is assumed that \(P_{ij}(\mathbf{k}, \Delta t)\) factorises into a magnetostatic correlation tensor \(P_{ij}(\mathbf{k}) \equiv P_{ij}(\mathbf{k}, 0)\) independent of time and a time-dependent dynamical correlation function \(\Gamma(\mathbf{k}, \Delta t)\), \[ P_{ij}(\mathbf{k}, \Delta t) = P_{ij}(\mathbf{k}) \Gamma (\mathbf{k}, \Delta t) \tag{23.26}\]
In the magnetostatic approximation, we ignore any time-dependence altogether, so \[ \Gamma \equiv 1 \]
While in reality \(P_{ij}\) may be arbitrarily complicated, three turbulence geometries have dominated much of the literature, both in analytical studies of transport coefficients and numerical test particle simulations. These three geometries are conceptually simple and particularly amenable to analytical computations of the components of the diffusion tensor and the other Fokker-Planck coefficients:
- 3D isotropic turbulence
- slab turbulence
- a composition of slab and 2D isotropic turbulence
23.1.4.1 3D isotropic turbulence
For 3D isotropic turbulence the magnetostatic correlation tensor takes the form (Batchelor 1982) \[ P^{\text{3D}}_{ij}(\mathbf{k}) = g^{\text{3D}}(k) \left( \delta _{ij} - \frac{k_{i} k_{j}}{k^{2}} + \imath \sigma(k) \epsilon_{ijm} \frac{k_{m}}{k} \right) \tag{23.27}\]
The k-dependent real functions \(g^{\text{3D}}(k)\) and \(\sigma(k)\) allow modelling of the overall spectrum and of a wavenumber-dependent helicity, respectively. Note that for linearly polarised waves \(\sigma(k)\equiv 0\). The normalization of \(g^{\text{3D}}(k)\) is fixed by requiring \[ \begin{aligned} \delta B^{2} &\equiv \langle \mathbf{\delta B}^{2}(\mathbf{x}) \rangle \! = \!\! \int \mathrm {d} ^{3} k \left ( P^{\text{3D}}_{xx}(\mathbf{k}) + P^{\text{3D}}_{yy}( \mathbf{k}) + P^{\text{3D}}_{zz}(\mathbf{k}) \right ) \\ &= 2 \int \mathrm {d} ^{3} k \, g^{\text{3D}}(k) = 8 \pi \int _{0}^{\infty } \mathrm {d} k \, k^{2} g^{\text{3D}}(k) \end{aligned} \tag{23.28}\]
23.1.4.2 Slab turbulence
In slab turbulence, it is assumed that all quantities are independent of the coordinates perpendicular to the background field (x and y) and that the turbulent field has no z-component. Consequently, the wave vectors \(\mathbf{k}\parallel \hat{z}\) and if we further demand turbulence to be axisymmetric, the turbulence correlation tensor reads \[ P^{\text{slab}}_{ij}(\mathbf{k}) = g^{\text{slab}}(k_{\parallel }) \frac{\delta (k_{\perp })}{k_{\perp }} \left ( \delta _{ij} + \imath \sigma (k_{\parallel }) \epsilon _{ijz} \right ) \tag{23.29}\] for \(i, j \in x, y\) and zero otherwise. In our case \(k_\parallel = k_z\) and \(k_{\perp } = \sqrt{k_{x}^{2} + k_{y}^{2}}\). Again, \(\sigma(k_\parallel)\) allows for wavenumber-dependent helicity, but vanishes for linear polarisation. The normalisation is then \[ \begin{aligned} \delta B^{2} &\equiv \int \mathrm{d} ^{3} k \left( P^{\text{slab}}_{xx}( \mathbf{k}) + P^{\text{slab}}_{yy}(\mathbf{k}) \right) \\ &= 4 \pi \int _{-\infty}^{\infty} \mathrm{d} k_{\parallel} \, g^{ \text{slab}}(k_{\parallel}) \end{aligned} \tag{23.30}\]
While slab turbulence might seem rather restrictive a turbulence model, it is quite attractive due to its simplicity. In addition, it could be argued that it is of physical relevance in situations where the turbulence is self-generated by anisotropies in the distribution of CRs (Kulsrud and Pearce 1969; Skilling 1975): It has been shown (e.g. Tademaru 1969) that the modes with wavevectors along the background magnetic field grow fastest.
23.1.4.3 Composite (slab + 2D isotropic) turbulence
23.1.4.4 Turbulence spectra
Having reviewed three simple turbulence geometries, we need to specify the spectral shapes \(g(k)\) in order to compute transport coefficients. In cascade models of turbulence, energy is injected on the largest scales in the so-called injection (energy) range. Non-linear interactions transfer energy to smaller scales over the so-called inertial range. At very small scales, the turbulent energy is dissipated in the so-called dissipation range. The scale separating the injection and inertial ranges is called the outer scale of turbulence and the scale separating the inertial and the dissipation range is called the dissipation scale.
Most studies deal with one of two spectra. The first one is a simple power law with spectral index \(q\) and low wavenumber cut-off \(k_0\), corresponding to the outer scale \((2\pi/k_0)\), \[ g_{\text{PL}}(k) = \left \{ \textstyle\begin{array}{l@{\quad }l} g_{0} (k/k_{0})^{-q} & \text{for } k \geq k_{0} \, , \\ 0 & \text{otherwise.} \end{array}\displaystyle \right . \]
The alternative is a broken power law with a flat spectrum below the wavenumber \(k_0\) and a power law slope \(q\) above, \[ g_{\text{BPL}}(k) = g_{0} \left ( 1 + \left ( \frac{k}{k_{0}} \right )^{1/s} \right )^{-q s} \, . \] Here, \(s\) is parametrising the softness of the break and \(s\to 0\) corresponds to a sharp break. It is assumed that the broken power law form can potentially also capture turbulence in the injection range, that is for \(k<k_0\).
23.1.4.5 Slab turbulence with broken power law spectrum
The result for the pitch-angle diffusion coefficient in slab turbulence and for the broken power law spectrum (Shalchi 2009) is given as an example. \[ g^{\text{slab}}(k_{\parallel }) = \frac{C \left (q, \frac{1}{2} \right )}{2 \pi k_{0}} \delta B^{2} \left ( 1 + \left ( \frac{k}{k_{0}} \right )^{2} \right )^{-q/2} \, . \tag{23.31}\]
The function \(C(q, s)\) is fixed by the normalisation condition, see Equation 23.30, \[ \begin{aligned} \frac{1}{C(q, s)} & \equiv \frac{4}{k_{0}} \int_{0}^{\infty} \mathrm {d} k \, \left ( 1 + \left( \frac{k}{k_{0}} \right )^{1/s} \right)^{-q s} \\ &= \frac{4}{k_{0}} \frac{\Gamma (s(q-1)) \Gamma (1+s)}{\Gamma (q \, s)} \end{aligned} \tag{23.32}\]
where \(\Gamma (\cdot)\) denotes the gamma function. We have assumed that \(q>1\) in order for the k-integral in Equation 23.32 to converge. For \(q=1\), instead, we need to assume a cut-off, that is \(g^{\text{slab}} = 0\) for \(k_{\parallel} > k_{\text{max}}\).
Substituting Equation 23.31 into Equation 23.29 and Equation 23.22, one encounters the resonance function (Schlickeiser 2002) \[ R^{\text{slab}} = \pi \delta (k_{\parallel } \mu v \pm \Omega ) \tag{23.33}\]
Eventually, this simplifies to \[ D_{\mu \mu } = \frac{\pi }{2} C \! \left ( \! q, \frac{1}{2} \! \right ) q k_{0} \frac{\delta B^{2}}{B_{z}^{2}} \frac{(1 - \mu ^{2}) \mu ^{q-1} (r_{\text{L}} k_{0})^{q-2}}{(1 + \mu ^{2} (r_{\text{L}} k_{0}))^{q/2}} \tag{23.34}\]
23.1.5 Field-line random walk
The computation of the pitch-angle diffusion coefficient in Equation 23.20 and Equation 23.22 is based on an evaluation of the turbulent part of the Lorentz force along trajectories around the homogeneous background field. As long as perturbations are small, this gives the dominant contribution to the parallel diffusion coefficient, Equation 23.17.
For perpendicular transport, however, there is another important contribution due to the fact that the field line is not perfectly homogeneous. Instead, the large-scale magnetic field evaluated for a particle along a field line changes direction with distance along this field line. Under certain conditions, this movement can be shown to be diffusive, see below. If the movement of the particle due to this effect is included in the computation of the mean-square displacements (or equivalently through the Taylor-Green-Kubo approach), this gives the so-called field-line random walk (FLRW) contribution to perpendicular transport. The contribution without this is oftentimes called the microscopic contribution.
For slab turbulence, the microscopic diffusion coefficient vanishes (the transport is in fact sub-diffusive), hence FLRW gives the only contribution. For other turbulence geometries, FLRW can also contribute, but might not be dominating.
Let’s again assume the regular background field \(\langle \mathbf{B} \rangle = B_{z} \hat{z}\) to be dominating over the perturbations \(\delta\mathbf{B}\). The equation determining the field line \(\{ x(z), y(z) \}\) is \[ \frac{\mathrm{d}x}{\mathrm{d}z} = \frac{\delta B_{x}}{B_{z}} \, ,\frac{\mathrm{d}y}{\mathrm{d}z} = \frac{\delta B_{y}}{B_{z}} \] This can formally be integrated to obtain the mean square displacement in the perpendicular directions, e.g. \[ \langle (\Delta x)^{2} \rangle \! = \! \frac{1}{B_{z}^{2}} \int _{0}^{z} \!\! \! \mathrm {d} z' \!\! \int _{0}^{z} \!\!\! \mathrm {d} z'' \langle \delta B_{x}(\mathbf{r}(z')) \delta B_{x}(\mathbf{r}(z'')) \rangle \tag{23.35}\]
In slab turbulence, the integrand only depends on z and it is easy to show that the perpendicular mean-square displacement \(\langle (\Delta r_{\perp})^{2} \rangle\) is ballistic at small z and diffusive for large z, e.g. \[ \langle (\Delta x)^{2} \rangle = \left \{ \textstyle\begin{array}{l@{\quad } l} z^{2} \left ( \delta B_{x} / B_{z} \right )^{2} & \text{for } z \to 0 \, , \\ 2 \kappa _{\text{FLRW}} |z| & \text{for } z \to \infty \, , \end{array}\displaystyle \right . \] with the FLRW diffusion coefficient \[ \kappa _{\text{FLRW}} = \frac{2 \pi ^{2}}{B_{z}^{2}} g^{\text{slab}}(0) \]
In other turbulence geometries, the integrand in Equation 23.35 also depends on x and y, such that an explicit solution is not possible without further assumptions.
If particles are assumed to diffuse along field lines, \(\langle (\Delta z)^{2} \rangle \propto \Delta t\), FLRW leads to subdiffusive perpendicular transport, \(\langle (\Delta r_{\perp})^{2} \rangle \propto \sqrt{\Delta t}\), a phenomenon known as compound (sub)diffusion (Jokipii 1966). Compound subdiffusion has been applied to a variety of environments like laboratory plasmas, the heliosphere, near-source transport as well as shock acceleration.
23.1.6 Short-comings of QLT
Despite its popularity, QLT exhibits a number of issues. The most well-known pathology of magnetostatic QLT is its inability to scatter particles through \(90^\circ\). While present in a number of turbulence geometries, it is easiest illustrated in slab turbulence where the dependence of the pitch-angle diffusion coefficient \(D_{\mu\mu}\) on the spectrum \(g^{\text{slab}}(k)\) becomes very simple. In fact, inspecting Equation 23.34 we see that \(D_{\mu\mu}\to 0\) for \(\mu\to 0\).
The root cause for the \(90^\circ\) problem is the narrow resonance condition in magnetostatic QLT, \(k_{\parallel } \mu \, r_{\text{L}} = \pm 1\) (Equation 23.33). Particles at finite \(\mu\) are in resonance with waves of finite parallel wavenumber, \(k_{\parallel} = \pm 1 / (\mu r_{\text{L}})\). For \(\mu\to 0\), however, the resonant parallel wavenumber grows without bounds. With the turbulence spectra being falling power laws, however, there is only little energy at small scales and the pitch-angle scattering rate vanishes. In practice, there is of course no energy at all at scales below the dissipation scale.
23.2 Test Particle Simulation
In the absence of a widely accepted extension of quasi-linear theory, wave-particle interactions must also be studied in numerical simulations where the equations of motion are directly solved in a realisation of the turbulent magnetic field. The applications of such test particle simulations of cosmic rays are manifold: testing transport theories, computing parameters like diffusion coefficients or making predictions for phenomena beyond standard diffusion theories, e.g. for cosmic ray small-scale anisotropies.
Ideally, one would test these non-linear theories beyond QLT by comparing their predictions with data from observations. While this approach is being followed by the heliospheric community, a difficulty remains in that the actual turbulence (if known) turns out to be much more complex than what is routinely assumed in analytical transport theories. Alternatively, transport theories can be tested by comparing their predictions with those from numerical experiments.
Test particle simulations compute the transport of high-energy charged particles through prescribed EM fields without taking into account the effect of the high-energy particles on the EM fields. To this end, a realisation of the turbulent magnetic field is generated and the equations of motion are solved for test particles, that is the contributions of the CR particles to the EM fields are ignored. Given the trajectories of a large enough number of test particles, one can numerically compute the Fokker-Planck coefficients.
In addition to testing transport theories by comparing the analytically computed diffusion coefficients to simulated ones, there are at least two more applications of test particle simulations:
For sources at distances closer or similar to the scattering mean free path, the diffusive transport theory is not necessarily applicable. An often-cited pathology of computing solutions to the diffusion equation is superluminal propagation speeds. Lately, there has been increased interest in the transition between the ballistic and diffusive phases of transport and test particle simulations allow exploring this transition for given turbulent EM fields.
Analytical transport theories usually make predictions only for the ensemble-averaged phase-space density and it is usually assumed that the observed phase-space densities are close to the ensemble average. This has now been called into question, in particular in view of the observation of small-scale anisotropies observed in the arrival directions of TeV-PeV CRs. Test particle simulations naturally simulate CR distributions for individual realisations of the turbulent fields and thus provide direct access to such stochasticity effects.
In CGS units, the electric and magnetic fields have the same unit.↩︎
Strictly speaking, the correlation lengths and times are tensors because of the vector nature of the magnetic field in the two-point functions; here, however, we only require them for order of magnitude arguments, so we do not distinguish the different components.↩︎
In particle transport, ballistic refers to the movement of particles with minimal influence from their environment. Imagine a particle traveling on a frictionless, straight path, almost like a bullet. That’s ballistic transport!↩︎