2.1. Cyclotron cooling

The emission of cyclotron radiation by particles in a plasma gives rise to a reaction force which must be included in the kinetic equation. The theory of relativistic plasmas, accounting for this radiation reaction by including the Abraham–Lorentz reaction force in the kinetic equation, has been developed by Andersson, Helander & Eriksson (Reference Andersson, Helander and Eriksson2001), Hazeltine & Mahajan (Reference Hazeltine and Mahajan2004) and subsequent authors. Here, we begin by specialising these results to a non-relativistic plasma.

The non-relativistic Abraham–Lorentz reaction force is given by

(2.1)\begin{equation} \boldsymbol{K} = \frac{{e}^{2}}{6{\rm \pi}\epsilon_{0}c^{3}}\dot{\boldsymbol{a}}, \end{equation}

where $\boldsymbol {a}$ is the acceleration vector and $e$ is the charge. Following the arguments of the aforementioned authors, it can be seen that, to leading order in a small-gyroradius expansion, the change in perpendicular energy, $w_{\perp }$, of a non-relativistic point charge as it accelerates in a magnetic field is given by Larmor's formula

(2.2a,b)\begin{equation} \left(\frac{{\mathrm{d}} w_{{\perp}}}{{\mathrm{d}} t}\right)_{\text{rad}} ={-} \frac{e^{2}a^{2}}{6{\rm \pi} \epsilon_{0}c^{3}}, \quad w_{{\perp}} = \frac{mv_{{\perp}}^{2}}{2}, \end{equation}

where $\boldsymbol {v}$ is the particle velocity vector and $m$ is the rest mass. Throughout this work, the subscripts $\perp$ and $\parallel$ indicate directions perpendicular and parallel to the magnetic field, respectively. The centripetal acceleration of a particle undergoing gyromotion is given by

(2.3)\begin{equation} a = \frac{eB}{m}v_{{\perp}}, \end{equation}

and hence it follows that

(2.4)\begin{equation} \left(\frac{{\mathrm{d}} w_{{\perp}}}{{\mathrm{d}} t}\right)_{\text{rad}} =- \frac{e^{4}B^{2}}{3{\rm \pi}\epsilon_{0}(mc)^{3}} w_{{\perp}}. \end{equation}

Simply put, the plasma will radiate its perpendicular kinetic energy on a time scale given by the radiation time:

(2.5)\begin{equation} \tau_{r} = \frac{3{\rm \pi}\epsilon_{0}(mc)^{3}}{e^{4}B^{2}}. \end{equation}

It is perhaps easiest to understand the influence of radiative emission when the change of energy is rewritten in terms of the perpendicular and parallel components of the velocity. We thus obtain

(2.6a,b)\begin{equation} \left(\frac{{\mathrm{d}} v_{{\perp}}}{{\mathrm{d}} t}\right)_{\text{rad}} ={-}\frac{v_{{\perp}}}{2\tau_{r}}, \quad \left(\frac{{\mathrm{d}} v_{{\parallel}}}{{\mathrm{d}} t}\right)_{\text{rad}} = 0. \end{equation}

Equations (2.6a,b) summarise the average effect of radiation reaction on the plasma. Fundamentally, it is these equations which will give rise to anisotropy in the distribution function.

It is pertinent to point out at this stage that collisions, whose effect is to isotropise the plasma, can mediate this cooling process due to the scattering of the velocity vector, as a result allowing the conversion of energy between perpendicular and parallel components. We will return to this discussion in § 3, following a discourse on the simpler collisionless case.

2.2. Collisionless kinetic equation

The kinetic equation governing the evolution of the distribution function $f(\boldsymbol {r},\boldsymbol {v},t)$ in the absence of collisions is given by

(2.7)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{v}f) + \boldsymbol{\nabla}_{v}\boldsymbol{\cdot} \left[\frac{q}{m}(\boldsymbol{E} + \boldsymbol{v}\times\boldsymbol{B})\,f + \left(\frac{{\mathrm{d}} \boldsymbol{v}}{{\mathrm{d}} t}\right)_{\text{rad}}f\right] = 0 , \end{equation}

where $\boldsymbol {\nabla }_{v}$ stands for the velocity derivative:

(2.8)\begin{equation} \boldsymbol{\nabla}_{v} = \hat{{\boldsymbol{x}}}\frac{\partial }{\partial v_{x}} + \hat{{\boldsymbol{y}}}\frac{\partial }{\partial v_{y}} + \hat{{\boldsymbol{z}}}\frac{\partial }{\partial v_{z}}. \end{equation}

Cyclotron emission is included in the kinetic equation through the inclusion of (2.6a,b) in the usual collisionless Vlasov equation.

We will immediately specialise to the case where the magnetic field is constant, $\boldsymbol {B} = B\hat {{\boldsymbol {b}}}$, and there is zero electric field, $\boldsymbol {E} = \boldsymbol {0}$. Note that, in this limit, we have assumed the plasma is homogeneous and there is no spatial dependence of the distribution function. We then write (2.7) in cylindrical $\boldsymbol {v}$-space coordinates $(v_{\perp },\alpha ,v_{\parallel })$ to obtain

(2.9)\begin{equation} \frac{\partial f}{\partial t} + \frac{1}{v_{{\perp}}}\frac{\partial }{\partial v_{{\perp}}}(v_{{\perp}}\dot{v}_{{\perp}}\,f) + \frac{\partial }{\partial \alpha} (\dot{\alpha}f) = 0, \end{equation}

