2.1 Thermodynamic equations for compression

In this model, we consider a cylindrical pure electron plasma in thermal equilibrium, confined to a PM trap with an axial magnetic field $\boldsymbol {B} = B\hat {z}$. The state variables of the plasma are its temperature $T$ and rotational (angular) velocity $\omega _r$. We neglect finite-length effects by modelling a plasma of physical length $L_p$ as infinite and periodic, with periodicity length $\mathcal {L}_p = 2L_p$ (Eggleston & O'Neil Reference Eggleston and O'Neil1999). All integrals along the $z$-axis in this paper are understood to be over one periodic cell, $-L_p \leq z \leq L_p$.

Electron dynamics follows from the single-particle (sp) Hamiltonian

(2.1)\begin{equation} \mathcal{H}_{\text{sp}}(r, \boldsymbol{p}; T, \omega_r) = K + q\phi_0\left(r; T, \omega_r\right), \end{equation}

with

(2.2)\begin{equation} K = \frac{1}{2m_e}\left(\,p_r^2 + p_z^2\right) + \frac{1}{2m_e r^2}\left(\,p_{\theta} - \frac{1}{2}m_e\varOmega_c r^2\right)^2, \end{equation}

where $m_e$ is the electron mass, $q=-e$ is the electron charge, $\varOmega _c = {qB}/{m_e c}$ is the (signed) cyclotron frequency, $c$ is the speed of light and $\phi _0(r; T, \omega _r)$ is the mean-field electrostatic potential of the plasma. The electrostatic potential depends on the radial density profile, which, given a fixed number of electrons in the trap, is uniquely determined, and in equilibrium, on the plasma temperature and rotation frequency. The explicit dependence on the thermodynamic state variables following the semicolon (e.g. $\phi _0(r;T,\omega _r)$) is a reminder that we are considering equilibrium quantities that evolve by virtue of the plasma transitioning between thermal equilibrium states during compression.

Though we consider only pure electron plasmas in this paper, our analysis naturally extends to positron plasmas. To facilitate this extension, we use a sign convention that readily encompasses both signs of charge. For a magnetic field pointing along $+\hat {z}$, a positively charged plasma column will rotate in the $-\hat {\theta }$ direction under its $E\times B$ drift. This makes the rotational angular velocity negative, $\omega _r<0$, for charges with $q>0$. Other angular rotation frequencies follow the same convention. To summarize our sign conventions, for $q>0$, we have $\varOmega _c>0$, $\omega _r<0$ and $p_\theta >0$ (when strongly magnetized), while for $q<0$, we have $\varOmega _c<0$, $\omega _r>0$ and $p_\theta <0$ (when strongly magnetized).

The distribution function of the plasma is (using the above sign conventions)

(2.3)\begin{equation} f_0(r, \boldsymbol{p}; T, \omega_r) = \frac{1}{\mathcal{Z}}\exp \left[-\frac{1}{T}\left(\mathcal{H}_{\text{sp}}(r, \boldsymbol{p}; T, \omega_r) -\omega_r p_{\theta}\right)\right]. \end{equation}

The constant $\mathcal {Z}$ is determined by normalizing the distribution function to the total number of electrons in a periodic cell,

(2.4)\begin{equation} \frac{1}{\mathcal{L}_p}\int f_0(r, \boldsymbol{p}; T, \omega_r)\,{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p} = \mathcal{N}_e, \end{equation}

where $\mathcal {N}_e$ is the line density of electrons along the axis.

The equilibrium density $n_0(r; T, \omega _r)$ is given by

(2.5)\begin{equation} n_0\left(r; T, \omega_r\right) =\int f_0(r,\boldsymbol{p};T,\omega_r)\,{\rm d}^3\boldsymbol{p} = n_0(0; T, \omega_r) \exp\left({\psi\left(r; T,\omega_r\right)}\right), \end{equation}

where $n_0(0; T, \omega _r)$ is the on-axis density. The function $\psi (r;T)$ (see Dubin & O'Neil Reference Dubin and O'Neil1999) combines the electrostatic and centrifugal potentials:

(2.6)\begin{equation} \psi(r;T,\omega_r) =-\frac{q}{T}\left(\phi_0(r;T,\omega_r)-\phi_0(0;T,\omega_r)- \frac{m_e\omega_r(\varOmega_c+\omega_r)}{2q}r^2\right). \end{equation}

The Poisson equation can then be written as

(2.7)\begin{equation} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\psi(r; T,\omega_r)}{\partial r}\right) = \frac{4{\rm \pi} q^2n_0(0;T,\omega_r)}{T}\left(\exp\left({\psi(r; T,\omega_r)}\right) - 1 - \delta\right), \end{equation}

with

(2.8)\begin{equation} \delta =-\frac{2m_e\omega_r(\varOmega_c+\omega_r)}{4{\rm \pi} q^2 n_0(0;T,\omega_r)}-1. \end{equation}

The rotating wall perturbs the plasma out of equilibrium:

(2.9a)\begin{gather} \phi_0(r;T,\omega_r)\mapsto\phi_0(r;T,\omega_r)+\delta\phi(t,r,\theta,z), \end{gather}

(2.9b)\begin{gather}f_0(r,\boldsymbol{p};T,\omega_r)\mapsto f_0(r,\boldsymbol{p};T,\omega_r)+\delta f(t,r,\theta,z,\boldsymbol{p}). \end{gather}

