2.1 Drift kinetic equation

Let a long, straight, current carrying wire be coincident with the z-axis of a cylindrical coordinate system ($r,\theta, z$). The magnetic field in the vicinity of the wire is then given by the expression

(2.1)\begin{align} \boldsymbol{B} = \frac{2I}{cr} \hat{\boldsymbol{\theta}}, \end{align}

and the vector potential by the expression

(2.2)\begin{equation} \boldsymbol{A} = A_z(r)\hat{\boldsymbol{z}} ={-}\frac{2I}{c} \ln (r) \hat{\boldsymbol{z}}. \end{equation}

The confinement region is bounded by the conducting surface of the wire at radius $r_{\rm wire}$ and an outer conducting cylinder at radius $r_{{\rm wall}}$, where $r_{{\rm wire}} < r_{{\rm wall}}$. A schematic of this configuration is given by figure 1.

Figure 1. Schematic of the configuration with the wire in red, the plasma in blue and the wall in grey. The current through the wire is indicated by the black arrow and the magnetic field by the blue arrow.

Following Taylor (Reference Taylor1964), the guiding-centre-drift Lagrangian of a particle moving in the field of the wire can be written in the form

(2.3)\begin{equation} L_j = \frac{m_jr^2 \dot{\theta}^2}{2} + \frac{q_j}{c} A_z \dot{z} - q_j\phi - \mu_j |B|, \end{equation}

where $j = {\rm e}$ refers to electrons and $j = {\rm p}$ refers to positrons. The respective particle has charge $q_j$, mass $m_j$ and magnetic moment $\mu _j = m_jv_{\perp }^2/2|B|$. Since the mass for electrons and positrons is the same we will drop that index from now on. The first term in the Lagrangian is the kinetic energy associated with the velocity component parallel to the magnetic field line, the second is the vector potential term, the third is the electrostatic potential energy and the last term is the potential energy associated with the magnetic moment. Both potential energy terms enter the Lagrangian with a minus sign. The canonical momenta for the species $j$ are given by the expressions

(2.4a,b)\begin{align} p_{\theta} = \frac{\partial L_j}{\partial \dot{\theta}} = mr^2 \dot{\theta} \quad \textrm{and} \quad p_{zj} = \frac{\partial L_j}{\partial \dot{z_j}} = \frac{q_j}{c} A_z. \end{align}

The corresponding guiding-centre-drift Hamiltonian is given by the expression

(2.5)\begin{equation} H_j = p_{\theta} \dot{\theta} + p_{zj} \dot{z} - L_j = \frac{p_{\theta}^2}{2mr^2} + \mu |B| + q_j \phi. \end{equation}

The canonically conjugate variables are ($\theta, p_{\theta }$) and ($z,p_{zj}$), where $p_{zj}$ represents any dependence on the radius $r$. Neglecting collisions for the relatively short time scales relevant to mode growth, the adiabatic invariant $\mu$ can be treated as a constant of the motion.

The distribution of guiding centres is given by the expression $f_j=f_j(\theta, p_{\theta },z,p_{zj},\mu,t)$ and satisfies the kinetic equation

(2.6)\begin{equation} \frac{\partial f_j}{\partial t} + [f_j, H_j] = 0. \end{equation}

The number of guiding centres within a phase space volume is given by the expression

(2.7)\begin{equation} {\rm d}N_j = f_j(\theta, p_{\theta},z,p_{zj},\mu,t)|\,{\rm d}\theta \,{\rm d}p_{\theta} \,{\rm d}z\,{\rm d}p_{zj} \,{\rm d}\mu| \end{equation}

and the density of guiding centres in physical space is given by the expression

(2.8)\begin{equation} n_j = \int \frac{{\rm d}N_j}{|{\rm d}z r \,{\rm d}\theta \,{\rm d}r|} = \left| \frac{2Iq_j}{c^2 r^2} \right| \int_0^\infty {\rm d}\mu \int_{-\infty}^\infty {\rm d}p_{\theta}\,f_j, \end{equation}

where use has been made of the relation $|{\rm d}p_{zj}| = |{2Iq_j}/{c^2 r}|\,{\rm d}r$. The drift in the z-direction follows from the Hamiltonian in (2.5)

(2.9)\begin{align} \dot{z_j} = \frac{\partial H_j}{\partial p_{zj}} ={-} \frac{p_{\theta}^2}{mr^3}\frac{\partial r}{\partial p_{zj}} + \mu \frac{\partial |B|}{\partial r} \frac{\partial r}{\partial p_{zj}} + q_j\frac{\partial \phi}{\partial r} \frac{\partial r}{\partial p_{zj}} = \frac{c^2}{2q_jI}\left( mv_{{\parallel}}^2 + \frac{mv_{{\perp}}^2}{2} \right) - \frac{c^2r}{2I}\frac{\partial \phi}{\partial r}, \end{align}

where ${\partial r}/{\partial p_{zj}} = -{c^2r}/{2q_jI}$. The first two terms on the right-hand side of (2.9) are the curvature and gradient-B drifts and the last term is the $E\times B$ drift. For a given electric potential, the $E\times B$ drift is independent of the sign of charge in contrast to the curvature and gradient-B drift.

2.2 Linear mode equation

In linear theory the distribution for species $j$ can be written in the form

(2.10)\begin{equation} f_j = f_j^{(0)}(p_{\theta}, p_{zj}, \mu) + f_j^{(1)}(p_{\theta}, p_{zj}, \mu) \exp[{\rm i}(kz-\omega t)], \end{equation}

where the first term is the equilibrium distribution and the second term is the perturbation. The corresponding electric potential is written in the form

(2.11)\begin{equation} \phi = \phi^{(0)}(p_{zj}) + \phi^{(1)}(p_{zj}) \exp[{\rm i}(kz-\omega t)]. \end{equation}

We neglect the finite resistivity of the bounding conductors. Consequently, the mode potential vanishes at both boundaries. The $p_{zj}$ in the argument of the potential functions represents radial dependence and must be used in the corresponding Hamiltonian $H_j$. Assume an equilibrium distribution of the form

(2.12)\begin{equation} f_j^{(0)} = g(p_{zj}) \exp\left[ -\frac{1}{T}(H_j^{(0)} - q_j\phi^{(0)})\right], \end{equation}

with the equilibrium Hamiltonian $H_j^{(0)} = {p_{\theta }^2}/{2mr^2} + \mu |B| + q_j \phi ^{(0)}$. The second factor in (2.12) is a Maxwellian velocity distribution. We consider the case where the temperature is independent of the radius and is the same for both species. Carrying out the integrals in (2.8) shows that $g_j(p_{zj})$ is the equilibrium density up to a multiplicative constant.

(2.13)\begin{equation} n_j^{(0)}(r) = \frac{q_jT^{3/2}}{c} \sqrt{\frac{{\rm \pi} m}{2}} g_j(p_{zj}). \end{equation}

Using the functional dependence of $f_j^{(0)}(p_{\theta }, p_{zj}, \mu )$ and $H_j^{(0)}(p_{\theta }, p_{zj}, \mu )$ and linearising the kinetic equation in the perturbation yields the relation

(2.14)\begin{equation} \left(\omega - k \frac{\partial H_j^{(0)}}{\partial p_{zj}} \right) f_j^{(1)} ={-} kq_j\phi^{(1)} \frac{\partial f_j^{(0)}}{\partial p_{zj}}, \end{equation}

with the solution

(2.15)\begin{equation} f_j^{(1)} ={-} \frac{kq_j\phi^{(1)} \dfrac{\partial f_j^{(0)}}{\partial p_{zj}}}{\omega - k\dfrac{\partial H_j^{(0)}}{\partial p_{zj}}}. \end{equation}

Poisson's equation for the mode potential then takes the form