where we have used that $\dot {{v}}_{\parallel } = 0$. The final term in this equation can then be eliminated by averaging over the gyroangle $\alpha$.

It is fruitful to write the resulting equation in terms of the perpendicular energy of the particles, so as to obtain

(2.10)\begin{equation} \frac{\partial f}{\partial t}- \frac{w_{{\perp}}}{\tau_{r}}\frac{\partial f}{\partial w_{{\perp}}} = \frac{f}{\tau_{r}}. \end{equation}

This equation can then be easily solved via the method of characteristics; the general solution is given by

(2.11)\begin{equation} f(w_{{\perp}},w_{{\parallel}},t) = F_{0}(w_{{\perp}}\,\mathrm{e}^{t/\tau_{r}},w_{{\parallel}}) \mathrm{e}^{t/\tau_{r}}. \end{equation}

That is, the distribution function will be a function of the parallel and perpendicular kinetic energies, the latter of which will decay exponentially quickly on the radiation time scale.

As an illustrative example, if the initial distribution is Maxwellian, then

(2.12)\begin{equation} f(w_{{\perp}},w_{{\parallel}},t) = n \left(\frac{m}{2{\rm \pi} T}\right)^{3/2} \exp \left[ - \frac{1}{T} (w_{{\perp}}\,\mathrm{e}^{t/\tau_{r}} + w_{{\parallel}}) + \frac{t}{\tau_{r}}\right], \end{equation}

which corresponds to a bi-Maxwellian

(2.13)\begin{equation} f(w_{{\perp}},w_{{\parallel}},t) = n \left(\frac{m}{2{\rm \pi} T_{{\perp}}^{2/3}T_{{\parallel}}^{1/3}}\right)^{3/2} \exp \left(-\frac{w_{{\perp}}}{T_{{\perp}}} - \frac{w_{{\parallel}}}{T_{{\parallel}}}\right) \end{equation}

with $T_{\parallel } = T = \text {const.}$ and

(2.14)\begin{equation} T_{{\perp}}(t) = T\, \mathrm{e}^{{-}t/\tau_{r}}. \end{equation}

At this stage of our investigation, we have been able to accomplish the first fundamental aim of this paper; namely, we now understand the following. Firstly, in the most simple case, cyclotron emission is responsible for an exponential decay of the perpendicular energy on the radiation time scale whilst the parallel energy is kept constant. Secondly, the inclusion of radiative emission, and the absence of collisions, results in a strongly anisotropic distribution function on the radiation time scale.

2.3. Validity of the collisionless approach

A caveat here is that, of course, many plasmas are in fact collisional. It is therefore expected that eventually collisions will come into play and must be taken into account.

Following Braginskii (Reference Braginskii1965), the conventional definitions of the electron–ion and electron–electron collision times are given by

(2.15a,b)\begin{equation} \tau_{ei} = \frac{6\sqrt{2}{\rm \pi}^{3/2}\epsilon_{0}^{2}{m_{e}^{1/2}}T_{e}^{3/2}}{ Ze^{4} n_{e} \ln \varLambda}, \quad \tau_{ee} = Z\tau_{ei} \end{equation}

respectively, where $\ln \varLambda$ is the Coulomb logarithm. We will also assume bulk neutrality throughout this work so that $n_{e} = Zn_{i}$. From the above expressions, we can see that the collision time will decrease as the plasma cools.

We might therefore envision the following scenario. A plasma might begin in a regime where the radiation time $\tau _{r}$ is smaller than the initial collision time $\tau _{c} = \text {min}(\tau _{ee},\tau _{ei})$ and hence the collisionless theory will be applicable at first. However, the collision time itself will decrease as the plasma cools and hence this assumption will be violated after sufficient time has elapsed. As a result, collisional effects must be taken into consideration.

We thus imminently turn our attention to the collisional problem. However, before doing so, let us allow ourselves a small digression to say something about the nature of cyclotron radiation.

2.4. Cyclotron cooling and plasma waves

Before we embark on solving the collisional kinetic equation, let us first turn our attention to the nature of the electromagnetic waves in play. The cyclotron radiation that is being generated has to propagate through our optically thin plasma as some electromagnetic plasma wave (Krall & Trivelpiece Reference Krall and Trivelpiece1986).

For the reasons outlined in the previous section, we declare an interest in optically thin plasmas where the radiation time is (at least initially) much shorter than the collision time. That is, we declare an interest in the limit

(2.16)\begin{equation} \frac{\tau_{r}}{\tau_{c}} = \frac{mn\ln \varLambda}{{\rm \pi} \epsilon_{0}B^{2}}\left(\frac{c}{v_{\text{th}}}\right)^{3} = \frac{\ln \varLambda}{\rm \pi} \left(\frac{\omega_{p}}{\varOmega_{c}}\right)^{2}\left(\frac{c}{v_{\text{th}}}\right)^{3} \ll 1, \end{equation}

where have introduced the two characteristic frequencies in the plasma, namely the plasma frequency and the cyclotron frequency:

(2.17a,b)\begin{equation} \omega_{p}^{2} = \frac{ne^{2}}{m\epsilon_{0}}, \quad \varOmega_{c} = \frac{eB}{m}. \end{equation}

Thus, it is easy to see from this ordering that electromagnetic waves with frequencies comparable to the cyclotron frequency are almost light waves:

(2.18)\begin{equation} \omega^{2} ({\sim} \varOmega_{c}^{2}) = \omega_{p}^{2} + k^{2}c^{2} \implies \omega \simeq kc. \end{equation}

That is, our cyclotron radiation is simply photon emission which can immediately escape the plasma because of the optical thinness.

We now turn our attention to the collisional problem.

5.1. Limit (A): neglecting the contribution from electron–electron collisions

In limit (A), we are justified in neglecting the contribution of $I_{ee}$ whilst retaining that of $I_{ei}$, noting that this is still consistent with the derivation of (3.19) providing the required ordering is satisfied. In this instance, (4.5a–c) becomes

(5.5)\begin{equation} \frac{\partial N}{\partial \tau} = \frac{\partial}{\partial v_{{\parallel}}}\left(\frac{N}{v_{{\parallel}}^{2}}\right). \end{equation}

The equation has been reduced to a quasilinear partial differential equation and so is amenable to solution by the method of characteristics. We thus obtain the general solution

(5.6)\begin{equation} N(v_{{\parallel}},\tau) = v_{{\parallel}}^{2}F\left(\frac{v_{{\parallel}}^{3}}{3} + \tau\right). \end{equation}

So, if we start from a Maxwellian distribution initially

(5.7)\begin{equation} N(v_{{\parallel}},0) = \left(\frac{m}{2{\rm \pi} T_{{\parallel}}}\right)^{1/2}\exp\left(-\frac{mv_{{\parallel}}^{2}}{2T_{{\parallel}}}\right), \end{equation}

then

(5.8)\begin{equation} N(v_{{\parallel}},\tau) = \left(\frac{m}{2{\rm \pi} T_{{\parallel}}}\right)^{1/2} \frac{v_{{\parallel}}^{2}}{(v_{{\parallel}}^{3} + 3\tau)^{2/3}} \exp\left(-\frac{m}{2T_{{\parallel}}}(v_{{\parallel}}^{3} + 3\tau)^{2/3}\right). \end{equation}

From this equation, we can also obtain an expression for the perpendicular temperature:

(5.9)\begin{equation} T_{{\perp}}(v_{{\parallel}},\tau) = 4\sigma \tau_{r}m|\ln \epsilon|\left(\frac{m}{2{\rm \pi} T_{{\parallel}}}\right)^{1/2} \frac{v_{{\parallel}}^{2}}{(v_{{\parallel}}^{3} + 3\tau)^{2/3}} \exp\left(-\frac{m}{2T_{{\parallel}}}(v_{{\parallel}}^{3} + 3\tau)^{2/3}\right). \end{equation}

It is important to recall that this solution was for $v_{\parallel } > 0$, which means that, despite appearances, both $T_{\perp }$ and $N$ are continuous. For $v_{\parallel } < 0$ we obtain the same solution as above with $v_{\parallel } \rightarrow -v_{\parallel }$.

5.1.1. Rate of energy loss

A quantity of direct experimental interest is the amount of power radiated by the plasma. This can be calculated by looking at the rate of change of energy loss:

(5.10)\begin{equation} P ={-}\frac{{\mathrm{d}} W}{{\mathrm{d}} t}, \end{equation}

where we have defined the total thermal energy stored in the plasma (per unit volume):

(5.11)\begin{equation} W = \int \tfrac{1}{2} mnv^{2} f \, {\mathrm{d}}^{3}\boldsymbol{v}. \end{equation}

The energy stored can be seen to satisfy the equation

(5.12a,b)\begin{equation} \frac{{\mathrm{d}} W}{{\mathrm{d}} t} ={-}\frac{W_{{\perp}}}{\tau_{r}}, \quad W_{{\perp}} = \int \tfrac{1}{2}mnv_{{\perp}}^{2} f \, {\mathrm{d}}^{3}\boldsymbol{v}. \end{equation}

Thus, it suffices to calculate $W_{\perp }$ in order to find the power radiated by the plasma. This is given by evaluating

(5.13)\begin{align} W_{{\perp}} &= \int \frac{1}{2}mnv_{{\perp}}^{2} f \, {\mathrm{d}}^{3} \boldsymbol{v} \nonumber\\ &= \int \frac{mn}{8\sigma\tau_{r}|\ln\epsilon|} v_{{\perp}}^{3} \exp \left(-\frac{mv_{{\perp}}^{2}}{2T_{{\perp}}(v_{{\parallel}},t)}\right) {\mathrm{d}} v_{{\parallel}} \,{\mathrm{d}} v_{{\perp}} \nonumber\\ &= \int \frac{nT_{{\perp}}(v_{{\parallel}},t)^{2}}{4\sigma\tau_{r}m|\ln\epsilon|} \, {\mathrm{d}}v_{{\parallel}}, \end{align}

where we have used (3.19) to carry out the $v_{\perp }$ integral. We arrive at

(5.14)\begin{equation} \frac{W_{{\perp}}(y)}{W_{{\perp}}(0)} = \frac{2}{\sqrt{\rm \pi}}\int_{0}^{\infty} \frac{x^{4}}{(x^{3}+y)^{4/3}}\exp (-(x^{3} + y)^{2/3})\, {\mathrm{d}} x, \end{equation}

where we have normalised $v_{\parallel }$ and $\tau$ to the electron thermal velocity by writing