Our model avoids explicit calculations involving collisional operators even though the time scale for cooling is much longer than the collision time. This is accomplished by assuming that thermalization via collisions occurs sufficiently rapidly so that the plasma evolves through thermal equilibrium states defined by the average energy $\langle E\rangle$ and average angular momentum $\langle \,p_\theta \rangle$. The temporal evolution of these states is found by matching the same quantities computed for the perturbed plasma after the RW perturbation has been active for a small time $\Delta t$.

Explicitly, consider the energy and angular momentum functionals

(2.10a)\begin{gather} \langle E\rangle [\phi,f] =\frac{1}{\mathcal{N}_e\mathcal{L}_p}\int\left[K+\frac{1}{2}q\phi(\boldsymbol{r})\right] f(t,\boldsymbol{r},\boldsymbol{p})\,{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}, \end{gather}

(2.10b)\begin{gather}\langle \,p_\theta\rangle[\phi,f] = \frac{1}{\mathcal{N}_e\mathcal{L}_p}\int p_\theta f(t,\boldsymbol{r},\boldsymbol{p})\,{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}, \end{gather}

where, since $\phi$ is the self-field, the factor of 1/2 in the energy is necessary to avoid double-counting. In thermal equilibrium, these functionals become ordinary functions of the state variables:

(2.11a)\begin{gather} \langle E\rangle [\phi_0(\boldsymbol{r};T,\omega_r),f_0(\boldsymbol{r},\boldsymbol{p};T,\omega_r)]\equiv \langle E\rangle (T,\omega_r), \end{gather}

(2.11b)\begin{gather}\langle \,p_\theta\rangle[\phi_0(\boldsymbol{r};T,\omega_r),f_0(\boldsymbol{r},\boldsymbol{p};T,\omega_r)]\equiv \langle \,p_\theta\rangle (T,\omega_r). \end{gather}

Our equilibration assumption is

(2.12a)\begin{gather} \langle E\rangle (T+\Delta T,\omega_r+\Delta\omega_r) = \langle E\rangle[\phi(t_0+\Delta t),f(t_0+\Delta t,\boldsymbol{r},\boldsymbol{p})], \end{gather}

(2.12b)\begin{gather}\langle \,p_\theta\rangle (T+\Delta T,\omega_r+\Delta\omega_r) = \langle \,p_\theta\rangle[\phi(t_0+\Delta t),f(t_0+\Delta t,\boldsymbol{r},\boldsymbol{p})], \end{gather}

where $f(t_0+\Delta t)$ is the equilibrium distribution at $t_0$, $f(t_0)\equiv f_0$, evolved with the Vlasov equation to time $t=t_0+\Delta t$ and $\phi (t_0+\Delta t)$ is the corresponding total self-consistent potential. Expanding to linear order and dividing through by $\Delta t$ gives

(2.13a)\begin{gather} \frac{{\rm d}\langle E\rangle}{{\rm d} t} =\frac{1}{\mathcal{N}_e\mathcal{L}_p}\int\left\{\left[K+\frac{1}{2}q\phi_0(r;T,\omega_r)\right] \frac{\partial f}{\partial t}+\frac{1}{2}q\frac{\partial\phi}{\partial t} f_0\right\}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}, \end{gather}

(2.13b)\begin{gather}\frac{{\rm d}\langle \,p_\theta\rangle}{{\rm d} t} = \frac{1}{\mathcal{N}_e\mathcal{L}_p}\int p_\theta \frac{\partial f}{\partial t}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}, \end{gather}

where, on the left-hand side, ${{\rm d}}/{{\rm d} t} = \dot {T}({\partial }/{\partial T}) + \dot {\omega }_r({\partial }/{\partial \omega _r})$.

Integrals of the form

(2.14)\begin{equation} \int g(\boldsymbol{r},\boldsymbol{p}) \frac{\partial f}{\partial t}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p} \end{equation}

can be simplified with the Vlasov equation,

(2.15)\begin{equation} \frac{\partial f}{\partial t} =- \frac{\partial f}{\partial \boldsymbol{r}}\boldsymbol{\cdot} \frac{\partial H}{\partial\boldsymbol{p}} + \frac{\partial f}{\partial \boldsymbol{p}}\boldsymbol{\cdot} \frac{\partial H}{\partial \boldsymbol{r}}, \end{equation}

where

(2.16)\begin{equation} H = K + q\phi_0(r;T,\omega)+q\delta\phi \end{equation}

is the full Hamiltonian. Inserting (2.15) into (2.14) and integrating by parts to pull derivatives off of $f$ yields:

(2.17)\begin{align} \int g(\boldsymbol{r},\boldsymbol{p}) \frac{\partial f}{\partial t}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p} = \int f\left(\frac{\partial g}{\partial \boldsymbol{r}}\boldsymbol{\cdot}\frac{\partial H}{\partial\boldsymbol{p}}+g \frac{\partial^2 H}{\partial\boldsymbol{r}\boldsymbol{\cdot}\partial\boldsymbol{p}} - \frac{\partial g}{\partial\boldsymbol{p}}\boldsymbol{\cdot} \frac{\partial H}{\partial\boldsymbol{r}}-g\frac{\partial^2 H}{\partial\boldsymbol{r}\boldsymbol{\cdot}\partial\boldsymbol{p}}\right){\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}. \end{align}

The second derivatives of $H$ cancel and we proceed to evaluate the remaining expressions with $g=p_\theta$ and $g=K + \frac {1}{2}q\phi _0$.

For $g=p_\theta$, we have