(2.16)\begin{align} \frac{1}{r} \frac{\partial}{\partial r} r \frac{\partial \phi^{(1)}}{\partial r} -k^2 \phi^{(1)} & ={-}4{\rm \pi} \sum_j q_j n_j^{(1)} \nonumber\\ & = \sum_j 8{\rm \pi} q_j^2 \left|\frac{ q_j I}{c^2 r^2}\right| k\phi^{(1)} \int_0^\infty {\rm d}\mu \int_0^\infty {\rm d}p_{\theta} \frac{\dfrac{\partial f_j^{(0)}}{\partial p_{zj}}}{\omega - k\dfrac{\partial H_j^{(0)}}{\partial p_{zj}}}, \end{align}

where use has been made of (2.8).

To solve (2.14) we divided by the resonance factor $\omega - k ({\partial H^{(0)}}/{\partial p_z})$, which can be zero for some value of $p_z$ if ${\rm Im}(\omega )=0$. This is the same problem as faced by Landau when solving the Vlasov equation (Landau Reference Landau1946). The underlying difficulty is that (2.15) is a proper initial value solution to (2.14) only for a growing wave where ${\rm Im}(\omega )>0$ is positive. As ${\rm Im}(\omega )$ approaches zero the contour of integration must be deformed to stay on the same side of any pole or branch cut. We avoid this complication by solving the mode equation only for growing modes. This is sufficient to identify domains of instability. However, it is important to keep in mind that a growing mode is not a general solution for an arbitrary initial density perturbation. The eigenfunction of such a mode does not necessarily match the radial dependence of an arbitrary initial perturbation. However, if this mode is the fastest-growing mode and has a non-zero overlap integral with the initial perturbation, it is the time-asymptotic limit of the solution for the initial perturbation. Likewise, the plasma modes in Landau's solution dominate only in the time-asymptotic limit.

The integral on the right-hand side can be rewritten in the form

(2.17)\begin{align} & \int_0^\infty {\rm d}\mu \int_0^\infty {\rm d}p_{\theta} \frac{\partial f_j^{(0)} / \partial p_{zj}}{\omega - k\dfrac{\partial H_j^{(0)}}{\partial p_{zj}}} \nonumber\\ & \quad = \frac{\sqrt{m}T^{{3}/{2}}r^2c}{4I} \int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]}{\omega - kq_j\dfrac{\partial \phi}{\partial p_{zj}} - \dfrac{kc^2T}{2Iq_j}\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} \nonumber\\ & \qquad\times \left[ \frac{\partial g_j(p_z)}{\partial p_{zj}} - \frac{c^2 g(p_{zj})}{2q_jI} \left(\frac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right) \right], \end{align}

with the scaled velocities $x_{\parallel } = \sqrt {{m}/{T}} \,v_{\parallel }$ and $x_{\perp } = \sqrt {{m}/{T}} \, v_{\perp }$. Here, use had been made of ${\rm d}\mu \,{\rm d}p_{\theta } = ({\sqrt {m}T^{{3}/{2}}r}/{2B})\,{\rm d}\kern0.7pt x_{\perp }^2\,{\rm d}\kern0.7pt x_{\parallel }$. Gauss's law and the relation ${\partial r}/{\partial p_{zj}} = -{c^2r}/{2q_jI}$ imply the result

(2.18)\begin{align} kq_j\frac{\partial \phi}{\partial p_{zj}} = \frac{2{\rm \pi} q_{\rm e} c^2 k (1-\eta)}{I} \int_0^r {\rm d}r'\,r' n_{\rm e}^{(0)}, \end{align}

where $\eta = {n_{\rm p}^{(0)}(r)}/{n_{\rm e}^{(0)}(r)}$ is the ratio of positron to electron density. We have assumed that the two species are uniformly mixed. Expressing $g_j(p_{zj})$ in terms of the number density $n_j$ according to (2.13) and plugging (2.17) back into Poisson's equation (2.16) yields the equation

(2.19)\begin{align} & \frac{1}{r} \frac{\partial}{\partial r} r \frac{\partial \phi^{(1)}}{\partial r} -k^2 \phi^{(1)}\nonumber\\ & \quad ={-}\sum_j \frac{\sqrt{2 {\rm \pi}}q_jkc^2\phi^{(1)}}{I} \nonumber\\ & \qquad \times\int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right] \left[ r\dfrac{\partial n_j^{(0)}}{\partial r} + n_j^{(0)}\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right) \right]}{\omega + \dfrac{2{\rm \pi} q_{\rm e} c^2 k (1-\eta)}{I} \displaystyle\int\nolimits_0^r {\rm d}r'\,r' n_{\rm e}^{(0)} - \dfrac{kc^2T}{2Iq_j}\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)}. \end{align}

We now want to bring (2.19) into a dimensionless form. To do so, we introduce the scaled radius $\xi = kr$ and the scaled density $\tilde {n}=n/n_0$ where $n_0$ is a constant characteristic value of the density. As an example, consider the Gaussian density profile

(2.20)\begin{equation} n_{\rm e}^{(0)}(r) = n_{0{\rm e}} \exp\left[ -\frac{(r-r_0)^2}{2\Delta r^2}\right], \end{equation}

which when scaled takes the form

(2.21)\begin{equation} \tilde{n}(\xi/k) = \exp\left[ -\frac{(\xi-\xi_0)^2}{2\Delta \xi^2}\right], \end{equation}

where $\xi _0 = kr_0$ is the scaled peak and $\Delta \xi = \Delta r k$ is the scaled width. For notational simplicity, we denote $\tilde {n}(r) = \tilde {n}(\xi /k)$ by $\tilde {n}(\xi )$ remembering that all lengths are scaled by $k$. The mode equation (2.19) can be rewritten in the form

(2.22)\begin{align} & \frac{1}{\xi} \frac{\partial}{\partial \xi} \xi \frac{\partial \phi^{(1)}}{\partial \xi} -\phi^{(1)} \nonumber\\ & \quad ={-}\frac{\phi^{(1)}}{\sqrt{2 {\rm \pi}}} \int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]\left[ \xi\dfrac{\partial \tilde{n}^{(0)}_{\rm e}}{\partial \xi} + \tilde{n}^{(0)}_{\rm e}\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right) \right]}{\varOmega - (1-\eta)\displaystyle\int\nolimits_0^{\xi} {\rm d}\xi'\, \xi'\tilde{n}^{(0)}_{\rm e} - k^2\lambda_{\rm D}^2\left(\frac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} \nonumber\\ & \qquad +\frac{\eta \phi^{(1)}}{\sqrt{2 {\rm \pi}}} \int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]\left[ \xi\dfrac{\partial \tilde{n}^{(0)}_{\rm e}}{\partial \xi} + \tilde{n}^{(0)}_{\rm e}\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right) \right]}{\varOmega - (1-\eta)\displaystyle\int\nolimits_0^{\xi} {\rm d}\xi' \,\xi'\tilde{n}^{(0)}_{\rm e} + k^2\lambda_{\rm D}^2\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)}. \end{align}

The scaled frequency is given by $\varOmega = {\omega }/{\alpha } = -{I k \omega }/{2{\rm \pi} q_{\rm e} c^2 n_{0{\rm e}}}$. The Debye length $\lambda _{\rm D} = \sqrt {T/4{\rm \pi} q_{\rm e}^2 n_{0{\rm e}}}$ depends on the characteristic density $n_{0{\rm e}}$ that was introduced before. The importance of finite-temperature effects depends on the ratio $a = {k^2\lambda _{\rm D}^2}/{\varOmega - (1-\eta )\int _0^{\xi } {\rm d}\xi ' \,\xi '\tilde {n}^{(0)}}$ in regions where the plasma density is non-zero. We can simplify (2.22) further