(5.15a–c)\begin{equation} v_{{\parallel}} = v_{\text{th}}x, \quad 3\tau = v_{\text{th}}^{3}y, \quad v_{\text{th}} = \sqrt{\frac{2T_{{\parallel}}}{m}}. \end{equation}

The integral (5.14) cannot be expressed in terms of standard mathematical functions. Instead, we must turn to numerically evaluating this function. Firstly, however, we can determine the behaviour of (5.14) in the long-time and short-time limits.

In the limit $y \gg 1$ we obtain (see appendix 1)

(5.16)\begin{equation} \frac{W_{{\perp}}(y)}{W_{{\perp}}(0)} \simeq \left(\frac{2}{3{\rm \pi}}\right)^{8/9} \varGamma \left(\frac{5}{3}\right) \left(\frac{t}{\tau_{ei}}\right)^{{-}7/9} \exp\left( -\left(\frac{9\sqrt{\rm \pi}}{4}\frac{t}{\tau_{ei}}\right)^{2/3}\right). \end{equation}

In the limit $y \ll 1$ we obtain (see appendix 2)

(5.17)\begin{equation} \frac{W_{{\perp}}(t)}{W_{{\perp}}(0)} \simeq 1 - \frac{\varGamma(5/3) \varGamma(2/3)}{\varGamma(7/3)} \left(\frac{128}{3{\rm \pi}}\right)^{1/6} \left(\frac{t}{\tau_{ei}}\right)^{1/3}. \end{equation}

In figure 1 we show the full numerical solution of (5.14) as well as the analytic solution for the long-term and short-term behaviour, given by (5.16) and (5.17), respectively. Some caution must be taken in interpreting this figure, and one must remember that the solution of (5.14) is only valid on time scales longer than the radiation time. Note that the quantity plotted on the vertical axis, $W_{\perp }(y)/W_{\perp }(0),$ is directly proportional to the emitted power as can be seen from (5.12a,b). The abscissa is simply a scaled time coordinate: $y = (9\sqrt {{\rm \pi} }/4)t/\tau _{ei}$.

Figure 1.

The ratio of perpendicular energy at normalised time $y = (9\sqrt {{\rm \pi} }/4)t/\tau _{ei}$ to the initial perpendicular energy, a quantity that is proportional to the rate of energy loss in the plasma. The blue curve is the calculated rate of energy loss from (5.14). The red curve is the asymptotic solution for $y \gg 1$ given by (5.16). The purple curve is the asymptotic solution for $y \ll 1$ given by (5.17). This figure is for the perpendicular energy loss in limit (A).

Although at a first glance it may seem as though the energy loss rate would be faster than exponential, this is of course not the case when the solution is restricted to the region where the orderings are valid.

5.2. Limit (B): neglecting the contribution from electron–ion collisions

In limit (B), when $n_{i} Z^{2} / n_{e} \ll | \ln \epsilon |$, we are justified in neglecting the contribution from $I_{ei}$ and obtain a relatively simple equation:

(5.18)\begin{equation} \frac{\partial N}{\partial \tau} ={-}\alpha \frac{\partial ^{2} N^{2}}{\partial v_{{\parallel}}^{2}}. \end{equation}

This is a nonlinear parabolic partial differential equation which essentially describes backwards diffusion in $v_{\parallel }$ space. Indeed, this is simply the backwards heat equation in a medium where the diffusion coefficient varies linearly with temperature. As such, this problem is ill-posed.

A simple solution exists (given suitable initial data) but the general solution will not depend continuously on the auxiliary data as it describes a reverse diffusion process meaning that arbitrarily small perturbations in the initial conditions will be amplified and can lead to even infinitesimal-time singularities. This means that it is possible to prepare arbitrarily small perturbations to the initial conditions (likely in any given Sobolev norm), which are sufficiently wiggly in higher derivatives to ensure the solution blows up before any given time $t > 0$.

In order to find this singular solution, we consider the problem

(5.19a,b)\begin{equation} \frac{\partial N}{\partial \tau} - \beta \frac{\partial}{\partial v_{{\parallel}}}\left(N\frac{\partial N}{\partial v_{{\parallel}}}\right) = 0, \quad \beta ={-} 2\alpha \end{equation}

with

(5.20)\begin{equation} \int_{-\infty}^{\infty} N(v_{{\parallel}}, \tau) \, {\mathrm{d}} v_{{\parallel}} = 1. \end{equation}

Equation (5.19a,b) admits a similarity solution. To see this, we introduce a dilation transformation:

(5.21a–c)\begin{equation} z = \delta^{a}v_{{\parallel}}, \quad s = \delta^{b}\tau, \quad v(s,z) = \delta^{c}N (\delta^{{-}a}z,\delta^{{-}b}s). \end{equation}

Our problem then becomes

(5.22)\begin{equation} \delta^{b-a} \frac{\partial v}{\partial s} - \beta \delta^{2a} \frac{\partial ^{2}}{\partial z^{2}}\left(\frac{1}{2}\delta^{{-}2c}v^{2}\right) = 0. \end{equation}

From this, we see that (5.19a,b) is invariant under the dilation transformation provided $b -c = 2(a-c)$, i.e. $c = 2a - b$.

We therefore pose the self-similar ansatz

(5.23a,b)\begin{equation} N(v_{{\parallel}},\tau) = \tau^{(2a-b)/c}y(\xi), \quad \xi = \frac{v_{{\parallel}}}{\tau^{a/b}}. \end{equation}

The condition (5.20) then gives