(2.18)\begin{align} & \int p_\theta \frac{\partial f}{\partial t}{\rm d}^3\boldsymbol{r}{\rm d}^3\boldsymbol{p} \nonumber\\ & \quad =\int f\left( - \frac{\partial p_\theta}{\partial\boldsymbol{p}}\boldsymbol{\cdot} \frac{\partial H}{\partial\boldsymbol{r}}\right){\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p} =-\int f \frac{\partial H}{\partial \theta}{\rm d}^3\boldsymbol{r}{\rm d}^3\boldsymbol{p}=-\int f \frac{\partial (q\delta\phi)}{\partial\theta}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}. \end{align}

For $g=K + \frac {1}{2}q\phi _0$, the terms ${\partial K}/{\partial \boldsymbol {r}}\boldsymbol {\cdot } {\partial K}/{\partial \boldsymbol {p}}$ cancel. Furthermore, none of the potential terms depend on $\boldsymbol {p}$, so that only remaining terms are

(2.19)\begin{equation} \int f\left(\frac{1}{2}\frac{\partial (q\phi_0)}{\partial\boldsymbol{r}}\boldsymbol{\cdot} \frac{\partial K}{\partial\boldsymbol{p}}-\frac{\partial K}{\partial\boldsymbol{p}}\boldsymbol{\cdot} \frac{\partial (q\phi_0+q\delta\phi)}{\partial\boldsymbol{r}}\right){\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}. \end{equation}

Since $\phi _0$ only has radial dependence, we find

(2.20)\begin{equation} \int \left(K + \frac{1}{2}q\phi_0\right) \frac{\partial f}{\partial t}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}=-\int f \left(\frac{p_r}{2m_e}\frac{\partial (q\phi_0)}{\partial r}+ \frac{\partial K}{\partial\boldsymbol{p}}\boldsymbol{\cdot}\frac{\partial (q\delta\phi)}{\partial\boldsymbol{r}}\right){\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}. \end{equation}

To find a closed set of equations, we use linear theory to approximate the evolution of $f$ and $\delta \phi$. With $f=f_0+\delta f$ in (2.20), terms with only a single perturbed quantity vanish upon azimuthal integration because the equilibrium is azimuthally symmetric and the perturbations are linear combinations of non-zero azimuthal modes. The same argument applies to the $({\partial \phi }/{\partial t})f_0=({\partial (\delta \phi )}/{\partial t})f_0$ term in (2.13a), which thus also vanishes. Terms that are the product of two perturbations (i.e. both $\delta \phi$ and $\delta f$) can be non-zero after azimuthal integration. There is also one term without any perturbations: $\int f_0(({p_r}/{2m_e})({\partial (q\phi _0)}/{\partial r}))\,{\rm d}^3\boldsymbol {r}\,{\rm d}^3\boldsymbol {p}$; this terms vanishes because $f_0$ is Gaussian in $p_r$, making the full integrand odd in $p_r$.

Putting everything together, we have

(2.21a)\begin{gather} \frac{{\rm d}\langle E\rangle}{{\rm d} t} =-\frac{1}{\mathcal{N}_e\mathcal{L}_p}\int \delta f\frac{\partial K}{\partial\boldsymbol{p}}\boldsymbol{\cdot} \frac{\partial (q\delta\phi)}{\partial\boldsymbol{r}}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}, \end{gather}

(2.21b)\begin{gather}\frac{{\rm d}\langle \,p_\theta\rangle}{{\rm d} t} =-\frac{1}{\mathcal{N}_e\mathcal{L}_p}\int\delta f\frac{\partial (q\delta\phi)}{\partial\theta}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p} . \end{gather}

Note that ${\partial K}/{\partial \boldsymbol {p}} = {\partial \mathcal {H}_{{\rm sp}}}/{\partial \boldsymbol {p}}=\dot {\boldsymbol {r}}$ since the potential term has no momentum dependence. Thus, ${\partial K}/{\partial \boldsymbol {p}}\boldsymbol {\cdot } {\partial (q\delta \phi )}/{\partial \boldsymbol {r}} = \dot {r}({\partial (q\delta \phi )}/{\partial r}) + \dot {z}({\partial (q\delta \phi )}/{\partial z}) + \dot {\theta }({\partial (q\delta \phi )}/{\partial \theta })$. We neglect the first term because it is proportional to the inverse of the compression time scale ($\dot {r}\partial \delta \phi /\partial r\sim \delta \phi /\tau _{{\rm comp}}\ll \dot {\theta }\partial \delta \phi /\partial \theta \sim \omega _r\delta \phi$). We neglect the second term because it describes the heating of particles resonant with the RW wave as their distribution flattens in a process akin to nonlinear Landau damping, which will not occur assuming that thermalization is rapid enough to maintain thermal equilibrium during compression. Since the plasma rotates coherently in equilibrium, we have $\dot {\theta }=\omega _r={\rm const.}$ and we can move $\omega _r$ outside the integral.

Electrons cool via cyclotron emission and are heated by a thermal reservoir at temperature $T_{{\rm res}.}$. We account for this by subtracting the term $({1}/{\tau _R})(T - T_{{\rm res}.})$ in the EOM for $\langle E\rangle$, where $\tau _R = {9 m_e c^3}/{8q^2\varOmega _c^2}$ is the time scale for cyclotron radiation (Beck, Fajans & Malmberg Reference Beck, Fajans and Malmberg1996).

Lastly, we add a phenomenological drag torque $\tau _{{\rm dr}}$ and its corresponding heating or cooling term $P_{{\rm dr}}$. If the drag torque can be modelled with a potential perturbation in a manner similar to the RW torque, then $P_{{\rm dr}}=\tau _{{\rm dr}}\omega _r$, but in more general cases (such as collisions with background gas particles deGrassie & Malmberg Reference deGrassie and Malmberg1980), the two quantities are independent.