(2.23)\begin{align} & \frac{1}{\xi} \frac{\partial}{\partial \xi} \xi \frac{\partial \phi^{(1)}}{\partial \xi} -\phi^{(1)} \nonumber\\ & \quad =\frac{\phi^{(1)}}{k^2\lambda_{\rm D}^2} \left[- \tilde{n}^{(0)} + \left(a\xi\frac{\partial \tilde{n}^{(0)}}{\partial \xi} + \tilde{n}^{(0)} \right)\frac{1}{\sqrt{2 {\rm \pi}}}\int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]}{1-a\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} \right] \nonumber\\ & \qquad -\frac{\eta \phi^{(1)}}{k^2\lambda_{\rm D}^2} \left[ \tilde{n}^{(0)} + \left(a\xi\frac{\partial \tilde{n}^{(0)}}{\partial \xi} - \tilde{n}^{(0)} \right)\frac{1}{\sqrt{2 {\rm \pi}}}\int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]}{1+a\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} \right]. \end{align}

The remaining integral can be carried out analytically by following Zocco et al. (Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018)

(2.24)\begin{equation} \frac{1}{\sqrt{2 {\rm \pi}}}\int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]}{1+a\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} =\frac{{\rm \pi} \exp\left[\dfrac{1}{a}\right] {\rm erfc}^2\left[\dfrac{1}{\sqrt{2a}}\right]}{2a}, \end{equation}

where erfc is the complementary error function for a complex argument. Note that the parameter $a$ does not depend on the scaled velocity variables.

3.1 Reduced description

This equation can be obtained directly by noting that in the zero-temperature approximation the two kinetic energy terms in the Hamiltonian (2.5) may be neglected and we are left with $H_{\rm e} = q_{\rm e}\phi (z, p_z, t)$. Because the velocity variables are not present in this reduced Hamiltonian, a reduced electron distribution can be used

(3.3)\begin{equation} h_{\rm e}(z, p_z, t) = \int_0^\infty {\rm d}\mu \int_0^\infty {\rm d}p_{\theta} \,f_{\rm e}(p_{\theta}, z, p_{z{\rm e}}, \mu, t). \end{equation}

The reduced distribution satisfies the reduced kinetic equation $\partial h_{\rm e} / \partial t + [h_{\rm e}, q_{\rm e}\phi ] = 0$. Linearising the reduced kinetic equation and combining it with Poisson's equation yields (3.2).

A necessary condition for instability can be established by following the same procedure as that used for a long, cylindrical, non-neutral plasma in a homogeneous magnetic field (Briggs et al. Reference Briggs, Daugherty and Levy1970). Let the scaled frequency be written as a real and an imaginary part, $\varOmega = \varOmega _r + {\rm i}\varGamma$, multiply mode equation (3.2) by the mode complex conjugate, integrate over ${\rm d}\xi \xi$ from the inner conductor to the outer conductor, use integration by parts on the partial derivative term and take the imaginary part of the resulting equation to obtain the requirement

(3.4)\begin{equation} \int_{\xi_{\rm wire}}^{\xi_{\rm wall}} {\rm d}\xi \frac{\varGamma |\phi^{(1)}|^2}{\left(\varOmega_r - \displaystyle\int\nolimits_0^{\xi} {\rm d}\xi' \,\xi'\widetilde{n_{\rm e}}^{(0)}\right)^2 +\varGamma^2} \frac{\partial (\widetilde{n_{\rm e}}^{(0)}\xi^2)}{\partial \xi} = 0. \end{equation}

Equation (3.4) can be satisfied for positive $\varGamma$ only when $\tilde {n}_{\rm e}^{(0)}\xi ^2$ is non-monotonic in $\xi$. Since the density is zero at the inner cylinder and at the outer cylinder but finite in between this requirement is always satisfied. The inner conductor is the surface of the cylindrical wire and the outer surface is a cylindrical boundary wall. Thus, it is not surprising that we find diocotron instabilities for several sample density profiles. The corresponding condition for a plasma in a Penning–Malmberg trap is that $\partial n_{\rm e}^{(0)}(r)/\partial r$ be non-monotonic. The $r^2$ difference between the two criteria is due to the $E\times B$ flow being incompressible in the canonical phase space of the reduced distribution, not in configuration space. This is the same factor that enters the ratio between the phase space density and the physical density $r^2n_{\rm e} = |2Ie/c^2|h_{\rm e}$, where $|2Ie/c^2|$ is a constant. For the uniform axial magnetic field of a Penning–Malmberg trap, the two spaces are the same. Even in the limit where the curvature and gradient-B drifts are negligible, the difference in the magnetic field is important.

Since the diocotron instability involves the requirement that the plasma cylinder is hollow, one might think that the plasma can be stabilised by charging the wire. However, this does not work for the case of a long straight wire, since both the electric field from the charged wire and the magnetic field fall off as $1/r$. Hence, the $E\times B$ drift from the charged wire is independent of $r$ and can be removed by a transformation to the drifting frame. In the lab frame, the drift produces a Doppler shift in the real part of the mode frequency, but does not change the imaginary part, that is, the growth rate.

3.2 Van Kampen method

An alternative to mode equation (3.2) for the electric potential is the Van Kampen eigenvalue equation for the density perturbation (Van Kampen Reference Van Kampen1955; Case Reference Case1959; Schecter et al. Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000). To this end, first note that the scaled Poisson equation for the mode potential takes the form

(3.5)\begin{align} \left(\frac{1}{\xi} \frac{\partial}{\partial \xi} \xi \frac{\partial}{\partial \xi}-1\right) \phi^{(1)} = \frac{1}{k^2} \nabla^2 \phi^{(1)} ={-}\frac{4{\rm \pi} q_{\rm e}}{k^2} n_{\rm e}^{(1)}. \end{align}

One can identify the left-hand side in (3.2) as the scaled, perturbed density $-4{\rm \pi} q n^{(1)}/k^2$. Hence, we can write (3.2) as

(3.6)\begin{align} \left(\varOmega - \int_0^{\xi} {\rm d}\xi' \,\xi'\tilde{n}_{\rm e}^{(0)}\right)\frac{4{\rm \pi} q_{\rm e}}{k^2} n_{\rm e}^{(1)} = 2\phi^{(1)}\frac{\partial (\tilde{n}_{\rm e}^{(0)}\xi^2)}{\partial \xi^2}. \end{align}

In Appendix A, we obtain Green's function