(5.24)\begin{equation} \int_{-\infty}^{\infty} N(v_{{\parallel}},\tau) \, {\mathrm{d}} v_{{\parallel}} = 1 \implies \tau^{3(a/b)-1} \int_{-\infty}^{\infty} y(\xi) \, {\mathrm{d}} \xi = 1. \end{equation}

This condition must hold for all $\tau > 0$ and hence we have $3a=b$. So we can write

(5.25a,b)\begin{equation} N(v_{{\parallel}},\tau) = \tau^{{-}1/3}y(\xi), \quad \xi = \frac{v_{{\parallel}}}{\tau^{1/3}}. \end{equation}

It is of course no coincidence to have arrived at this particular self-similar ansatz. We could have also posited this solution simply by noting that the quantity $\xi = v_{\parallel }^{3}/\tau$ is dimensionless, which we could have deduced from (5.6).

The derivatives then transform as

(5.26)\begin{gather} \frac{\partial N}{\partial \tau} ={-}\frac{1}{3}\tau^{{-}4/3} (y + \xi y^{\prime}), \end{gather}

(5.27)\begin{gather}\frac{\partial ^{2}}{\partial v_{{\parallel}}^{2}}\left(\frac{1}{2}N^{2}\right) = \frac{1}{2}\tau^{{-}2/3} (2yy^{\prime\prime} + 2 y^{\prime2} ) \tau^{{-}2/3}, \end{gather}

where the prime notation now denotes differentiation with respect to $\xi$.

Equation (5.19a,b) becomes

(5.28)\begin{equation} \beta y y^{\prime\prime} + \beta y^{\prime 2} + \tfrac{1}{3}y + \tfrac{1}{3}\xi y^{\prime} = 0, \end{equation}

or equivalently

(5.29)\begin{equation} 3\beta (yy^{\prime})^{\prime} + (\xi y)^{\prime} = 0. \end{equation}

Direct integration gives

(5.30)\begin{equation} 3\beta yy^{\prime} + \xi y = k_{0}. \end{equation}

The integration constant $k_{0}$ is chosen such that the flux of particles $N ({{\mathrm {d}} N}/{{\mathrm {d}} v_{\parallel }})$ is continuous at the edge of the support where $N$ goes to zero. In our case, this is the same constant that one gets from imposing zero flux at the symmetry point because the flux has to vanish either when $N=0$ (at the edge) or when ${{\mathrm {d}} N}/{{\mathrm {d}} v_{\parallel }} = 0$ (at a symmetry point).

This equation can then be integrated once again:

(5.31)\begin{equation} y(\xi) ={-}\frac{1}{6\beta}\xi^{2} + k_{1}. \end{equation}

In order to satisfy the requirement that $N \rightarrow 0$ as $| v_{\parallel } | \rightarrow \infty ,$ we write the solution as

(5.32)\begin{equation} N(v_{{\parallel}},\tau) = \begin{cases} \dfrac{\tau^{{-}1/3}}{6\beta}\left(A^{2} - \left(\dfrac{v_{{\parallel}}}{\tau^{1/3}}\right)^{2}\right) & |v_{{\parallel}}| < A\tau^{1/3} \\ 0 & |v_{{\parallel}}| > A\tau^{1/3}, \end{cases} \end{equation}

where $A = \sqrt [3]{9\beta /2}$ so as to satisfy (5.20).

We can thus obtain a weak, in the sense that the first derivative is discontinuous, solution of (4.5a–c) in the limit where electron–ion collisions are neglected. It is imperative to mention that this solution holds only when $\tau < 0.$ In essence, what we have done in this section is to solve an ill-posed backwards diffusion equation (5.18) by solving instead a well-posed forward diffusion equation (5.19a,b) and running the solution backwards in time. The rationale in adopting this approach is to motivate the correct boundary conditions to apply to the equation.

As this problem is ill-posed, our self-similar solution requires very specific initial conditions. An arbitrarily small perturbation (in any given Sobolev norm) can cause the solution to fail and develop an infinitesimal-time singularity. Of course the self-similar solution is still of value, and the formation of singularities is somewhat fictitious since our assumption that $\boldsymbol {\nabla }_{{\parallel }} \ll \boldsymbol {\nabla }_{{\perp }}$ will certainly break down before any singularity arises. The solution obtained is of the correct general form (i.e. it attains a maximum for $v_{\parallel } = 0,$ and satisfies $N \rightarrow 0$ as $|v_{\parallel }|\rightarrow 0$) as to possibly be approached from more general initial conditions. Such solutions are not uncommon in plasma physics and the generation of sharp structures by resistive or ambipolar diffusion has also been investigated (Low Reference Low1973).

It is also fitting to remark here that including higher powers of the expansion parameter $\epsilon$ in our expansion of the collision operator (i.e. in the calculation of $I_{ee}$) would usually bring higher-order derivatives of $N$ into (5.18) which would regularise the partial differential equation. This is a possible avenue of further exploration with the caveat that the underlying equations would then likely be rendered too difficult to solve with analytical techniques.

As before, one can calculate the perpendicular temperature:

(5.33)\begin{equation} T_{{\perp}}(v_{{\parallel}},\tau) = \begin{cases} 2\sigma\tau_{r}m|\ln\epsilon|\dfrac{\tau^{{-}1/3}}{3\beta}\left(A^{2} - \left(\dfrac{v_{{\parallel}}}{\tau^{1/3}}\right)^{2}\right) & |v_{{\parallel}}| < A\tau^{1/3} \\ 0 & |v_{{\parallel}}| > A\tau^{1/3}. \end{cases} \end{equation}