Our full system of equations is then

(2.22a)\begin{gather} \frac{{\rm d}\langle E\rangle}{{\rm d} t} = \langle\tau\rangle\omega_r - \frac{1}{\tau_R} \left(T - T_{{\rm res}.}\right)+P_{{\rm dr}}, \end{gather}

(2.22b)\begin{gather}\frac{{\rm d}\langle \,p_\theta\rangle}{{\rm d} t} = \langle\tau\rangle+\tau_{{\rm dr}}, \end{gather}

(2.22c)\begin{gather}\langle\tau\rangle =-\frac{1}{\mathcal{N}_e\mathcal{L}_p}\int\delta f\frac{\partial (q\delta\phi)}{\partial\theta}{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}, \end{gather}

where the last line introduces the average torque acting on the plasma.

We pause here to remark about the signs of various terms. Recall from our sign conventions in § 2.1 that electrons ($q<0$) have $\omega _r>0$ and $\langle \,p_\theta \rangle \sim qB\sum _i r_i^2<0$. For the plasma to compress, we need $\langle \,p _{\theta }\rangle$ to become less negative (smaller $\sum _i r_i^2$), which requires a positive torque according to (2.22b). A positive torque and rotation frequency in turn means the first term of (2.22a) is positive, leading to energy gain and heating. This is the expected behaviour since the RW torque does work on the plasma to compress it. For $q>0$ (e.g. positrons), the signs are reversed: $\omega _r<0$, $\langle \,p_\theta \rangle >0$, $\langle \tau \rangle <0$, but as expected, the first term of (2.22a) remains positive since the RW still does work to compress the plasma.

Next, we rewrite (2.22a) and (2.22b) using the dimensionless variables

(2.23a)\begin{gather} D = 1 - \frac{\omega_r}{\omega_{{\rm RW}}}, \end{gather}

(2.23b)\begin{gather}\varTheta = \frac{T}{m_e v_{\phi, 0}^2}, \end{gather}

(2.23c)\begin{gather}\hat{t} = \frac{1}{\tau_R}t, \end{gather}

where $\omega _{{\rm RW}}$ is the angular rotation frequency of the rotating wall (the sign of $\omega _{{\rm RW}}$ will depend on the direction of rotation, but to compress electrons, we need $\omega _{{\rm RW}}>0$), $v_{\phi, 0} = {|\omega _{{\rm RW}}|}/{k_0}$ is the phase velocity of the fundamental axial RW mode and

(2.24)\begin{equation} k_0= \frac{2{\rm \pi}}{\mathcal{L}_p} \end{equation}

is the wavenumber of the fundamental mode. The variable $D$ is a normalized detuning between the plasma rotation and RW frequencies. Variable $D$ is bounded above by unity so long as angular velocity of the plasma is parallel to that of the RW. Variable $D$ is bounded below by $0$ so long as $|\omega _r| < |\omega _{{\rm RW}}|$, which is the case analysed here. When $D\approx 0$, the plasma has (practically) fully compressed because $\omega _r\approx \omega _{{\rm RW}}$. The variable $\varTheta = {T}/{m_e v_{\phi, 0}^2}$ is the ratio between the thermal energy of an electron and the kinetic energy of an electron that is in-phase with the fundamental RW mode. The variables $D$ and $\varTheta$ are in bijective correspondence with $\langle \mathcal {H}\rangle$ and $\langle \,p _{\theta }\rangle$ as we show in the next section. Hereafter, we take $D$ and $\varTheta$ to be the state variables.

In terms of these dimensionless variables, our EOM governing the energy and angular momentum dynamics of the plasma are

(2.25)\begin{equation} \begin{bmatrix} \dot{\varTheta}\\ \dot{D} \end{bmatrix} = {\mathcal{M}}^{-1}\begin{bmatrix} (1-D)\hat{\tau}(D,\varTheta)-(\varTheta-\varTheta_{{\rm res}.})+\hat{P}_{{\rm dr}}\\ \hat{\tau}+\hat{\tau}_{{\rm dr}} \end{bmatrix}, \end{equation}

where

(2.26)\begin{gather} \varTheta_{{\rm res}.} = \frac{T_{{\rm res}.}}{m_e v_{\phi, 0}^2}, \end{gather}

(2.27)\begin{gather}\hat{P}_{{\rm dr}} = \frac{\tau_R}{m_e v_{\phi, 0}^2}P_{{\rm dr}}, \end{gather}

(2.28)\begin{gather}\hat{\tau} = \left(\frac{\tau_R \omega_{{\rm RW}}}{m_e v_{\phi, 0}^2}\right)\langle\tau\rangle, \end{gather}

with the last line defining the dimensionless torque (and a similar relation holds between $\hat {\tau }_{{\rm dr}}$ and $\tau _{{\rm dr}}$). The partial derivative matrix needed to change variables is

(2.29)\begin{equation} {\mathcal{M}}=\begin{bmatrix} \dfrac{\partial \hat{E}}{\partial\varTheta} & \dfrac{\partial\hat{E}}{\partial D}\\ \dfrac{\partial \hat{p}_\theta}{\partial\varTheta} & \dfrac{\partial\hat{p}_\theta}{\partial D} \end{bmatrix}, \end{equation}

where

(2.30)\begin{gather} \hat{E}=\frac{\langle E\rangle}{m_e v_{\phi, 0}^2}, \end{gather}