(3.7)\begin{align} & G(\xi, \xi') \nonumber\\ & \quad = \begin{cases} \dfrac{({I}_0(\xi_{\rm wall}){K}_0(\xi')-{I}_0(\xi'){K}_0(\xi_{\rm wall}))({I}_0(\xi_{\rm wire}){K}_0(\xi)-{I}_0(\xi){K}_0(\xi_{\rm wire}))}{{I}_0(\xi_{\rm wall}){K}_0(\xi_{\rm wire})-{I}_0(\xi_{\rm wire}){K}_0(\xi_{\rm wall})} & {\rm for} \ \xi_{\rm wire} < \xi < \xi' , \\[10pt] \dfrac{({I}_0(\xi_{\rm wall}){K}_0(\xi)-{I}_0(\xi){K}_0(\xi_{\rm wall}))({I}_0(\xi_{\rm wire}){K}_0(\xi')-{I}_0(\xi'){K}_0(\xi_{\rm wire}))}{{I}_0(\xi_{\rm wall}){K}_0(\xi_{\rm wire})-{I}_0(\xi_{\rm wire}){K}_0(\xi_{\rm wall})} & {\rm for} \ \xi' < \xi < \xi_{\rm wall}, \end{cases} \end{align}

where ${I}_0(\xi )$ and ${K}_0(\xi )$ are Bessel functions of an imaginary argument. The perturbed potential can be written in terms of the perturbed density by using Green's function

(3.8)\begin{equation} \phi^{(1)} ={-}4{\rm \pi} q_{\rm e}\int_{\xi_{\rm wire}}^{\xi_{\rm wall}} {\rm d}\xi' \,\xi' G(\xi, \xi')n_{\rm e}^{(1)}(\xi'), \end{equation}

yielding the Van Kampen eigenvalue equation for the density perturbation

(3.9)\begin{equation} \left(\varOmega - \int_0^{\xi} {\rm d}\xi' \,\xi'\tilde{n}_{\rm e}^{(0)}\right)n_{\rm e}^{(1)} + 2\frac{\partial (\tilde{n}_{\rm e}^{(0)}\xi^2)}{\partial \xi^2} \int_{\xi_{\rm wire}}^{\xi_{\rm wall}} {\rm d}\xi' \,\xi' G(\xi, \xi')n_{\rm e}^{(1)}(\xi')=0. \end{equation}

The derivation of (3.9) does not involve division by $\varOmega - \int _0^{\xi } {\rm d}\xi ' \,\xi '\tilde {n}^{(0)}$. There is no restriction that $\varGamma$ be greater than zero. The Van Kampen method systematically finds all modes, stable and unstable, without the need for analytic continuation. As we show, the eigenmodes include discrete modes as well as continuum modes (Case Reference Case1959). As mentioned earlier, our focus here is on unstable modes.

For completeness, we note that the different eigenmodes can be combined to provide a solution to an arbitrary initial value problem. Let $n_1^{(1)} (\xi )$ and $n_2^{(1)} (\xi )$ be two eigenfunctions corresponding to frequencies $\varOmega _1$ and $\varOmega _2$ and define the inner product

(3.10)\begin{equation} \langle n_1^{(1)} (\xi) | n_2^{(1)} (\xi) \rangle = \int_{\xi_{\rm wire}}^{\xi_{\rm wall}} {\rm d}\xi \,\xi \frac{n_1^{(1)} (\xi) n_2^{(1)} (\xi)^{*}}{\dfrac{\partial (\tilde{n}^{(0)}\xi^2)}{\partial \xi^2}}. \end{equation}

From (3.9) one finds the orthogonality relation

(3.11)\begin{equation} (\varOmega_1 -\varOmega_2^{*})\langle n_1^{(1)} (\xi) | n_2^{(1)} (\xi) \rangle = 0. \end{equation}

This relation is useful in expanding an arbitrary initial density perturbation in eigenmodes

(3.12)\begin{equation} n^{(1)} (\xi, t) = \sum_d a_d n_d^{(1)}(\xi) \exp[-{\rm i}\varOmega_d t] + \int {\rm d}\varOmega \,a(\varOmega) n_{\varOmega}^{(1)}\exp[-{\rm i}\varOmega_d t], \end{equation}

where $a_d$ and $a(\varOmega )$ are time-independent constants given in terms of the initial perturbation $n^{(1)} (\xi, t)$ by the relations

(3.13)\begin{gather} a_d = \frac{\langle n_1^{(1)} (\xi, t=0) | n_d^{(1)} (\xi) \rangle}{\langle n_d^{(1)} (\xi) | n_d^{(1)} (\xi) \rangle} \end{gather}

(3.14)\begin{gather}a(\varOmega) = \frac{\langle n_1^{(1)} (\xi, t=0) | n_{\varOmega}^{(1)} (\xi) \rangle}{\langle n_{\varOmega}^{(1)} (\xi) | n_{\varOmega}^{(1)} (\xi) \rangle}. \end{gather}

the sum over $a_d$ covers the discrete modes and the integral over $a(\varOmega )$ covers the continuum modes.

3.3 Analytic solution

Analytic solutions to mode equation (3.2) and Van Kampen equation (3.9) are possible for the special case where the equilibrium density profile is given by the expression

(3.15)\begin{equation} \tilde{n}_{\rm e}^{(0)} = \frac{\xi_{\rm in}^2}{\xi^2} \left[\varTheta (\xi^2-\xi_{\rm in}^2)-\varTheta (\xi^2-\xi_{\rm out}^2) \right], \end{equation}

where $\xi _{\rm in}$ and $\xi _{\rm out}$ are the inner and outer edge of the plasma density, which is generally not the same as the inner and outer boundary $\xi _{\rm wire}$ and $\xi _{\rm wall}$ and $\varTheta$ is a step function. This density profile corresponds to a flat-top profile for the equilibrium distribution function $h_{\rm e}^{(0)}(\xi )$. For this density distribution the radial derivative in (3.9)

(3.16)\begin{equation} \frac{\partial (\tilde{n}^{(0)}\xi^2)}{\partial \xi^2} = \xi_{\rm in}^2 \left[\delta (\xi^2-\xi_{\rm in}^2)-\delta (\xi^2-\xi_{\rm out}^2) \right], \end{equation}

is the sum of two delta functions. This feature simplifies both the mode equation and the Van Kampen equation, allowing analytic solutions. We will solve the Van Kampen equation since Green's function has already solved part of the problem. Equation (3.9) implies that the density perturbation is a sum of two delta functions as well:

(3.17)\begin{equation} n^{(1)} = C_{\rm in}\delta (\xi^2-\xi_{\rm in}^2) - C_{\rm out}\delta (\xi^2-\xi_{\rm out}^2). \end{equation}

Evaluating (3.9) at the inner and outer edge of the plasma respectively yields

(3.18)\begin{gather} \varOmega C_{\rm in} + \xi_{\rm in}^2 [C_{\rm in}G(\xi_{\rm in}, \xi_{\rm in}) + C_{\rm out}G(\xi_{\rm in}, \xi_{\rm out})] = 0, \end{gather}

(3.19)\begin{gather}C_{\rm out}\left[ \varOmega - \xi_{\rm in}^2 \ln\left( \frac{\xi_{\rm out}}{\xi_{\rm in}}\right) \right] - \xi_{\rm in}^2\left[ C_{\rm in} G(\xi_{\rm out}, \xi_{\rm in} )+ C_{\rm out} G(\xi_{\rm out}, \xi_{\rm out} )\right] = 0. \end{gather}

Eliminating $C_{\rm in}$ and $C_{\rm out}$ between the two equations gives a quadratic equation for $\varOmega$

(3.20)\begin{equation} \left[ \frac{\varOmega}{\xi_{\rm in}^2} - \ln\left(\frac{\xi_{\rm out}} {\xi_{\rm in}} \right) - G(\xi_{\rm out}, \xi_{\rm out}) \right]\left[ \frac{\varOmega}{\xi_{\rm in}^2} + G(\xi_{\rm in}, \xi_{\rm in}) \right] - G(\xi_{\rm in}, \xi_{\rm out})^2 = 0. \end{equation}

where the relation $G(\xi _{\rm in}, \xi _{\rm out}) = G(\xi _{\rm out}, \xi _{\rm in})$ has been used.

Physically the dispersion relation describes the interaction between two charge density perturbations. There is one at $\xi = \xi _{\rm in}$, where the step function $n^{(0)} (\xi ) \xi ^2$ rises, which we labelled $C_{\rm in}$. The second charge density perturbations is labelled $C_{\rm out}$ and is located at $\xi = \xi _{\rm out}$, where the step function falls back to zero. In the dispersion equation (3.20), the interaction between the two density perturbations is represented by $G(\xi _{\rm in}, \xi _{\rm out})$. When $\xi _{\rm out} - \xi _{\rm in}$ is large compared with unity, $G(\xi _{\rm in}, \xi _{\rm out}) \approx \exp [-(\xi _{\rm out} - \xi _{\rm in})]$ is exponentially small and the interaction between the two density perturbations is negligible. Equation (3.20) then simplifies yielding the two modes

(3.21)\begin{gather} \varOmega_{\rm in} ={-}\frac{\xi_{\rm in}^2}{2} G(\xi_{\rm in}, \xi_{\rm in}) \quad \text{for} \ C_{\rm in} \neq 0 ,C_{\rm out} = 0 , \end{gather}

(3.22)\begin{gather}\varOmega_{\rm out} = \frac{\xi_{\rm in}^2}{2} \left[\ln \left(\frac{\xi_{\rm out}}{\xi_{\rm in}}\right) + G(\xi_{\rm out}, \xi_{\rm out}) \right] \quad \text{for} \ C_{\rm in} = 0 ,C_{\rm out} \neq 0. \end{gather}

When $\xi _{\rm out}-\xi _{\rm in}$ is not large compared with unity, $G(\xi _{\rm in},\xi _{\rm out})$ is not exponentially small and must be retained in (3.20). The two density perturbations interact significantly and each mode involves both perturbations. The solution to (3.20) can be expressed in terms of the frequencies $\varOmega _{\rm in}$ and $\varOmega _{\rm out}$:

(3.23)\begin{equation} \varOmega_{{\pm}} = \tfrac{1}{2}\left( \varOmega_{\rm in} + \varOmega_{\rm out} \pm \sqrt{(\varOmega_{\rm in} - \varOmega_{\rm out})^2 - 4\xi_{\rm in}^4G(\xi_{\rm in}, \xi_{\rm out})^2} \right). \end{equation}

The right-hand side of (3.23) contains an imaginary part if $(\varOmega _{\rm in} - \varOmega _{\rm out})^2 < 4\xi _{\rm in}^4G(\xi _{\rm in}, \xi _{\rm out})^2$. In this case, the two solutions $\varOmega _{\pm }$ correspond to a growing and a damped diocotron mode. The growth rate is largest when the two frequencies $\varOmega _{\rm in}$ and $\varOmega _{\rm out}$ are nearly resonant and the interaction strength $4\xi _{\rm in}^4G(\xi _{\rm in}, \xi _{\rm out})^2$ is large.

An interesting property of these two modes is the respective sign of the wave energies. If the frequencies are both positive, the energy associated with density perturbation $C_{\rm in}$ has positive energy, and the second mode, density perturbation $C_{\rm out}$, has negative energy. In Appendix B we calculate the energy required to create the two density perturbations and find the signs of the energy associated with the two perturbations to be opposite. For the growing mode, the energy flow is from the negative energy perturbation to the positive energy perturbation, and for the damped mode, the energy flow is reversed. Such an interaction, also known as reactive instability, is known to happen both in plasma physics and fluid dynamics (Cairns Reference Cairns1979). The direction in which the energy flows is associated with the phase shift between the two perturbations. This phase shift is given by the expression

(3.24)\begin{align} \tan^{{-}1} \left[ \frac{{\rm Im}(C_{\rm in}/C_{\rm out})}{{\rm Re}(C_{\rm in}/C_{\rm out})} \right] ={-}\frac{{\rm Im}(\varOmega)}{{\rm Re}(\varOmega) + \xi_{\rm in}^2 G(\xi_{\rm in},\xi_{\rm in})}. \end{align}

If the mode is growing ($\varGamma = {\rm Im}(\varOmega )>0$), the outer perturbation is ahead of the inner perturbation and vice versa. The ratio between the amplitudes of the two density perturbations for one distinct mode is given by the expression

(3.25)\begin{equation} \frac{C_{\rm in}}{C_{\rm out}} = \frac{\xi_{\rm in}^2 G(\xi_{\rm in},\xi_{\rm out})}{\varOmega + \xi_{\rm in}^2 G(\xi_{\rm in},\xi_{\rm in})}. \end{equation}

Figure 2 shows the changing mode structure for a decreasing sequence of the plasma thickness $\xi _{\rm out}-\xi _{\rm in}$. The inner and outer conductors are located at $\xi _{\rm wire} = 1$ and $\xi _{\rm wall} = 8$. The thickest plasma in figure 2(a) supports two modes consisting of separate density perturbations at the respective plasma edges. For the slightly thinner plasma in figure 2(b) a small density perturbation on the opposite edge in addition to the dominant one becomes apparent. This trend continues in figure 2(c). Note also that the difference between the frequencies of the two modes decreases. Finally, figure 2(d) shows the case of an unstable mode. The frequency becomes complex and the two solutions turn into their respective complex conjugate, a growing and a damped solution.

Figure 2. The normalised initial density distribution for a pure electron plasma (grey solid line) and the density perturbation at the inner plasma edge (red dashed line) as well as the outer edge (blue dotted line) plotted over the scaled radius $\xi = kr$. The thickness of the initial density distribution decreases from the upper left to the lower right panel. For the thinnest case in panel (d), the frequency becomes complex. Both the damped (purple dashed line) and the growing mode (blue dotted line) are shown.

To understand the dependence on the plasma geometry, that is, on $r_{\rm in}$ and $r_{\rm out}$, we allow these quantities to vary, holding $r_{\rm wire}$, $r_{\rm wall}$ and $k$ constant. This translates to holding $\xi _{\rm wire}$ and $\xi _{\rm wall}$ constant and varying $\xi _{\rm in}$ and $\xi _{\rm out}$ independently. Figure 3 shows a contour plot of ${\rm Im}[\varOmega (\xi _{\rm in}, \xi _{\rm out}, \xi _{\rm wire}, \xi _{\rm wall})]$, where the x-axis is $\xi _{\rm in}$ and the y-axis is the difference $\xi _{\rm out}-\xi _{\rm in}$. Again, the boundaries are fixed at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$.

Figure 3. Contour plot of the growth rate ${\rm Im}(\varOmega )$ of discrete modes, where the abscissa is $\xi _{\rm in}$ and the ordinate is $\xi _{\rm out}-\xi _{\rm in}$. The inner and outer cylindrical conductors are located at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$, respectively. The white line marks the region in which the outer radius would be smaller than the inner radius $\xi _{\rm out}<\xi _{\rm in}$.

As was discussed previously, good interaction between the two density perturbations requires that $\xi _{\rm out}-\xi _{\rm in}$ is not too large. This is in agreement with the vertical range of instability in figure 3 up to $\xi _{\rm out}-\xi _{\rm in} \approx 1$. Likewise, we expect that neither density perturbation can be too near to one of the conducting boundaries. Otherwise, the boundaries would shield the interaction potential. This is consistent with the gap between the last contour of finite growth rate and the ordinate on the left as well as the white solid line on the right in figure 3.

The unscaled growth rate is given by the expression

(3.26)\begin{equation} \gamma = \frac{2{\rm \pi} c^2q_{\rm e}n_{0{\rm e}}}{Ik} {\rm Im}[\varOmega(\xi_{\rm in}, \xi_{\rm out}, \xi_{\rm wire}, \xi_{\rm wall})]. \end{equation}

The growth rate varies linearly with the density and inversely with the current. The scaling of the growth rate with the wavelength is not directly apparent since the scaled distance $\xi =kr$ is a function of the wavelength. Consider a constant wavenumber $k_0$ with the corresponding growth rate $\gamma (k_0) = \alpha \,{\rm Im}[\varOmega (k_0r_{\rm in}, k_0r_{\rm out}, k_0r_{\rm wire}, k_0r_{\rm wall})]$. Changing the wavenumber $k$ while keeping the position of the boundary and the plasma edges fixed results in a change of the growth rate according to the ratio

(3.27)\begin{equation} \frac{\gamma(k)}{\gamma(k_0)} = \frac{k_0}{k} \frac{{\rm Im}[\varOmega(kr_{\rm in}, kr_{\rm out}, kr_{\rm wire}, kr_{\rm wall})]}{{\rm Im}[\varOmega(k_0r_{\rm in}, k_0r_{\rm out}, k_0r_{\rm wire}, k_0r_{\rm wall})]}. \end{equation}

The pre-factor of $k_0/k$ on the right-hand side of (3.27) comes from the $k$-dependence in the scaling factor $\alpha$. The change of the growth rate with wavelength is shown in figure 4.

Figure 4. Growth rate for different wavenumbers $k$ normalised by the growth rate for the reference wavenumber $k_0 = r_{\rm wire}^{-1}$.

Instability is limited to wavenumbers with $k/k_0 < 1.3$, where $k_0 = r_{\rm wire}^{-1}$. Again this restriction comes from the requirement that $\xi _{\rm out}-\xi _{\rm in} = k(r_{\rm out}-r_{\rm in})$ is not large compared with unity. The variation of this quantity in figure 4 comes from the variation of $k$, rather than the variation of $r_{\rm out}-r_{\rm in}$.

For an infinitely long plasma, the allowed wavenumbers $k$ lie on a continuum. That is not the case for the model of a plasma in a large-aspect-ratio dipole trap which implies a periodic boundary condition. Consequently, the allowed wavenumbers are restricted to the discrete values $k=n/R$, where $n$ is a positive integer and $R$ is the radius of the coil. This restriction suggests the possibilities of eliminating unstable modes by choosing $R < (r_{\rm out}-r_{\rm in})/1.3$. Unfortunately, this is inconsistent with our model of a large-aspect-ratio dipole. Nevertheless, this observation does suggest that the opposite limit, a dipole trap with near-unity aspect ratio, might have interesting stability properties.

3.4 Numerical results

For a general density profile, a numerical solution of the mode equation or of the Van Kampen equation is necessary. For solutions of the mode equation, we use the shooting method. This method relies on an initial guess for the complex scaled frequency $\varOmega$. The boundary condition requires that the mode potential is zero at the inner boundary $\xi _{\rm wire}$ and we choose a finite initial slope of the potential at that point. The specific value of the slope does not matter since the mode equation is linear. These two initial conditions allow the equation to be integrated outwards to the conductor at $\xi _{\rm wall}$. The outer boundary condition requires the potential to vanish which is generally not the case for the first guess of $\varOmega$. A root finder that is based on Powell's method (Powell Reference Powell1964) is used to find the complex value of $\varOmega$, for which the boundary condition is satisfied. This value is the correct complex frequency of the mode.

There is a subtlety in the integration of the mode equation from the inner to the outer boundary. To solve (2.14), we divided by the resonance factor $\omega - k ({\partial H^{(0)}}/{\partial p_z})$, which can be zero for some value of $\xi$ if $\omega$ is real. The contour of integration must then be deformed to stay on the same side of the singularity. We avoid this complication by solving the mode equation only for growing modes. In this case the equation can be integrated along the real $\xi$-axis. In subsequent plots of growth rates we use zero value whenever the shooting method does not return a positive value.

In order to solve the eigenvalue problem (3.9) for a general density profile we follow Schecter et al. (Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000). The radius is discretised such that (3.9) can be rewritten in matrix form. For every grid point we obtain one eigenvalue, the scaled frequency, and one eigenfunction, the density perturbation. As discussed in § 3.2 this method yields a discretised representation of all the modes without recourse to analytic continuation. We focus on the discrete unstable mode, which corresponds to the mode obtained by solving the mode equation.

Consider the Gaussian distribution in (2.21). Solving the mode equation (2.22) yields the scaled growth rate $\varGamma = {\rm Im}(\varOmega (\xi _{\rm in}, \xi _{\rm out}, \xi _{\rm wire}, \xi _{\rm wall})$. Figure 5 shows a contour plot of this growth rate in the plane $(\Delta \xi, \xi _{0})$ for fixed values $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$. This figure illustrates the dependence of the growth rate on the plasma geometry for fixed wavenumber and fixed boundaries. The same results are obtained when solving the Van Kampen equation (3.9), selecting only the discrete unstable mode. The result for a Gaussian density profile agrees qualitatively with the analytic result for the profile given by (3.15). Some of the results in figure 5 resemble density profiles with a non-negligible value in the close proximity of a boundary. If a large fraction of the plasma is close to either boundary we expect a quick loss of particles and a consequent truncation of the density profile.

Figure 5. Contour plot of the growth rate ${\rm Im}(\varOmega )$ of discrete modes for a Gaussian equilibrium density profile that is centred around $\xi _{0}$ and has the width $\Delta \xi$ as defined in (2.21). The inner and outer cylindrical conductors are at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$, respectively. The hatched areas highlight the density profiles for which the density in close proximity of a boundary is beyond 1 % of the maximum density.

3.5 Finite-temperature effects

In this section, we allow the parameter $k^2\lambda _{\rm D}^2$ to be finite but retain zero value for the parameter $\eta = 0$ in mode equation (2.22). This is the case of a warm pure electron plasma. Using the shooting method, we solve the mode equation for the Gaussian equilibrium density distribution given in (2.21), with the scaled values $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$. Figure 6 shows a plot of the scaled growth rate vs $k^2\lambda _{\rm D}^2$ for four values of the Gaussian thickness $\Delta \xi$. In all four cases, the scaled growth rate decreases to zero as $k^2\lambda _{\rm D}^2$ increases towards unity, demonstrating a stabilising influence of the curvature and gradient-B drifts on the diocotron mode. As mentioned before, only solutions with ${\rm Im}(\varOmega )>0$ can be accepted.

Figure 6. The growth rate of the diocotron mode is plotted over $k^2\lambda _{\rm D}^2$. The density profile of the plasma is given by the Gaussian in (2.21), centred in between the conducting surfaces. The thickness of the lines corresponds to the thickness $\Delta \xi$ of the density profile. The velocity distribution is a Maxwellian. The inner and outer cylindrical conductors are at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$, respectively.

The stabilisation might appear to be due to Debye screening of the potential perturbation, but note that the mode propagates perpendicular to the magnetic field. The plasma cannot rearrange itself freely such that it shields the potential perturbation. The factor $k^2\lambda _{\rm D}^2$ appears only in front of the gradient-B and curvature drift. Instead, we believe that the stabilisation is due to the spread in these drifts, not to the mere existence of a drift velocity.

The spread comes from the Maxwellian velocity distribution. Thus, we ask what happens if the equilibrium velocity distribution is replaced by $f^{(0)} = g(p_z)\delta (v_{\parallel }-v_{\parallel 0})\delta (v^2_{\perp }-v^2_{\perp 0})$. This distribution is a valid equilibrium, since the variables $v_{\parallel } = {p_{\theta }}/{mr(p_z)}$ and $v_{\perp } = \mu |{2I}/{mcr(p_z)}|$ are both constants of the motion in that $p_{\theta }$, $p_z$ and $\mu$ are all constants of the motion. In Appendix C, we obtain the mode equation for this distribution and solve for the growth rate assuming the same Gaussian density distribution as used for figure 6. The velocity integral in the mode equation is evaluated trivially because of the delta function velocity dependence. As can be seen in figure 7 we find no evidence of stabilisation for large $k^2\hat {\lambda }_{\rm D}^2$ when there is no spread in the drift velocities. Here, we define the characteristic length scale $\hat {\lambda }_{\rm D} = \sqrt {W/4{\rm \pi} q_{\rm e}^2 n_{0{\rm e}}}$ which replaces the Debye length. Instead of the temperature $T$, it depends on the total kinetic energy $W$. Accordingly, the scaled velocities are defined as $x_{\parallel } = \sqrt {{m}/{W}} \, v_{\parallel }$ and $x_{\perp } = \sqrt {{m}/{W}} \,v_{\perp }$.

Figure 7. The growth rate of the diocotron mode is plotted over $k^2\lambda _{\rm D}^2$. The density profile of the plasma is given by the Gaussian in (2.21), centred in between the conducting surfaces. The thickness of the lines corresponds to the thickness $\Delta \xi$ of the density profile. The velocity distribution is a delta distribution that peaks at $x_{\perp } = x_{\parallel } = 1$. The inner and outer cylindrical conductors are at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$ respectively.

The linear stability analysis does not reveal the mechanism behind the stabilisation. However, figure 6 indicates that the stability threshold for $k^2\lambda _{\rm D}^2$ is a strictly increasing function of the growth rate in the cold limit $\varGamma (k^2\lambda _{\rm D}^2 = 0)$. We suspect that, if the spread in the drifts is comparable to $\varGamma (k^2\lambda _{\rm D}^2 = 0)$, the interaction between the negative and positive energy mode described in § 3.3 is mitigated.

4.1 Neutral warm pair plasma

To identify the interchange mode in its simplest form, we first consider the limit of a neutral pair plasma, where $\eta =1$. There is no diocotron mode in this limit. The mode equation reduces to the form

(4.1)\begin{align} & \frac{1}{\xi} \frac{\partial}{\partial \xi} \xi \frac{\partial \phi^{(1)}}{\partial \xi} -\phi^{(1)} \nonumber\\ & \quad ={-}\frac{\phi^{(1)}}{\sqrt{2 {\rm \pi}}} \int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]}{\varOmega - k^2\lambda_{\rm D}^2\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} \left[ \xi\frac{\partial \tilde{n}^{(0)}_{\rm e}}{\partial \xi} + \tilde{n}^{(0)}_{\rm e}\left(\frac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right) \right] \nonumber\\ & \qquad +\frac{\phi^{(1)}}{\sqrt{2 {\rm \pi}}} \int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \frac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]}{\varOmega + k^2\lambda_{\rm D}^2\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)} \left[ \xi\frac{\partial \tilde{n}^{(0)}_{\rm e}}{\partial \xi} + \tilde{n}^{(0)}_{\rm e}\left(\frac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right) \right]. \end{align}