5.2.1. Rate of energy loss

Following the previous subsection, one can calculate the rate at which energy is radiated from the plasma by first calculating $W_{\perp }$ and then appealing to equation (5.12a,b). This is given by

(5.34)\begin{equation} W_{{\perp}} = \int \frac{nT_{{\perp}}(v_{{\parallel}},t)^{2}}{4\sigma\tau_{r}m|\ln \epsilon|} \, {\mathrm{d}}v_{{\parallel}} = 4 mn \sigma \tau_{r} |\ln \epsilon| \left(\frac{3}{4\tau A^{3}}\right)^{2} \int_{-\infty}^{\infty} (A^{2} \tau^{2/3} - v_{{\parallel}}^{2})^{2}\, {\mathrm{d}}v_{{\parallel}}, \end{equation}

which can be evaluated exactly to give

(5.35)\begin{equation} W_{{\perp}} = 8 m n\sigma \tau_{r}|\ln \epsilon| \left(\frac{3}{4\tau A^{3}}\right)^{2} \int_{0}^{A\tau^{1/3}}(A^{2}\tau^{2/3} - v_{{\parallel}}^{2}) \, {\mathrm{d}}v_{{\parallel}} = \frac{12}{5} \frac{mn \sigma \tau_{r} |\ln \epsilon|}{A\tau^{1/3}}. \end{equation}

This leads to an energy loss rate of

(5.36)\begin{equation} -\frac{{\mathrm{d}} W}{{\mathrm{d}} t} = \frac{12}{5} \left(\frac{1}{2Z}\right)^{1/3} \frac{mn\sigma^{2/3}|\ln \epsilon|}{At^{1/3}}, \end{equation}

which can be written in terms of the Braginskii electron–electron collision time $\tau _{ee}$ as

(5.37)\begin{equation} -\frac{{\mathrm{d}} W}{{\mathrm{d}} t} ={-} \frac{6}{5}{\rm \pi}^{1/3} n_{e}T_{{\parallel} e} |\ln \epsilon|^{2/3} \left(\frac{1}{t\tau_{ee}^{2}}\right)^{1/3}. \end{equation}

At first glance, this equation might appear somewhat peculiar. We must first of course recall that our solution was valid only when time is run backwards, and hence the term $t$ appearing on the right-hand side is negative; as a result, the total energy does indeed decrease as expected.

5.3. Validity of the collisional approach; a self-accelerating process

We have found the solution to (4.5a–c) in two different limits based on the plasma parameters, thus elucidating the plasma distribution function in the collisional regime (3.19). Throughout this work, a number of orderings have been invoked and we thus turn our attention to the validity of our solution; specifically, when it is valid to treat $\epsilon$ as constant. Thus far we have focused on the evolution of the plasma on time scales $t$ with

(5.38)\begin{equation} 0 < \tau_{r} < t \sim \tau_{c}. \end{equation}

In this regime, we treat $\epsilon$ as constant seeing as though any variation is very small. However, as per our earlier remarks in § 2.3, we must take care to recall that $\tau _c$ is not constant and will itself decrease as the plasma cools via cyclotron emission. Our entire process (radiative cooling and collisional scattering) is therefore self-accelerating. Collisional scattering converts parallel kinetic energy to perpendicular kinetic energy. This perpendicular energy is then radiated, lowering the temperature of the plasma and in turn increasing the collisional scattering rate. As a result, the plasma will only remain in this regime for a few collision times, before entering a regime where

(5.39)\begin{equation} \tau_{c} < \tau_{r} < t. \end{equation}

In this final regime, any remaining parallel kinetic energy will be converted to perpendicular kinetic energy and then radiated. The total energy of the plasma will thus decay exponentially quickly on the radiation time scale.

In essence, we can partition the evolution of the plasma energy into three distinct regimes:

1. (i) Initially, $0 < \tau _{r} < t < \tau _{c}.$ In this regime the behaviour is collisionless as described in § 2. The perpendicular energy decays exponentially quickly on the radiation time scale.

2. (ii) After some time, $\tau _{r} < t \sim \tau _{c}$ and collisional scattering becomes important. This regime is what we have studied in §§ 3–5. This is a self-accelerating process and the plasma will only remain in this regime for a few collision times.

3. (iii) Eventually, $\tau _{c} < \tau _{r} < t.$ In this regime, the distribution function is isotropic (Maxwellian), and any remaining parallel kinetic energy will be converted to perpendicular kinetic energy via collisional scattering and then radiated.

We have now accomplished our second aim; namely, we have been able to calculate the evolution of the distribution function under collisional scattering. Unusually for plasma physics, this was possible even though the distribution function is far from Maxwellian. Thus far, we have made the restrictive assumption of straight-field-line geometry. We will now ask ourselves whether we can relax this assumption.

6.1. The bounce-averaged kinetic equation

In order to remove the dependence on $f_{1}$, we define the bounce average of a function $Q(\mu ,w,l)$ by

(6.8)\begin{equation} \bar{Q}(\mu,w) = \int Q(\mu,w,l) \left.\frac{{\mathrm{d}} l}{v_{{\parallel}}} \right/\int \frac{{\mathrm{d}} l}{v_{{\parallel}}}, \end{equation}

with