(2.31)\begin{gather}\hat{p}_\theta = \frac{\langle \,p_{\theta}\rangle\omega_{{\rm RW}}}{m_e v_{\phi, 0}^2}. \end{gather}

To arrive at a closed set of equations, we (i) express $\hat {E}$, $\hat {p}_\theta$ and their partial derivatives in terms of $\varTheta$ and $D$, and (ii) explicitly evaluate the perturbations $\delta \phi$ and $\delta f$ and substitute them into the torque (2.22c).

2.2 Expressions for the energy and angular momentum

We begin with the expression for energy in equilibrium:

(2.32)\begin{equation} \langle E\rangle (T,\omega_r) =\frac{1}{\mathcal{N}_e\mathcal{L}_p} \int\left[K+\frac{1}{2}q\phi_0(r)\right] f_0(r,\boldsymbol{p})\,{\rm d}^3\boldsymbol{r}\,{\rm d}^3\boldsymbol{p}. \end{equation}

Completing the square on $p_{\theta }$ in the exponent of $f_0$, (2.3) allows us to perform the momentum integrals as simple Gaussian expectations. We take advantage of the usual Gaussian expectation results ${\langle \,p_r^2\rangle _G}/{2m_e}={\langle \,p_z^2\rangle _G}/{2m_e}={\langle (\,p_{\theta }-\langle \,p _{\theta }\rangle )^2\rangle _G}/{2m_e r^2}={T}/{2}$ and $\langle \,p _{\theta }\rangle _G = \frac {1}{2}m_e (\varOmega _c+2\omega _r)r^2$, where the $G$ subscript reminds us that we average over a Gaussian distribution here. Simplifying the resulting expression leaves

(2.33)\begin{equation} \langle E\rangle = \frac{3T}{2} + \frac{\displaystyle\int \left(\frac{1}{2}m_e \omega_r^2 r^2 + \frac{1}{2}q\phi_0\right) {\rm e}^\psi r\,{\rm d} r}{\displaystyle\int {\rm e}^\psi r\,{\rm d} r}, \end{equation}

where $\psi$ was introduced in (2.6), and we performed the trivial $\theta$ and $z$ integrals.

We now introduce the Debye length, normalized radial coordinate, ratio of plasma rotation to cyclotron frequency and a convenient variable $\varDelta$:

(2.34a)\begin{gather} \lambda_D^2 = \frac{T}{4{\rm \pi} q^2n_0(r=0; T, \omega_r)}, \end{gather}

(2.34b)\begin{gather}\rho \equiv r/\lambda_D, \end{gather}

(2.34c)\begin{gather}\chi \equiv \omega_r/\varOmega_c,\quad |\chi|\ll 1, \end{gather}

(2.34d)\begin{gather}\varDelta \equiv \frac{2\omega_r(\varOmega_c+\omega_r)m_e\lambda_D^2}{T}=-(1+\delta). \end{gather}

The normalization $\int n_0\,{\rm d}^3\boldsymbol {r} = N_e\implies 2{\rm \pi} n_0(r=0)\int {\rm e}^\psi r\,{\rm d} r = \mathcal {N}_e$ can be rewritten with the help of the Debye length and normalized radial coordinate

(2.35)\begin{equation} \int {\rm e}^\psi\rho\,{\rm d}\rho = 2q^2\mathcal{N}_e/T. \end{equation}

Substituting all of these in (2.33) and exchanging the remaining $\phi _0$ in favour of $\psi$ gives

(2.36)\begin{align} \langle E\rangle = \frac{3T}{2} + \frac{q\phi_0(r=0; T, \omega_r)}{2} + \frac{T^2}{4q^2\mathcal{N}_e}\left(-\int \psi \,{\rm e}^\psi \rho\,{\rm d}\rho + \frac{\varDelta}{4}\left[\frac{1+3\chi}{1+\chi}\right]\int \rho^3 \,{\rm e}^\psi \,{\rm d}\rho\right). \end{align}

To calculate the potential at the centre of the plasma column, note that the electric field can be written as