This equation again can be solved by using the shooting method. Figure 8 shows the scaled growth rate $\varGamma ={\rm Im}(\varOmega )$ vs $k^2 \lambda _{\rm D}^2$. Both the electron and the positron density distributions are given by the Gaussian from (2.21), centred in between the conducting cylindrical boundaries at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$ with different widths $\Delta \xi = \xi _{\rm out}-\xi _{\rm in}$. All curves start from $k^2 \lambda _{\rm D}^2=0$ and reach a peak for $k^2 \lambda _{\rm D}^2$ of order unity, but then decrease with further increase in $k^2 \lambda _{\rm D}^2$. In the next section, analytic explanations for the rapid rise and fall will be obtained.

Figure 8. The growth rate is plotted over $k^2\lambda _{\rm D}^2$ for a quasi-neutral pair plasma ($\eta =1$) with different Gaussian density profiles and a Maxwellian velocity distribution. The Gaussian density profiles have different widths in the range $0.4 < \Delta \xi < 0.7$ and they are centred in between the conducting surfaces.

In the scaled mode equation (2.22) the parameter $k^2 \lambda _{\rm D}^2$ enters as the coefficient of the curvature and gradient-B drifts. The rapid rise in growth rate in figure 8 with increasing $k^2 \lambda _{\rm D}^2$ reflects the central role played in the instability by the curvature and gradient-B drifts. We identify the instability as an interchange mode in that its mechanism involves the curvature and gradient-B drifts acting transverse to a density gradient. The magnetic field has bad curvature, and there is a negative density gradient on the outside of the Gaussian. The growth rate decreases for $k^2\lambda _{\rm D}^2 \gtrsim 1$. Again, this stabilisation for large $k^2 \hat {\lambda }_{\rm D}^2$ disappears if we consider a delta distribution instead of a Maxwellian as shown in figure 9. The definition of the parameter $\hat {\lambda }_{\rm D}$ is the same as in § 3.5. Similar to the diocotron mode, the interchange mode appears to be stabilised by a large spread in the drift velocities.