(6.9)\begin{equation} v_{{\parallel}} = \sqrt{\frac{2}{m}(w-\mu B)}. \end{equation}

The integral in (6.8) is taken between consecutive bounce points, defined by $v_{\parallel } = 0$, for trapped particles. For circulating particles, the integral is taken once around a field line if the field line is closed. If the field line is not closed, but instead traces out a magnetic surface, as in a stellarator or tokamak, then the bounce average for circulating particles is given by

(6.10)\begin{equation} \bar{Q}(\mu,w) =\lim_{L \rightarrow \infty} \int_{{-}L}^{L} Q(\mu,w,l) \left.\frac{{\mathrm{d}} l}{v_{{\parallel}}} \right/\int_{{-}L}^{L} \frac{{\mathrm{d}} l}{v_{{\parallel}}}. \end{equation}

Taking the bounce average of (6.7) then yields

(6.11)\begin{equation} \frac{\partial f_{0}}{\partial t} - \overline{\frac{1}{\tau_{r}}\frac{\partial }{\partial \mu}(\mu f_{0})} - \mu \overline{\left(\frac{B}{\tau_{r}}\right)}\frac{\partial f_{0}}{\partial w} = \overline{C_{ee}(\,f_{0})} + \overline{C_{ei}(\,f_{0})}. \end{equation}

Now, here we must be careful when bounce-averaging as $\tau _{r} \propto B^{-2}$. Let us write

(6.12)\begin{equation} \tau_{r} = \tau_{0}\left(\frac{B}{B_{0}}\right)^{{-}2}, \end{equation}

where $B_{0}$ and $\tau _{0}$ appearing on the right-hand side of this equation are constants. We thus arrive at the bounce-averaged kinetic equation

(6.13)\begin{equation} \frac{\partial f_{0}}{\partial t} - \frac{1}{\tau_{0}} \overline{\frac{B^{2}}{B_{0}^{2}}}\frac{\partial }{\partial \mu}\left(\mu f_{0}\right) - \frac{\mu}{\tau_{0}} \overline{\frac{B^{3}}{B_{0}^{2}}} \frac{\partial f_{0}}{\partial w} = \overline{C_{ee}(\,f_{0})} + \overline{C_{ei}(\,f_{0})}. \end{equation}

Since $\overline {B^{2}}$ and $\overline {B^{3}}$ are, in general, complicated functions of $\mu /w,$ one cannot hope to solve this equation in full generality. However, we can deduce the distribution function on a significantly long time scale $\tau$ which is larger than the radiation time and comparable to the collision time $\tau _{r} \ll \tau \lesssim \tau _{c}$. In this instance, the perpendicular energy will have been radiated away to leading order and we have

(6.14)\begin{equation} \lambda B \ll 1, \end{equation}

where we have introduced $\lambda = \mu /w$. Trapped particles will thus be absent from the population. In this regime, one can neglect the third term in (6.13) and thus obtain

(6.15)\begin{equation} \frac{\partial f_{0}}{\partial t} - \frac{1}{\tau_{0}}\overline{\frac{B^{2}}{B_{0}^{2}}}\frac{\partial }{\partial \mu}\left(\mu f_{0}\right) = \overline{C_{ee}(\,f_{0})} + \overline{C_{ei}(\,f_{0})}. \end{equation}

Moreover, we can also evaluate

(6.16)\begin{equation} \overline{\frac{B^{2}}{B_{0}^{2}}} = \frac{1}{B_{0}^{2}}\oint \left.\frac{B^{2}(l)}{\sqrt{1-\lambda B}} \,{\mathrm{d}} l \right/ \oint\frac{{\mathrm{d}} l}{\sqrt{1-\lambda B}} = \left.\frac{1}{B_{0}^{2}}\oint B^{2} \,{\mathrm{d}} l \right/\oint {\mathrm{d}} l \equiv 1, \end{equation}

where the last equivalence follows from the choice to simply define $B_{0}^{2}$ to be the average of $B^{2}$ over a field line.

Hence, we arrive at

(6.17)\begin{equation} \frac{\partial f_{0}}{\partial t} - \frac{1}{\tau_{0}}\frac{\partial }{\partial \mu}(\mu f_{0}) = \overline{C_{ee}(\,f_{0})} + \overline{C_{ei}(\,f_{0})}. \end{equation}

6.2. The bounce-averaged collision operator

We have already shown that at leading order, the contribution from electron–ion collisions is formally smaller than that from electron–electron collisions provided $n_{i}Z^{2}/n_{e} = O(1).$ Let us thus restrict our attention to the electron–electron collision operator.

We know that at leading order

(6.18)\begin{equation} C_{ee}(\,f_{0}) \simeq \frac{\sigma g(v_{{\parallel}})}{v_{{\perp}}}\frac{\partial }{\partial v_{{\perp}}}v_{{\perp}} \frac{\partial f_{0}}{\partial v_{{\perp}}}. \end{equation}

We must be cautious as $v_{\perp }$ varies along the orbit and should not be used as a coordinate in the bounce-averaged equation. When the collision operator is instead written using the magnetic moment $\mu$ as a coordinate, we obtain

(6.19)\begin{equation} \overline{C_{ee}(\,f_{0})} \simeq 2\sigma m\overline{\frac{g(v_{{\parallel}})}{B}} \frac{\partial }{\partial \mu}\left(\mu\frac{\partial f_{0}}{\partial \mu}\right) \end{equation}

and thus