(2.37)\begin{equation} E(r) = \frac{4{\rm \pi} q}{r}\int_0^rn_0(r')r'\,{\rm d} r'. \end{equation}

Setting $\phi _0(R_W)=0$ (conducting wall at the boundary), we have

(2.38)\begin{equation} \phi_0(r=0; T, \omega_r) = \int_0^{R_w}\frac{4{\rm \pi} q}{r}\int_0^r n_0(r')r'\,{\rm d} r'\,{\rm d} r. \end{equation}

This can be rewritten in terms of $\psi$ and $\rho$ as

(2.39)\begin{equation} \frac{q\phi_0(r=0; T, \omega_r)}{2} = \frac{T}{2}\int_0^{\rho_w}\frac{{\rm d}\rho}{\rho} \int_0^\rho \exp({\psi(\rho')})\rho'\,{\rm d}\rho', \end{equation}

where $\rho _w=R_w/\lambda _D$. This completes our expression for the energy in terms of $T$ and $\omega _r$ (and hence in terms of $\varTheta$ and $D$). The remaining implicit dependence of $\psi$ on $T$ can be uncovered by searching for $\delta$ such that the solution of the Poisson equation (2.7) satisfies the normalization $\int {\rm e}^\psi \rho \,{\rm d}\rho = 2q^2\mathcal {N}_e/T$.

The calculation of the angular momentum proceeds in a similar manner. We can take advantage of our previous result for the expectation of $p_{\theta }$ under Gaussian integration:

(2.40)\begin{equation} \langle \,p_{\theta}\rangle = \frac{\displaystyle\int\frac{1}{2}m_e (\varOmega_c+2\omega_r)r^2 \,{\rm e}^\psi r\,{\rm d} r}{\displaystyle\int {\rm e}^\psi r\,{\rm d} r}. \end{equation}

Substituting our dimensionless variables gives

(2.41)\begin{equation} \langle \,p_{\theta}\rangle = \frac{T^2\varDelta}{8q^2\mathcal{N}_e} \frac{1+2\chi}{\varOmega_c\chi(1+\chi)}\int\rho^3 \,{\rm e}^\psi\,{\rm d}\rho. \end{equation}

Our final expressions for the dimensionless energy and angular momentum are thus

(2.42a)\begin{align} \hat{E} & = \frac{3}{2}\varTheta + \frac{1}{2}\varTheta \int_0^{\rho_w}\frac{{\rm d}\rho}{\rho}\int_0^\rho \exp({\psi(\rho')})\rho'\,{\rm d}\rho' \nonumber\\ & \quad +\frac{\alpha\varTheta^2}{4}\left(-\int \psi \,{\rm e}^\psi \rho\,{\rm d}\rho +\frac{\varDelta}{4}\left[\frac{1+3\chi}{1+\chi}\right]\int \rho^3 \,{\rm e}^\psi \,{\rm d}\rho\right), \end{align}

(2.42b)\begin{align} \hat{p} & = \frac{\alpha\varTheta^2}{8}\frac{\varDelta}{1-D}\frac{1+2\chi}{1+\chi}\int\rho^3 \,{\rm e}^\psi\,{\rm d}\rho, \end{align}

where $\alpha \equiv {mv_{\phi,0}^2}/{q^2\mathcal {N}_e}$. As noted above, $\varDelta$ and $\varTheta$ are implicitly linked by the Poisson equation. Equation (2.42) can be simplified further using the assumption $|\chi |\ll 1$, as is typically the case in experiments.

Despite the seeming complexity of the expressions (2.42) in dimensionless variables, the underlying physical quantities are straightforward averages over the plasma. So the average energy $\langle E\rangle$ (see (2.33)) is a sum of the average thermal energy $3T/2$, the average rotational kinetic energy $\langle \frac {1}{2}m_e \omega _r r^2\rangle$ and the self-energy $\langle \frac {1}{2}q\phi _0\rangle$. The average angular momentum $\langle \,p_\theta \rangle$ (see (2.40)) is a sum of the kinetic angular momentum $\langle m_e\omega _r r^2\rangle$ and the magnetic field correction $({qB}/{2c})\langle r^2\rangle$. The evolution is driven by the average torque $\langle \tau \rangle$, which is just the average of $\boldsymbol {r}\times \boldsymbol {F}=-q\boldsymbol {r}\times \boldsymbol {\nabla }\phi =-q({\partial \phi }/{\partial \theta })\hat {z}$.

We now turn to the final piece of our model: the derivation of the torque.

2.3 Derivation of the rotating wall torque

The RW torque is found by first solving for $\delta f$ and $\delta \phi$, and then inserting them into (2.22c). In the previous section, the Vlasov distribution on the full six-dimensional phase space was for the thermodynamic equations. In this section, we take advantage of the small Larmor radius and use the guiding centre approximation. Our distribution will thus be over a four-dimensional phase space $(r,\theta,z,v_z)$, with the first three variables denoting the position of the guiding centre and the last being the guiding centre velocity along the magnetic field. See Lifshitz & Pitaevski (Reference Lifshitz and Pitaevski1981) for details on transitioning from the particle to the guiding centre phase space. From now on, all variables are assumed to refer to the guiding centre. The full distribution is

(2.43)\begin{equation} f(t, r, \theta, z, v_z) = f_0(r, v_z; \varTheta, D) + \delta f(t, r, \theta, z, v_z). \end{equation}

The equilibrium distribution function, $f_0(r, v_z)$, can be written as the equilibrium density, $n_0(r; \varTheta, D)$, multiplied by a Maxwellian velocity distribution,

(2.44)\begin{equation} f_0(r, v_z; \varTheta, D) = n_0\left(r; \varTheta, D\right) \frac{1}{\sqrt{2{\rm \pi} v_{th}^2}}\exp\left[-\frac{1}{2}\left(\frac{v_z}{v_{th}}\right)^2\right], \end{equation}

with thermal velocity $v_{th} = v_{\phi,0}\sqrt {\varTheta } = \sqrt {T/m_e}$.

The distribution function, $f(t, r, \theta, z, v_z)$, evolves according to the drift kinetic equation (DKE) (Dubin & O'Neil Reference Dubin and O'Neil1999; Lifshitz & Pitaevski Reference Lifshitz and Pitaevski1981),

(2.45)\begin{equation} \frac{\partial f}{\partial t} + v_z\frac{\partial f}{\partial z} + \frac{c}{|\boldsymbol{B}|^2}\left(\boldsymbol{E}\times\boldsymbol{B}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} f - \frac{q}{m_e}\frac{\partial\phi}{\partial z}\frac{\partial f}{\partial v_z}= \epsilon\delta f, \end{equation}

where $\epsilon$ is an effective collision frequency, $\delta f$ is the perturbation of the distribution function, and $\boldsymbol {E}$ and $\boldsymbol {B}$ are the total electric and magnetic fields. Note that $\boldsymbol {E} = -\boldsymbol {\nabla } \phi (t, r, \theta, z)$ includes the RW potential.

We solve the linearized DKE together with the linearized Poisson equation

(2.46)\begin{equation} \nabla^2\delta\phi(t, r, \theta, z) =-4{\rm \pi} q\delta n(t, r, \theta, z), \end{equation}

where

(2.47)\begin{equation} \delta n = \int \delta f\,{\rm d} v_z, \end{equation}

by expanding in Fourier modes:

(2.48a)\begin{gather} \delta f(t, r, \theta, z, v_z) = \sum_{m, l, s}a_{m, l, s}\delta f_{m, l, s}(r, v_z) \exp\left({{\rm i}\left(m k_0 z + l\theta - st\omega_{{\rm RW}}\right)}\right), \end{gather}

(2.48b)\begin{gather}\delta n(t, r, \theta, z) = \sum_{m, l, s}a_{m, l, s}\delta n_{m, l, s}(r) \exp\left({{\rm i}\left(m k_0 z + l\theta - st\omega_{{\rm RW}}\right)}\right), \end{gather}

(2.48c)\begin{gather}\delta\phi(t, r, \theta, z) = \sum_{m, l, s}a_{m, l, s}\delta\phi_{m, l, s}(r)\exp\left({{\rm i}\left(m k_0 z + l\theta - st\omega_{{\rm RW}}\right)}\right). \end{gather}

The mode functions $\delta \phi _{m, l, s}(r)$ satisfy the boundary condition $\delta \phi _{m, l, s}(R_w)=1$. For convenience, we chose the same expansion coefficients $a_{m, l, s}$ for all three of $\delta \phi$, $\delta f$ and $\delta n$. The mode functions $\delta f_{m, l, s}$ and $\delta n_{m, l, s}$ are then uniquely determined by $\delta \phi _{m, l, s}$ as we show below. The coefficients $a_{m, l, s}$ are determined by requiring that the RW potential at the trap radius satisfies

(2.49)\begin{equation} \delta \phi(R_w) = V(z) \cos\left(\theta - \omega_{{\rm RW}} t\right), \end{equation}

where $V(z) = \varPhi _0$ on the RW electrodes and $V(z) = 0$ elsewhere. We picked this form to simplify our calculations, although one could modify our model to include the sectored nature of RW apparatuses. As such, only the $l = s = \pm 1$ modes are non-zero in our model. The coefficients, $a_{m, l, s}$, are thus given by the usual relation $a_{m, \pm 1, \pm 1} = ({1}/{2\mathcal {L}_p})\int \cos (m k_0 z)V(z)\,{\rm d} z$. Because $V(z)$ is a step function, the coefficients are

(2.50)\begin{equation} a_{m, \pm 1, \pm 1} = \varPhi_0\frac{\sin\left(m k_0 \mathcal{L}_e/2\right)}{m k_0\mathcal{L}_p}, \end{equation}

where $\mathcal {L}_e=2L_e$ is twice the physical electrode length.

Returning to the DKE, we insert the Fourier decompositions and linearize to find

(2.51)\begin{equation} \delta f_{m, l, s}(r, v_z) = \frac{q}{m_e v_{{\rm th}}^2} \left(\frac{mk_0v_z}{s\omega_{{\rm RW}} - l\omega_{E\times B} - m k_0 v_z - {\rm i}\epsilon}\right)f_0(r, v_z; \varTheta, D) \delta\phi_{m, l, s}(r). \end{equation}

In (2.51), we dropped the diamagnetic drift, $({m_e v_{{\rm th}}^2}/{q})({lc}/{rB})({\partial \ln n_0(r; \varTheta, D)}/{\partial r})$, which is one of the terms arising from $\boldsymbol {\nabla }_{\perp }f(t, r, \theta, z, v_z)$ in (2.45). This term is small for the regime of interest here (see, for example, Section IIB of Dubin & O'Neil Reference Dubin and O'Neil1999). In the absence of any diamagnetic drift,

(2.52)\begin{equation} \omega_r = \omega_{E\times B} = \frac{v_{E\times B}}{r} = \frac{|\boldsymbol{E}\times \boldsymbol{B}|}{B^2r}c, \end{equation}

and we use these interchangeably in the rest of the paper.

The density perturbation, $\delta n_{m, l, s}(r)$, is obtained by integrating $\delta f_{m, l, s}(r, v_z)$ over the axial velocity $v_z$. Using the plasma dispersion function $Z$ (Fried & Conte Reference Fried and Conte2015), we have

(2.53)\begin{equation} \delta n_{m, l, s} =-\frac{q}{T}n_0\left(r; T, \omega_r\right)\delta\phi_{m, l, s} \left[1+\frac{v_{{\rm res}}}{v_{{\rm th}}\sqrt{2}}Z \left(\frac{v_{{\rm res}}}{v_{{\rm th}}\sqrt{2}}\right)\right],\end{equation}

where the resonance between the rotating wall and the combined electron rotational and axial motion is given by

(2.54)\begin{equation} v_{{\rm res}} = \frac{s\omega_{{\rm RW}} - l\omega_{E\times B}}{m k_0}. \end{equation}

Inserting this density perturbation into the linearized Poisson equation (2.46) yields an equation for $\delta \phi _{m, l, s}$:

(2.55)\begin{align} & \left\{\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right) - \left(\frac{l^2}{r^2} + \left(m k_0 \right)^2 + \frac{4{\rm \pi} q^2 n_0\left(r; \varTheta, D\right)}{T} \right.\right. \nonumber\\ & \quad \left.\left.\times \left[1+\frac{v_{{\rm res}}}{v_{{\rm th}}\sqrt{2}}Z\left(\frac{v_{{\rm res}}}{v_{{\rm th}} \sqrt{2}}\right)\right]\right)\right\}\delta\phi_{m, l, s} = 0. \end{align}

This is a homogeneous Bessel-like equation with a source term proportional to density manifesting a Debye screening effect.

The solution to this equation together with the relationship between $\delta n_{m, l, s}$ and $\delta \phi _{m, l, s}$ is sufficient to evaluate the torque, but further simplifications can be made. Writing

(2.56)\begin{equation} \delta n_{m, l, s}(r) = b_{m, l, s}(r)\delta \phi_{m, l, s(r)}, \end{equation}

where $b_{m, l, s}$ is defined by (2.53), the torque (2.22c) becomes

(2.57)\begin{align} \langle\tau\rangle & =-\frac{q}{\mathcal{L}_p\mathcal{N}_e}\int \delta n \frac{\partial\delta\phi(t, r, \theta, z)}{\partial\theta}r\,{\rm d} r \,{\rm d}\theta \,{\rm d} z\nonumber\\ & =-\frac{2{\rm \pi} q}{\mathcal{N}_e}\sum_{m, l, s}il \int b_{m, l, s}^*(r)|a_{m, l, s} \delta\phi_{m, l, s}(r)|^2 r\,{\rm d} r\nonumber\\ & =-\frac{8{\rm \pi} q}{\mathcal{N}_e}\sum_{m>0}|a_{m, 1, 1}|^2 \int\Im\left[b_{m, 1, 1}\right]|\delta\phi_{m, 1, 1}(r)|^2 r\,{\rm d} r. \end{align}

In going from the second to the third line of (2.57), we summed over $s = l = \pm 1$, restricted the wavenumber sum to positive $m$ and used the fact that $\Im [b_{m,1,1}]=-\Im [b_{m,-1,-1}]$ (see below).

To explicitly evaluate $\Im [b_{m,1,1}]$, express the plasma dispersion function using the Sokhotski–Plemelj formula

(2.58)\begin{equation} \delta n_{m, l, s}(r) = P\int\frac{g(r, v_z)\delta\phi_{m, l, s}(r)}{v_z - \displaystyle\left(\frac{s\omega_{{\rm RW}} - \ell\omega_{E\times B}}{m k_0}\right)}{\rm d} v_z-{\rm i}{\rm \pi} \delta\phi_{m, l, s} g\left(r, v_z = \frac{s\omega_{{\rm RW}} - \ell\omega_{E\times B}}{m k_0}\right), \end{equation}

where $P$ denotes the Cauchy principal value and

(2.59)\begin{equation} g(r, v_z) =\frac{q}{m_e}\frac{\partial f_0}{\partial v_z}=-\frac{q v_z}{m_e v_{th}^2} f_0(r, v_z; \varTheta, D). \end{equation}

We can now read off $b_{m, l, s}(r)$ from (2.58):

(2.60)\begin{equation} b_{m, l, s} = P\int\frac{g(r, v_z)(r)}{v_z - \displaystyle \left(\frac{s\omega_{{\rm RW}} - \ell\omega_{E\times B}}{m k_0}\right)}{\rm d} v_z-{\rm i}{\rm \pi} g\left(r, v_z = \frac{s\omega_{{\rm RW}} - \ell\omega_{E\times B}}{m k_0}\right), \end{equation}

from which we see

(2.61)\begin{gather} \Im[b_{m,1,1}] =-{\rm \pi} g\left(r,v_z=\frac{\omega_{{\rm RW}}-\omega_{E\times B}}{mk_0}\right), \end{gather}

(2.62)\begin{gather}\Im[b_{m,-1,-1}] =-{\rm \pi} g\left(r,v_z=-\frac{\omega_{{\rm RW}}-\omega_{E\times B}}{mk_0}\right). \end{gather}

Since $g(r,v_z)$ is an odd function of $v_z$, we have $\Im [b_{m,1,1}]=-\Im [b_{m,-1,-1}]$ as stated above.

Recalling that $f_0(\boldsymbol {r}, \boldsymbol {v}; \varTheta, D)$ has the Maxwellian velocity distribution given by (2.44), we obtain the desired expression for the dimensionless torque

(2.63a)\begin{gather} \hat{\tau}\left(D, \varTheta\right) =-4\eta \left(\frac{D}{\varTheta^{3/2}}\right) \sum_{m>0} \hat{\tau}_{m}, \end{gather}

(2.63b)\begin{gather}\hat{\tau}_{m} = \frac{1}{m}\exp\left[-\frac{1}{2m^2}\frac{D^2}{\varTheta}\right] \left|\frac{a_{m,1,1}}{\varPhi_0}\right|^2\langle\delta\phi_{m, 1, 1}^2\rangle, \end{gather}

where we have defined the dimensionless torque strength

(2.64)\begin{equation} \eta = \sqrt{\frac{\rm \pi}{2}}\left(\frac{q\varPhi_0}{m_e v_{\phi, 0}^2}\right)^2 \tau_R\omega_{{\rm RW}}, \end{equation}

and introduced the notation

(2.65)\begin{align} \langle\delta\phi_{m, 1, 1}^2\rangle & = \frac{2{\rm \pi}}{\mathcal{N}_e}\int n_0(r)|\delta\phi_{m, 1, 1}|^2 r\,{\rm d} r \nonumber\\ & = \frac{\displaystyle\int n_0(r)|\delta\phi_{m, 1, 1}|^2 r\,{\rm d} r}{\displaystyle\int n_0(r)r\,{\rm d} r}. \end{align}

Now it is straightforward to substitute the solution of the Poisson equation (2.55) into the torque. This completes our model.