Figure 9. The growth rate of the interchange mode is plotted over $k^2\hat {\lambda }_{\rm D}^2$ for a quasi-neutral pair plasma ($\eta =1$) and different Gaussian density profiles. The Gaussian density profiles have different widths in the range $0.4 < \Delta \xi < 0.7$ and they are centred in between the conducting surfaces. The velocity distribution is a delta distribution that peaks at $x_{\perp } = x_{\parallel } = 1$.

The growth rate displayed in figure 8 might therefore be interpreted as a competition between the curvature and gradient-B drifts as drivers of the instability and the stabilisation due to the spread in those drifts leading to a maximum growth rate at $k^2 \hat {\lambda }_{\rm D}^2 \approx 1$ and stabilisation for $k^2 \hat {\lambda }_{\rm D}^2 > 5$.

However, the instability is different from the standard description of an interchange instability. In the standard description, the curvature and gradient-B drifts are modelled by a drift due to an artificial gravity. The interchange instability is presented as an analogue of the Rayleigh–Taylor instability for a heavy fluid on top of a light fluid (Guzdar et al. Reference Guzdar, Satyanarayana, Huba and Ossakow1982). The large ion mass produces an ion gravitational drift that is large compared with the electron gravitational drift. For the electron–positron plasma the drifts for the two species are equal in magnitude. In the case of an electron–ion plasma the large ion mass gives rise to a large polarisation drift for the ions in the mode field or, equivalently, a large plasma dielectric constant $\epsilon =1+4{\rm \pi} n M_i c^2/B^2 \gg 1$. Because of the small mass for the electrons and positrons, we have neglected the polarisation drift or, equivalently, assumed that the dielectric constant, $\epsilon =1+4{\rm \pi} n m_{\rm e} c^2/B^2$, is nearly unity. In addition, since the magnetic field is inhomogeneous, the $E \times B$ drift is not incompressible.