(6.20)\begin{equation} \frac{\partial f_{0}}{\partial t} - \frac{1}{\tau_{0}}\frac{\partial }{\partial \mu}(\mu f_{0}) = 2 \sigma m \overline{\frac{g(v_{{\parallel}})}{B}} \frac{\partial}{\partial \mu}\left(\mu \frac{\partial f_{0}}{\partial \mu}\right). \end{equation}

It is of course no accident that this equation looks remarkably similar to the equation obtained in the straight-field-line limit. One can make the analogy exact and state that the leading-order dynamics of both systems are governed by the differential equation provided we associate the $B^{2}$ appearing in Larmor's formula, and the function $g(v_{\parallel })$ arising from the electron–electron collision operator, with their averages over a magnetic field line.

That is, on time scales $\tau$ with

(6.21)\begin{equation} \tau_{b} \ll \tau_{r} \ll \tau \lesssim \tau_{c} \end{equation}

the distribution function will become strongly anisotropic in general magnetic geometry.

7.1. Cyclotron sources in the laboratory

Non-neutral plasmas, specifically plasmas consisting of charged particles with a single sign of charge can be confined in Penning–Malmberg traps (Dubin & O'Neil Reference Dubin and O'Neil1999; Danielson et al. Reference Danielson, Dubin, Greaves and Surko2015). Such a trap consists of a vacuum region inside an electrode structure consisting of a stack of hollow, metal cylinders. A uniform axial magnetic field is then applied to inhibit particle motion in the radial direction. Voltages must also be imposed on the end electrodes to prevent particle loss in the magnetic field direction.

It is frequently useful to compress plasmas radially; for instance, to increase the plasma density. This is usually accomplished by applying a torque on the plasma using rotating electric fields, the so-called ‘rotating wall technique’ (Anderegg, Hollmann & Driscoll Reference Anderegg, Hollmann and Driscoll1998). Very long confinement times (of the order of hour or days) can be achieved using these techniques, making their use highly desirable. Particle cooling is often necessary to maintain good confinement by mitigating the heating caused by the torque using the rotating wall method. In the case of electrons or positrons, if the magnetic field is sufficiently strong, the particles will cool by cyclotron radiation (O'Neil Reference O'Neil1985). An example of where cyclotron cooling is employed successfully to this end is in the production of anti-hydrogen, where this process is used cool pure electron plasma to sub-eV temperatures (Amoretti et al. Reference Amoretti, Carraro, Lagomarsino, Manuzio, Testera and Variola2002).

Another major application of the theory developed in this work will be for the first laboratory experiment, currently under development, to create and confine the first terrestrial electron–positron plasmas in the laboratory. This is done by first accumulating positrons from a powerful source and then injecting these into a pure electron plasma confined by the dipolar magnetic field of a current-carrying circular coil, so that a stationary, quasi-neutral electron–positron plasma is formed (Sunn Pedersen et al. Reference Sunn Pedersen, Danielson, Hugenschmidt, Marx, Sarasola, Schauer, Schweikhard, Surko and Winkler2012). This system will satisfy the necessary conditions of being optically thin and with a radiation time which should be initially shorter than the collision time. However, there is an additional complication in that the plasma will be so strongly magnetised that Coulomb collisions are no longer effectively described by the Landau operator. The theory of scattering in strongly magnetised anisotropic pair plasmas is developed in the companion paper (II).

7.2. Synchrotron sources in astrophysics

This work was built on the fact that a charged particle moving in a magnetic field radiates energy. At non-relativistic energies, the focus of this paper, this process is called cyclotron cooling. At relativistic velocities it is known as synchrotron radiation.

Synchrotron sources are ubiquitous and the emission of relativistic and ultra-relativistic electrons gyrating in a magnetic field is a process that dominates much of high-energy astrophysics. Indeed, it is known that synchrotron radiation is responsible for the non-thermal optical and X-ray emission observed in the Crab Nebula (Rees & Gunn Reference Rees and Gunn1974) and that pulsars are strong synchrotron sources (Sturrock Reference Sturrock1971).

Of course, the physics of these systems should also be treated carefully using a model taking relativistic effects into account in the modelling of the radiation term and the collision terms. This is an area of active research and a potential axis along which this current work could be further developed. There also likely exist analogous systems in astrophysics where the electrons are subrelativistic but nevertheless are still strongly radiating and weakly collisional in the sense discussed here. Such a system would of course still be required to satisfy

(7.1)\begin{equation} \frac{\tau_{r}}{\tau_{c}} = \left(\frac{\omega_{p}}{\varOmega_{c}}\right)^{2}\left(\frac{c}{v_{\text{th}}}\right)^{3} \ll 1. \end{equation}

The easiest way to satisfy this is of course to allow relativistic plasmas, but we might also envision some exotic astrophysical plasma which is sufficiently nebulous (very low density) but also sufficiently strongly magnetised (very strong magnetic fields) as to allow this condition to be met in a non-relativistic plasma. In such systems, the theory in this paper will be directly relevant.

However, these is another caveat here which must be carefully considered. We have presented here a mechanism through which plasmas can become strongly anisotropic in velocity space. In actuality, there is a plethora of instabilities which could act on astrophysical plasmas and restore isotropy on time scales shorter than the collision time. It seems to be that the question of (in)stability could depend sensitively on plasma parameters and that considerations of isotropy-restoring instabilities do not necessarily preclude the types of systems studied in this work.