An interesting feature of the neutral case is that the real part of the frequency vanishes for all values of $k^2\lambda _{\rm D}^2$ (Kennedy Reference Kennedy2020). The mode is purely growing and does not propagate. To understand this result, note that (4.1) is invariant under interchange of $\varOmega \rightarrow -\varOmega$. Consequently, there is nothing in the equation to choose a direction of propagation for the wave, so the real part of the frequency must be zero. We see this result explicitly in graphs of the ${\rm Re}(\varOmega )$ vs $\eta$ in the next section.

4.2 Higher-order radial modes

Up to this point, no mention has been made of higher-order radial modes. The growth rates in figure 8 are all for the fundamental radial eigenmodes. However, higher-order radial modes occur for $\eta \approx 1$ and sufficiently small $k^2\lambda _{\rm D}^2$. Figure 10 shows plots of the first three radial eigenmodes for $k^2\lambda _{\rm D}^2 = 0.2$ and $\eta = 1$. The eigenmodes are typically labelled by the number of zero crossings. One can see that the higher order mode potentials are oscillatory in the region where the slope of the density distribution is negative.

Figure 10. The initial density distribution at different scaled radii $\xi = kr$ is given by the black dashed line. The black solid curves correspond to the potential perturbation for different mode numbers. The thickness of the lines decreases with increasing mode number.

Likewise, figure 11 shows plots of the growth rates for these eigenmodes vs $k^2\lambda _{\rm D}^2$. Note that the fundamental radial eigenmode is the fastest growing and so is the most important for stability considerations.

Figure 11. The growth rates for different radial eigenmodes of the interchange mode are plotted over $k^2\lambda _{\rm D}^2$ assuming a quasi-neutral pair plasma ($\eta =1$). The unperturbed density profile is a Gaussian with the width $\Delta \xi = 0.7$ that is centred in between the conducting surfaces.

To understand these results intuitively, we rewrite (4.1) in a simpler form. Setting $\eta = 1$ and $\varOmega =i\varGamma$ yields

(4.2)\begin{equation} \frac{1}{\xi} \frac{\partial}{\partial \xi} \xi \frac{\partial \phi^{(1)}}{\partial \xi} -\phi^{(1)} = \frac{\phi^{(1)}}{k^2 \lambda_{\rm D}^2} \left[ g \tilde{n}^{(0)} + f \xi\frac{\partial \tilde{n}^{(0)}}{\partial \xi}\right], \end{equation}

where

(4.3)\begin{gather} g(|a|^2) = \sqrt{\frac{2}{\rm \pi}}|a|^2\int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \dfrac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)^2}{1+|a|^2\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)}, \end{gather}

(4.4)\begin{gather}f(|a|^2) = \sqrt{\dfrac{2}{\rm \pi}}|a|^2\int_0^\infty {\rm d}\kern0.7pt x_{{\perp}}^2 \int_0^\infty {\rm d}\kern0.7pt x_{{\parallel}} \dfrac{\exp\left[-\dfrac{1}{2}(x_{{\perp}}^2+x_{{\parallel}}^2)\right]\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)}{1+|a|^2\left(\dfrac{x_{{\perp}}^2}{2}+x_{{\parallel}}^2\right)}, \end{gather}

As can be seen in figure 12, the function $f$ rises monotonically from $0$ to the maximum value of ${\rm \pi}$ as $|a|^2 = {k^4\lambda _{\rm D}^4}/{\varGamma ^2}$ varies from $0$ to $\infty$. The function $g$ rises monotonically from $0$ to the maximum value of $2$ in the same interval.

Figure 12. Evaluating the two integrals in (4.3) for different values of $|a|^2$ yields the solid curve for $g$ and the dashed curve for $f$. The ratio $g/f$ is given by the dotted line.

First, we look at (4.2) in the limit where $|a|^2 = k^4\lambda _{\rm D}^4/\varGamma ^2 \ll 1$. The denominators in (4.3) can then be approximated by unity and performing the integrals yields $g(|a|^2) = 14|a|^2$ and $f(|a|^2) = 4|a|^2$. The mode equation (4.2) takes the form of a Sturm–Liouville problem (Goedbloed & Poedts Reference Goedbloed and Poedts2004).

(4.5)\begin{equation} \frac{1}{\xi} \frac{\partial}{\partial \xi} \xi \frac{\partial \phi^{(1)}}{\partial \xi} - \phi^{(1)} = \lambda_n \phi^{(1)} \left[ \frac{7}{2} \tilde{n}^{(0)} + \xi\frac{\partial \tilde{n}^{(0)}}{\partial \xi}\right]. \end{equation}

A Sturm–Liouville problem is an eigenvalue problem for which the eigenvalue $\lambda _n = 4k^2\lambda _{\rm D}^2/\varGamma _n^2$ increases with the number of zero-crossings $n$ of the eigenfunction. Hence, the lowest-order radial eigenmode ($n = 0$) has the largest scaled growth rate as was observed in figure 11. Likewise, the rapid increase of the growth rate in figures 8 and 11 as $k^2\lambda _{\rm D}^2 = \lambda _n \varGamma _n^2/4$ increases from zero can be understood in terms of this expression for the eigenvalue. All the modes in figure 8 are the fundamental, but with different values of the width for the Gaussian. The variation from curve to curve arises from changes in the eigenvalue as the width of the Gaussian is varied.

If $\frac {7}{2} \tilde {n}^{(0)} + \xi ({\partial \tilde {n}^{(0)}}/{\partial \xi }) > 0$ in the entire confinement region, the second derivative of the potential perturbation is purely positive (purely negative for $\phi ^{(1)} < 0$). The boundary conditions cannot be satisfied in that case. Consequently, the interchange instability occurs for density profiles with a steep, negative gradient only. The same result was found in the local limit (Mishchenko et al. Reference Mishchenko, Plunk and Helander2018). When $k^2\lambda _{\rm D}^2$ approaches unity we have to go back to (4.2). The condition for instability of the interchange mode then generalises to $\frac {g}{f} \tilde {n}^{(0)} + \xi ({\partial \tilde {n}^{(0)}}/{\partial \xi }) < 0$. The fraction ${g}/{f}$ monotonically decreases from $7/2$ to $2/{\rm \pi}$ as $|a|^2$ goes from zero to infinity.

Looking at regions where the density falls off steeply, such that $({\xi }/{\tilde {n}^{(0)}})({\partial \tilde {n}^{(0)}}/{\partial \xi }) \ll -1$, we can think of the quantity $f/k^2\lambda _{\rm D}^2$ as the eigenvalue $\lambda _n$. The function $f(|a|^2)$ monotonically increases to the maximum value of ${\rm \pi}$. Thus, as $\lambda _n k^2\lambda _{\rm D}^2$ approaches ${\rm \pi}$, the growth rate $\varGamma _n$ is forced to zero and there is no solution for $\lambda _n k^2\lambda _{\rm D}^2 > {\rm \pi}$. This behaviour is seen in figure 8 and is probably related to the spread in drift velocities for the Maxwellian velocity distribution that was discussed before.

4.3 Arbitrary non-neutrality

Next, we investigate the stability of pair plasmas with arbitrary non-neutrality. Figure 13 shows the imaginary part of the scaled growth rate ${\rm Im}(\varOmega )$ vs the density ratio $\eta$ for different values of $k^2\lambda _{\rm D}^2$. Figure 14 shows the corresponding real part of the scaled frequency. All of the curves assume that both the electrons and positrons have a Gaussian density distributions that is given by (2.21).

Figure 13. The growth rate is plotted over the density ratio of positrons to electrons $\eta$. The parameter $k^2\lambda _{\rm D}^2$ increases as the curves gradually change colour from blue to red. The density distribution is given by a Gaussian that is centred in between the conducting cylindrical boundaries at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$ and has the width $\Delta \xi = \xi _{\rm out}-\xi _{\rm in} = 0.4$.

Figure 14. The frequency plotted over the density ratio of positrons to electrons $\eta$. The parameter $k^2\lambda _{\rm D}^2$ increases as the curves gradually change colour from blue to red. The density distribution is given by a Gaussian that is centred in between the conducting cylindrical boundaries at $\xi _{\rm wire}=1$ and $\xi _{\rm wall}=8$ and has the width $\Delta \xi = \xi _{\rm out}-\xi _{\rm in} = 0.4$. We omit the values for which ${\rm Im}(\varOmega ) \leqslant 0$.

As we have seen, the diocotron instability depends crucially on the axial $E \times B$ drift, but is stabilised by the curvature and gradient-B drifts. In contrast, the interchange mode depends crucially on the curvature and gradient-B drifts, and shear associated with the $E \times B$ drift has a stabilising influence. Thus, the growing mode is diocotron-like for small values of $\eta$ and small values of $k^2\lambda _{\rm D}^2$. Likewise the interchange mode dominates for large values of $\eta$ and large values of $k^2\lambda _{\rm D}^2$. For intermediate values of $\eta$ and finite $k^2\lambda _{\rm D}^2$, a mode is formed by the competition between the driving and stabilising influences of the different drifts.

Consistent with this picture, the blue curve in figure 13, for which $k^2\lambda _{\rm D}^2 = 0$, exhibits the largest growth rate for small values of $\eta$ and completely vanishes for $\eta =1$. This is the limit of a pure diocotron mode. With a small but finite temperature corresponding to $k^2\lambda _{\rm D}^2 = 0.01$, the interchange mode becomes apparent for $\eta > 0.7$. The red curve, for which $k^2\lambda _{\rm D}^2 = 1$, illustrates the opposite limit where the diocotron mode is suppressed relative to the interchange mode. The plasma is stable for $\eta =0$ but the growth rate increases as $\eta$ approaches unity. For $k^2\lambda _{\rm D}^2>1$ we see a stabilising effect on the interchange mode due to the relative drift motion of particles with different thermal velocities.

As mentioned in the previous section, the real part of the frequency in figure 14 goes to zero as the plasma becomes quasi-neutral. The mode does not propagate in that limit.