Guiding centre motion for particles in a ponderomotive magnetostatic end plug

Abstract

The Hamiltonian dynamics of a single particle in a rotating plasma column, interacting with an magnetic multipole is perturbatively solved for up to second order, using the method of Lie transformations. First, the exact Hamiltonian is expressed in terms of canonical action-angle variables, and then an approximate integrable Hamiltonian is introduced, using another set of actions and angles, which describe the centre of oscillation for the particle. The perturbation introduces an effective ponderomotive potential, which to leading order is positive. At the second order, the pseudopotential consists of a sum of terms of the Miller form, and can have either sign. Additionally, at second order, the ponderomotive interaction introduces a modification to the particle effective mass, when considering the motion along the column axis. It is found that particles can be axially confined by the ponderomotive potentials, but acquire radial excursions which scale as the confining potential. The radial excursions of the particle along its trajectory are investigated, and a condition for the minimal rotation frequency for which the particle remains radially confined is derived. Last, we comment on the changes to the aforementioned solution to the pseudopotentials and particle trajectory in the case of resonant motion, that is, a motion which has the same periodicity as the perturbation.

1 Introduction

A charged particle interacting with an electromagnetic wave of slowly varying amplitude can experience a ponderomotive force of the Miller type (Gaponov & Miller 1958; Motz & Watson 1967). If the (generalized) particle has one or more internal degrees of freedom (e.g. cyclotron motion), the interaction can be attractive or repulsive, depending on the difference between the wave frequency and the natural frequencies of the internal degrees of freedom (Dodin, Fisch & Rax 2004; Dodin & Fisch 2006).

A particle can experience oscillating fields in its rest frame (or in its gyrocentre frame) if it moves through spatially varying corrugated static fields, which can be electrostatic (Anderegg et al. 1995) or magnetostatic (Rubin, Rax & Fisch 2023) or both. In (Rubin et al. 2023), the leading-order ponderomotive pseudopotential, resulting from the average magnetic field along the particle trajectory, was shown to be independent of the rotation frequency and the cyclotron frequency of the axial field. In addition, this pseudopotential is always positive, leading to a repulsive force away from the perturbation region. Particle motion in electromagnetic fields modelled by $\boldsymbol {E}\boldsymbol {\cdot } \boldsymbol {B}=0$ was investigated by (Ochs & Fisch 2023a), in which a particle drifting in a slab experiences only pseudopotentials of the Miller type, without the always-repulsive leading-order term.

Particle dynamics in these fields is oscillatory in all three degrees of freedom. These oscillations give rise to higher-order ponderomotive potentials, which can be of use as a ponderomotive end plug for open field line magnetic confinement devices, including mirror-type confinement schemes (Baldwin 1977; Gormezano 1979; Post 1987; Ryutov 1988). The end-plug concept considered here utilizes plasma rotation, which is useful in and of itself as a confinement strategy (Lehnert 1971; Bekhtenev et al. 1980; Volosov & Pekker 1981; Hassam 1997; Fetterman & Fisch 2008, 2010a,b; Teodorescu et al. 2010; Fowler, Moir & Simonen 2017; White, Hassam & Brizard 2018; Miller et al. 2023). Following Pastukhov (1974) and Schwartz et al. (2023), even a small additional confining potential may have an appreciable effect on the collisional energy loss rate due to the increased confinement of tail particles, in roughly Maxwellian particle distributions.

The imposition of a multipole field on top of an axisymmetric configuration breaks the axisymmetry and the associated Noether invariant (Noether 1918). Additionally, it modifies the existing adiabatic invariant of the perpendicular motion, $\mu$. Consideration of the system from a Hamiltonian framework allows us to derive these two adiabatic invariants, while deriving the ponderomotive pseudopotentials. Additionally, the same treatment brings out corrections to the particle trajectory. By nature of these multipole fields, the largest repulsive potential would occur near the outer radius of the device. Radial oscillations in the particle trajectory are also largest at the outer radius, which could cause a radial particle loss instead of axial particle confinement. This paper uses the Lie transformation method to change variables from the particle coordinates to guiding centre coordinates, describing the centre of oscillation, and generates an approximate Hamiltonian which determines the time evolution of these variables. The dynamics in these variables is easily solved for. The transformation can be inverted to determine the particle trajectory from the guiding centre time evolution. This procedure has been used in the plasma physics literature for several applications (Brizard 2022, 2023; Cary & Brizard 2009), in addition to works in celestial mechanics such as (Deprit 1981; Martinusi 2020).

Another effect of note in ponderomotive interactions is modification of the effective particle mass (Dodin & Fisch 2008; Zhmoginov, Dodin & Fisch 2011), when considering the dynamics in the perturbation region. This effect does not affect axial particle confinement, but is an important aspect of the particle dynamics, which was neglected in our previous work.

The purpose of this work is to formally explore the dynamics of the field configuration proposed in (Rubin et al. 2023) up to second order in a perturbed quantity that is related to the ratio of the energy of the multipole field to the energy in the axial magnetic field and radial electric field, in order to find the Miller pseudopotential, mass effects, evaluate the radial excursions and treat the case of resonant periodic motion. The second-order corrections to the potentials are particularly important near resonances, including the resonance corresponding to zero rotation. Small rotation frequencies are practically favourable, with lower energy content in the plasma and lower electric fields touching solid matter. We also suggest a method to utilize this field configuration for mass separation, which is useful for nuclear waste treatment (Litvak et al. 2003; Gueroult & Fisch 2014; Timofeev 2014; Gueroult, Hobbs & Fisch 2015; Vorona et al. 2015; Dolgolenko & Muromkin 2017).

This paper is organized as follows. In § 2, we present the action-angle form for the Hamiltonian of an axially magnetized and rotating plasma column, with a multipole magnetic field added to it, and identify several of its features. In § 3, we find the adiabatic invariants and ponderomotive potentials for reflection off of the multipole field, up to second order in the energy ratio. We find the radial excursions of these particles to that order, and consider axial reflection and radial deconfinement. We give the simple example of a particle with zero gyroradius. In § 4 we consider the case of resonant periodic motion, and investigate the confinement properties of such a configuration. In Appendix A we present the Lie transformation to the guiding centre coordinates.

2 Model

The non-relativistic Hamiltonian of a charged particle with charge $e$ and mass $m$, with a canonical momentum $\mathbf{p}$ interacting with the static fields $\boldsymbol {B} = \boldsymbol {\nabla }\times \boldsymbol {A}$ and $\boldsymbol {E} =-\boldsymbol {\nabla }\varPhi$, derived from the vector potenital $\mathbf{A}$ and the electrostatic potential $\Phi$ is

(2.1)\begin{equation} H(\boldsymbol{x}, \boldsymbol{p}) = \frac{\left(\boldsymbol{p}-e\boldsymbol{A}\right)^2}{2m} + e\varPhi.\end{equation}

We consider particle motion in a magnetized, rotating plasma column, with a multipole magnetic field added on. The axial magnetic field $\boldsymbol {B}_0 = B_z \boldsymbol {e}_z$, and a radial electric field, $\boldsymbol {E}_0 = -r\omega B_z \boldsymbol {e}_r$ describe the magnetization and solid-body rotation. Here, $B_z$, $\omega$ are constants, and $r, \alpha, z$ are polar coordinates, associated with the right-handed basis $( \boldsymbol {e}_{r},\boldsymbol {e}_{\alpha },\boldsymbol {e}_{z})$. A set of Cartesian coordinates $x,y,z$ is defined such that $x=r\cos \alpha$, $y=r\sin \alpha$ and $z=z$. The added multipole field

(2.2)\begin{equation} \boldsymbol{B}_1 = B_w f(z) \left(\frac{r}{R}\right)^{n-1}\left(\sin\left(n\alpha\right)\boldsymbol{e}_r+\cos\left(n\alpha\right)\boldsymbol{e}_\alpha\right),\quad r< R, \end{equation}

is used in order to both break axisymmetry, and to add an oscillating field in the (non-inertial) frame rotating with the plasma. Here, $n$ is an integer, $R$ is the radius of the cylindrical current sheet generating this field, $B_w$ is the amplitude of the perturbation and $f:\mathbb {R}\to [0,1]$ is a shape function, representing the ramp up of the multipole field.

A set of scalar and vector potentials generating these field is given by

(2.3)\begin{gather} \varPhi_0 = \tfrac{1}{2}r^2B_z\omega = \tfrac{1}{2}B\omega(x^2+y^2), \end{gather}

(2.4)\begin{gather}\boldsymbol{A}_0 = \tfrac{1}{2}rB\boldsymbol{e}_\alpha = \tfrac{1}{2}(x\boldsymbol{e}_y-y\boldsymbol{e}_x)B_z, \end{gather}

(2.5)\begin{gather}\boldsymbol{A}_1 =-B_w f(z)\frac{R}{n}\left(\frac{r}{R}\right) ^{n}\cos \left(n\alpha\right) \boldsymbol{e}_{z} , \quad r< R. \end{gather}

2.1 Action-angle variables

Substituting (2.3), (2.4) into (2.1) yields the Hamiltonian in Cartesian coordinates

(2.6)\begin{equation} H_0 = \frac{1}{2m}(\,p_x^2+p_y^2+p_z^2)-\frac{1}{2}\varOmega_c(xp_y-yp_x)+ \frac{1}{8}m(\varOmega_c^2+4\varOmega_c\omega)(x^2+y^2). \end{equation}

Following the work of Brillouin (1945) and Davidson (1990), particle motion in these fields is integrable. The Brillouin frequency $\varOmega _B$ is related to the cyclotron frequency $\varOmega _{c}$ and the $\boldsymbol {E}\times \boldsymbol {B}$ rotation frequency $\omega$ by

(2.7a,b)\begin{equation} \varOmega _{c}=\frac{eB_z}{m}, \quad \varOmega_B =\sqrt{\varOmega _{c}^{2}+4\omega \varOmega_c}. \end{equation}

In order for a particle to remain confined within these fields, rather than be accelerated radially by the electric field, $\varOmega _B$ must remain real, meaning $\omega /\varOmega _c > -1/4$.

Using a canonical transformation, generated by the generating function

(2.8)\begin{equation} F = \frac{1}{8\,m \varOmega_B}(2 p_x-m \varOmega_B y )^2\cot (\theta ) +\frac{1}{8\,m \varOmega_B} (2 p_x+m \varOmega_B y)^2\cot (\varphi), \end{equation}

we find new coordinates; actions $D, J$ and angles $\theta, \varphi$, related to the Cartesian coordinates $x, y$ and their conjugate momenta $p_x, p_y$, by

(2.9)\begin{gather} x =\sqrt{\frac{2}{m\varOmega_B}}\left(\sqrt{D}\cos \theta -\sqrt{J}\cos\varphi\right), \end{gather}

(2.10)\begin{gather}y=\sqrt{\frac{2}{m\varOmega_B}}\left(\sqrt{D}\sin \theta +\sqrt{J}\sin \varphi\right) , \end{gather}

(2.11)\begin{gather}p_{x} =\sqrt{\frac{1}{2}m\varOmega_B}\left(-\sqrt{D}\sin \theta +\sqrt{J}\sin \varphi\right) , \end{gather}

(2.12)\begin{gather}p_{y}=\sqrt{\frac{1}{2}m\varOmega_B}\left(\sqrt{D}\cos \theta +\sqrt{J}\cos \varphi\right). \end{gather}

The axial coordinate $z$ and momentum $p_z$ remain unchanged. The action $D$ is the orbital angular momentum, and the action $J$ is the spin angular momentum. In the fields $\boldsymbol {E}_0$ and $\boldsymbol {B}_0$, they are related to the gyrocentre position $R_G$ and the gyroradius $\rho$ by

(2.13a,b)\begin{equation} D = \tfrac{1}{2}m\varOmega_BR_G^2,\quad J = \tfrac{1}{2}m\varOmega_B\rho^2. \end{equation}

In these coordinates, the Hamiltonian for the particle interaction with $\boldsymbol {E}_0$ and $\boldsymbol {B}_0$ can be expressed as

(2.14)\begin{equation} H_0= \frac{p_z^2}{2m}+\varOmega_-{D} - \varOmega_+ {J}.\end{equation}

This motion is described by two uncoupled harmonic oscillators and a free degree of freedom. The two harmonic oscillators generate a cycloid motion with frequencies

(2.15)\begin{equation} \varOmega _{{\pm} }=-\tfrac{1}{2}(\varOmega _{c}\pm \varOmega_B). \end{equation}

2.2 Fourier expansion of the perturbation

The addition of the multipole field to the integrable Hamiltonian (2.14) consists of the two terms $e^2\boldsymbol {A}_1^2/2m$, and $-(\boldsymbol {p}-e\boldsymbol {A}_0)\boldsymbol {\cdot }e\boldsymbol {A}_1/m$. For these vector potentials $\boldsymbol {A}_0\boldsymbol {\cdot } \boldsymbol {A}_1=0$ and $\boldsymbol {p}\boldsymbol {\cdot } \boldsymbol {A}_1=p_z A_{1z}$. We write

(2.16)\begin{gather} H_1=-\frac{p_ze A_{1z}}{m}+\frac{e^2 A_{1z}^{2}}{2m}, \end{gather}

(2.17)\begin{gather}-\frac{p_ze A_{1z}}{m} = p_z \varOmega_w f(z) \frac{R}{n}\left(\frac{r}{R}\right) ^{n}\cos \left(n\alpha\right), \end{gather}

(2.18)\begin{gather}\frac{e^2 A_{1z}^{2}}{2m} =\frac{1}{4} m \varOmega_w ^2f^2(z)\frac{R^2}{n^2}\left(\frac{r}{R}\right) ^{2n}\left(1+\cos\left(2n\alpha\right)\right). \end{gather}

Where $\varOmega _w = e B_w/m$ is the cyclotron frequency related to the amplitude of the perturbation field.

We now require the axial momentum to be of the same order as the perturbation vector potential

(2.19)\begin{equation} p_z\sim e\boldsymbol{A}_{1z}, \end{equation}

so the inertial term $p_z^2/2m$ would be balanced by $e^2\boldsymbol {A}_1^2/2m$. This brings the two terms in the perturbation to be of the same order.

The term $p_z^2/2m$ still belongs in $H_0$, because neglecting it reduces the dimensionality of the motion. Furthermore, $H_0$ represents the integrable motion, and the motion of a free particle is integrable.

Substituting (2.9), (2.10) in $r^n \cos (n\alpha )$ yields the $r, \alpha$ dependence of $A_{1z}$ in action-angle form

(2.20)\begin{equation} r^n \cos (n\alpha)=\left(\frac{2}{m\varOmega_B}\right)^{n/2}\sum_{\ell=0}^{n}\mathcal{C}^{n}_{\ell} (-1)^\ell{D}^{(n-\ell)/2}{J}^{\ell/2} \cos(\ell(\theta+\varphi)-n\theta). \end{equation}

The $A_{1}^2$ term can be derived from the $\cos ^2(n\alpha )=(1+\cos (2n\alpha ))/2$ relation, or from squaring the sum in (2.20) and using the cosine product identity

(2.21)\begin{equation} P=mv_z+eA_{wz}=mv_z-eB_w f(z)\frac{R}{n}\sum_\ell\mathcal{C}^{n}_{\ell} (-1)^\ell\mathcal{D}^{n/2-\ell/2}\mathcal{J}^{\ell/2} \cos((\theta+\varphi)\ell-n\theta). \end{equation}

Substituting (2.9), (2.10), (2.11) and (2.12) into (2.17) and (2.18) yields the contribution of the $\boldsymbol {A}_1$ field to the energy in action-angle form

(2.22)\begin{align} H_1 & = p_z \frac{ \varOmega_w R }{n}f(z)\sum_{\ell=0}^{n} U_{\ell} \cos(\ell(\theta+\varphi)-n\theta)\nonumber\\ & \quad + \frac{1}{4} m\varOmega_w^2\frac{R^2}{n^2}f^2(z)\left(\sum_{\ell=0}^{n} V_{0\ell} \cos(\ell(\theta+\varphi)) + \sum_{\ell=0}^{2n} V_{2\ell} \cos(\ell(\theta+\varphi)-2n\theta)\right), \end{align}

with the coefficients $U_\ell, V_{0\ell }, V_{2\ell }$, defined in (2.30), (2.31) and (2.32).

The dimensionless variables of this system are achieved by the following scaling:

(2.23a,b)\begin{gather} \mathcal{D} = \frac{2D}{m\varOmega_B R^2}= \frac{R_G^2}{R^2},\quad \mathcal{J} = \frac{2J}{m\varOmega_B R^2}= \frac{\rho^2}{R^2}, \end{gather}

(2.24a)\begin{gather} \zeta = \frac{z}{R},\quad \mathcal{P} = \frac{2 p_z}{m\varOmega_B R},\quad \mathcal{H} = \frac{2H}{m\varOmega_B \varOmega_c R^2},\quad \tau = \varOmega_c t,\quad g(\zeta) = f(R\zeta). \end{gather}

2.3 Dimensionless treatment

The dimensionless Hamiltonian $\mathcal {H} = \mathcal {H}_0+\mathcal {H}_1$ for the motion in these fields can be written compactly as

(2.25)\begin{gather} \mathcal{H}_0= \tfrac{1}{4}\varOmega_b\mathcal{P}^2+\omega_-\mathcal{D} - \omega_+ \mathcal{J}, \end{gather}

(2.26)\begin{gather}\mathcal{H}_1 = \sum_{\sigma,\ell}\mathcal{V}_{\sigma, \ell} \cos \varTheta_{\sigma,\ell}, \end{gather}

with the angular dependence $\varTheta _{\sigma,\ell } = (\ell -\sigma n)\theta +\ell \varphi$, with $\ell, \sigma$ integers. Here, $\mathcal {H}_0$ contains the contributions of the $\boldsymbol {E}_0, \boldsymbol {B}_0$ fields and the inertia and $\mathcal {H}_1$ is the added contribution of $\boldsymbol {B}_1$. The small parameter $\epsilon$ is the ratio of the multipole field energy to the other electromagnetic field energies

(2.27a,b)\begin{equation} \epsilon =\frac{ \varOmega_w^2}{n^2 \varOmega_c^2 \varOmega_b}, \quad \mathcal{H}_1\sim \epsilon \mathcal{H}_0. \end{equation}

The frequencies are defined by

(2.28a,b)\begin{equation} \omega _{{\pm} }=-\frac{1}{2}(1\pm \varOmega_b),\quad \varOmega_b=\sqrt{1+4\frac{\omega}{\varOmega_c}}. \end{equation}

The coefficients $\mathcal {V}_{\sigma, \ell }$ of the perturbation are

(2.29)\begin{equation} \mathcal{V}_{\sigma, \ell}=\begin{cases} \sqrt{\varOmega_b} \mathcal{P} \sqrt{\epsilon} g(\zeta) U_{\ell}, & \sigma =1,\\ \tfrac{1}{2}\epsilon g^2(\zeta)V_{\sigma,\ell}, & \sigma\in\{0,2\},\\ 0 & \mathrm{otherwise}, \end{cases}\end{equation}

with the radial dependence enclosed in $U_\ell, V_{0\ell }, V_{2\ell }$, which depend only on the actions $\mathcal {D}, \mathcal {J}$, defined by

(2.30)\begin{gather} U_{\ell} = (-1)^\ell \mathcal{C}_{\ell}^{n} \mathcal{D}^{n/2-\ell/2} \mathcal{J}^{\ell/2}, \end{gather}

(2.31)\begin{gather}V_{0,\ell} =\begin{cases} \sum_{i=0}^{n/2}\mathcal{C}_{2i}^{n}\mathcal{C}_{i}^{2i}\left(\mathcal{D}+\mathcal{J}\right)^{n-2i}{\left(\mathcal{D}\mathcal{J}\right)^{i}}, & \ell = 0, \\ \sum_{i=1}^{n}\sum_{j=0}^{(i-1)/2} (-1)^{i} 2 \mathcal{C}_{i}^{n}\mathcal{C}_{j}^{i}(\mathcal{D}+\mathcal{J})^{n-i}(\mathcal{D}\mathcal{J})^{i/2} \delta_{\ell, i-2j}, & \ell\neq 0, \end{cases}\end{gather}

(2.32)\begin{gather}V_{2,\ell} = (-1)^\ell\mathcal{C}_{\ell}^{2n}\mathcal{D}^{n-\ell/2} \mathcal{J}^{\ell/2}. \end{gather}

The symbol $\delta _{\ell, i-2j}$ is the Kronecker delta, with indices $\ell$ and $i-2j$, and $\mathcal {C}_{\ell }^{n}= \left( \begin{smallmatrix}{n}\\ {\ell }\end{smallmatrix}\right)=n!/\ell !( n-\ell )!$ are the binomial coefficients, which are defined to be 0 for $\ell <0$ and $\ell >n$.

The sum over $\ell$ in (2.26) is implicitly defined by the binomial coefficients in $V_{\sigma,\ell }, U_\ell$; $\ell \in \{0,\ldots,n\}$ for $\sigma \in \{0,1\}$, and $\ell \in \{0,\ldots,2n\}$ for $\sigma = 2$.

The dynamics of the system is generated by the Hamiltonian using the Poisson bracket. Grouping the actions and angles as $\boldsymbol {P} = (\mathcal {P},\mathcal {D},\mathcal {J})$ and $\boldsymbol {Q} = (\zeta,\theta,\varphi )$, respectively. The Poisson bracket

(2.33)\begin{equation} \{F,G\} = \frac{\partial F}{\partial \zeta}\frac{\partial G}{\partial \mathcal{P}}-\frac{\partial F}{\partial \mathcal{P}}\frac{\partial G}{\partial \zeta}+\frac{\partial F}{\partial \theta}\frac{\partial G}{\partial \mathcal{D}}-\frac{\partial F}{\partial \mathcal{D}}\frac{\partial G}{\partial \theta}+\frac{\partial F}{\partial \varphi}\frac{\partial G}{\partial \mathcal{J}}-\frac{\partial F}{\partial \mathcal{J}}\frac{\partial G}{\partial \varphi}, \end{equation}

if defined for any functions $F(\boldsymbol {P},\boldsymbol {Q}), G(\boldsymbol {P},\boldsymbol {Q})$ of phase space. When substituting the Hamiltonian in place of $G$, this operator gives the time derivative of $F$.

It is clear that $\{\boldsymbol {P},\mathcal {H}_0\}=0$, the actions are invariants under interaction with the unperturbed fields. It is also apparent that this is no longer the case when the multipole field is introduced $\{\boldsymbol {P},\mathcal {H}_1\}\neq 0$. Also, the angles satisfy $\{\boldsymbol {Q},\mathcal {H}\}=\{\boldsymbol {Q},\mathcal {H}_0\}+O(\epsilon ) = (\varOmega _b\mathcal {P}/2, \omega _-,-\omega _+)+O(\epsilon )$.

Looking at $\mathcal {H}_1$, we identify the following;

1. (i) The term $\tfrac {1}{2}\epsilon g^2 V_{00}$ is independent of the angles $\theta$ and $\varphi$. This is the repulsive potential identified in (Rubin et al. 2023), and alone, it would reflect particles with $0<\mathcal {P} < \sqrt {2 \epsilon V_{00}/ \varOmega _b}$ entering into interaction with the multipole field from a region in which $g=0$. The three terms in (2.25) in addition to this fourth term constitute the integrable part of the Hamiltonian. This forth nonlinear term in the actions $D, J$ causes a shift in the frequencies of rotation.

2. (ii) Miller potentials can be derived from the oscillating terms in $\mathcal {H}_1$. The terms derived from $p_z A_{1z}$ are going to be proportional to $\epsilon \mathcal {P}^2$, which are a ponderomotive modification to the effective mass of the particle, when considering the axial dynamics. The terms derived from $A_{1z}^2$ are going to be proportional to $\epsilon ^2$ and may be attractive or repulsive ponderomotive terms, similar to the ones derived from radio frequency (RF) waves, for example in Dodin & Fisch (2005).

3. (iii) For a certain frequency ratio $\omega /\varOmega _c$, the motion can become periodic with the same periodicity as the multipole. For some choices of $\omega /\varOmega _c$, only one term in $\mathcal {H}_1$ becomes near constant, while other choices result in two terms becoming near constant at the same time. We treat periodic motion in § 4.

We proceed to investigate the phase space volume of confined particles. In § 3, we treat the adiabatic case, where the coordinates $\theta, \varphi$ do not affect the ponderomotive potentials, and can be averaged out. However, the motion in the $x$–$y$ plane is constrained to the previously mentioned limit of $\sqrt {\mathcal {D}}+\sqrt {\mathcal {J}}<1$. In general, the oscillating terms in $\mathcal {H}_1$ would introduce $O(\sqrt {\epsilon })$ variations in $\mathcal {D}, \mathcal {J}$. As such, some particles would be pushed into $r>R$ due to the interaction with the multipole field.

3 Ponderomotive potentials away from resonance

Having identified the components of the Hamiltonian, we approach the task of employing canonical perturbation theory in order to asymptotically solve for the particle trajectory. We use the Lie–Poisson method outlined in Deprit (1969) and Cary (1981) to perform a canonical transformation (symplectic, volume preserving) of the phase space variables, from the old action and angle variables $\boldsymbol {P},\boldsymbol {Q}$, which exhibit complex oscillatory behaviour, to a new set of ‘guiding centre’ actions and angles $\bar {\boldsymbol {P}} = (\bar {\mathcal {P}},\bar {\mathcal {D}},\bar {\mathcal {J}})$ and $\bar {\boldsymbol {Q}} = (\bar {\zeta },\bar {\theta },\bar {\varphi })$. The guiding centre variables are selected such that their evolution exhibits smaller oscillations (at the cost of higher oscillation frequency). The Hamiltonian generating their evolution is the guiding centre Hamiltonian $\mathcal {K}$, which we will pick to be independent of $\bar {\theta }$, $\bar {\varphi }$, and we approximate the evolution of the new variables up to $O(\epsilon ^2)$. The variable transformation is detailed in Appendix A. We consider the limits

(3.1)\begin{gather} \epsilon < 1, \end{gather}

(3.2)\begin{gather}\forall (\sigma,\ell)\neq (0,0): \frac{R}{L}\frac{\varOmega_b\mathcal{P}}{2}\frac{1}{\varOmega_{\sigma,\ell}} \ll1. \end{gather}

Equation (3.1) corresponds to a small multipole field compared with the axial and radial electromagnetic fields $\boldsymbol {E}_0, \boldsymbol {B}_0$, and (3.2) requires each term in the perturbed Hamiltonian $\mathcal {H}_1$ to vary smoothly along the particle trajectory. Here, $\varOmega _{\sigma,\ell }=\{\varTheta _{\sigma,\ell },\mathcal {H}_0\}=(\ell -\sigma n)\omega _--\ell \omega _+$, and $R/L = g'(\zeta )$.

The ponderomotive pseudopotentials are derived in Appendix A. The second-order effects appear in the sum of the Miller-type terms in the approximate guiding centre Hamiltonian

(3.3)\begin{equation} \mathcal{K} =\frac{1}{4}\varOmega_b\bar{\mathcal{P}}^2+\omega_-\bar{\mathcal{D}} - \omega_+ \bar{\mathcal{J}}+ \mathcal{V}_{00}-\frac{1}{4}\sum_{(\sigma,\ell)\neq (0,0)}\frac{\boldsymbol{\nabla}_{D,\sigma,\ell}\mathcal{V}_{\sigma,\ell}^2}{\varOmega_{\sigma,\ell}},\end{equation}

where $\boldsymbol {\nabla }_{D,\sigma,\ell } = (\ell -\sigma n)\partial _{\mathcal {D}}+\ell \partial _{\mathcal {J}}$, and all the coefficients $\mathcal {V}$ and their derivatives are evaluated at the new actions and angles $\bar {\boldsymbol {P}}, \bar {\boldsymbol {Q}}$. The derivatives with respect to the actions are written compactly. The first four terms in (3.3) are the integrable part of the motion, and the sum constitute the average contribution of the non-integrable part of the motion.

We identify $\bar {\mathcal {D}}, \bar {\mathcal {J}}$ as the two adiabatic invariants, which are the constants of motion for this approximate Hamiltonian. The axial dynamics is generated by

(3.4)\begin{equation} \mathcal{K}_{\text{axial}} = \tfrac{1}{4}\varOmega_b\bar{\mathcal{P}}^2\left(1- \Delta m^{-1}(\bar{\zeta})\right)+ V(\bar{\zeta}), \end{equation}

where we absorbed the terms independent of $\bar {\zeta }$ and $\bar {\mathcal {P}}$ into $\mathcal {K}_{\text {axial}}=\mathcal {K} -\omega _-\bar {\mathcal {D}} + \omega _+ \bar {\mathcal {J}}$. Particle reflection would occur if $\mathcal {K}_{\text {axial}}< V(\bar {\zeta })$. The modification to the mass is generated by the Miller potentials for $\sigma = 1$, due to the $\mathcal {P}$ dependence of $\mathcal {V}_{1\ell }$. The average potential is the sum of $\mathcal {V}_{00}$ and the Miller potentials for $\sigma =0$ and $\sigma =2$. The mass modification $\tilde {m}$ cannot be larger than one, or else the asymptotic expansion procedure fails.

Explicitly, the mass modification term and the ponderomotive potentials are

(3.5)\begin{gather} \Delta m^{-1}= \epsilon g^2\sum_{\ell=0}^{n}\frac{(\ell- n) \partial_{\mathcal{D}}+\ell \partial_{\mathcal{J}}}{\ell (\omega_--\omega_+)-n \omega_-}U_{\ell}^2, \end{gather}

(3.6)\begin{gather}V= \frac{1}{2}\epsilon g^2V_{00}+ \frac{1}{16}\epsilon^2 g^4\left[\sum_{\ell=1}^{n}\frac{\partial_{\mathcal{D}}+\partial_{\mathcal{J}}}{\omega_+-\omega_-}{V}^2_{0 \ell}+\sum_{\ell=0 }^{2n}\frac{(\ell-2 n) \partial_{\mathcal{D}}+\ell \partial_{\mathcal{J}}}{\ell (\omega_+-\omega_-)+2 n \omega_-}{V}^2_{2 \ell}\right]. \end{gather}

Considering particles entering into interaction with the multipole, from a region where the perturbation is zero, the value of $\bar {\mathcal {D}}, \bar {\mathcal {J}}$ is $\mathcal {D}, \mathcal {J}$ at the start of the motion.

3.1 Approximate solution for the particle motion

The approximate Hamiltonian describes the dynamics of the guiding centre of the particle oscillations. Using the two adiabatic invariants to reduce the dimensionality of the problem, this is a one-dimensional system with $\bar {\mathcal {D}}, \bar {\mathcal {J}}$ being constants of motion.

The axial momentum is solved as a function of position by inverting the Hamiltonian

(3.7)\begin{equation} \bar{\mathcal{P}}(\bar{\zeta})={\pm}\sqrt{\frac{\bar{\mathcal{P}}_0^2-\dfrac{4}{\varOmega_b} V(\bar{\zeta})}{1- \Delta m^{-1}(\bar{\zeta})}}, \end{equation}

with $\bar {\mathcal {P}}_0$ being the axial momentum outside of the multipole region.

The guiding centre is implicitly given by inverting the integral

(3.8)\begin{equation} \tau = \int^{\bar{\zeta}} \frac{\mathrm{d}s}{\sqrt{\varOmega_b(\varOmega_b\bar{\mathcal{P}}_0^2/4- V(s))(1- \Delta m^{-1}(s))}}, \end{equation}

where the integration is performed along the particle path.

Figure 1 presents a comparison between a numerical solution for the motion and some approximate Hamiltonians. The blue curve on the left plot presents the energy in the axial motion, $\varOmega _b \mathcal {P}^2/4$, as calculated numerically. The orange curve is the solution for the motion with $V(\bar {\zeta }) = \mathcal {V}_{0,0}$ and $\Delta m^{-1}=0$, as presented in Rubin et al. (2023). The green curve takes into account the mass shift term in (3.5). The red curve uses the mass shift term and the full potential in (3.6). In this case the red curve is the best fit amongst the three analytic expressions, even though none of them predicts the exact reflection point.

Figure 1. (a) Energy in the axial degree of freedom. In blue: numerical solution to the energy in the axial degree of freedom. Reflection occurs when the axial energy reaches zero (purple line). Orange: approximate solution with the potential being only $\mathcal {V}_{0,0}$. Green: approximate solution taking into account the $\mathcal {V}_{0,0}$ and the mass shift term. Red: approximate solution taking into account all terms in (3.3). (b) Numerical solution of the trajectory, projected on the $x$–$y$ plane. Thin black circles: inner and outer radii of unperturbed cycloid motion. Parameters: $\omega /\varOmega _c = -0.012$, $n=2$, $\epsilon = 0.01$, $\bar {\mathcal {P}}_0=0.068$, $\bar {\mathcal {D}}=0.65$, $\bar {\mathcal {J}}=0.00005$.

The solution for the motion was performed using a second-order volume preserving particle pusher (Zenitani & Umeda 2018), generalizing Boris’ method (Boris 1970; Qin et al. 2013), using the LOOPP code (Ochs & Fisch 2021, 2023b).

The old variables are related to the gyrocentre variables, up to the first order in $\mathcal {V}$, by the expressions

(3.9)\begin{gather} \mathcal{D}= \bar{\mathcal{D}}-\sum_{(\sigma,\ell)\neq (0,0)}\frac{\ell-\sigma n}{\varOmega_{\sigma,\ell}}\mathcal{V}_{\sigma,\ell}(\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}) \cos\bar{\varTheta}_{\sigma,\ell}, \end{gather}

(3.10)\begin{gather}\mathcal{J}= \bar{\mathcal{J}}-\sum_{(\sigma,\ell)\neq (0,0)}\frac{\ell}{\varOmega_{\sigma,\ell}}\mathcal{V}_{\sigma,\ell}(\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}) \cos\bar{\varTheta}_{\sigma,\ell}, \end{gather}

(3.11)\begin{gather}\mathcal{P}= \bar{\mathcal{P}}-\sum_{(\sigma,\ell)\neq (0,0)}\frac{1}{\varOmega_{\sigma,\ell}}\frac{\partial \mathcal{V}_{\sigma,\ell}}{\partial \zeta}(\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}) \cos\bar{\varTheta}_{\sigma,\ell}, \end{gather}

(3.12)\begin{gather}\zeta = \bar{\zeta} + \frac{1}{\varOmega_{1,\ell}}\sqrt{\varOmega_b} \sqrt{\epsilon} g(\bar{\zeta}) \sum_{\ell}U_{\ell}(\bar{\mathcal{D}},\bar{\mathcal{J}}) \sin\bar{\varTheta}_{1,\ell}, \end{gather}

(3.13)\begin{gather}\theta= \bar{\theta}+\sum_{(\sigma,\ell)\neq (0,0)}\frac{1}{\varOmega_{\sigma,\ell}}\frac{\partial \mathcal{V}_{\sigma,\ell}}{\partial \mathcal{D}}(\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}) \sin\bar{\varTheta}_{\sigma,\ell}, \end{gather}

(3.14)\begin{gather}\varphi= \bar{\varphi}+\sum_{(\sigma,\ell)\neq (0,0)}\frac{1}{\varOmega_{\sigma,\ell}}\frac{\partial \mathcal{V}_{\sigma,\ell}}{\partial \mathcal{J}}(\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}) \sin\bar{\varTheta}_{\sigma,\ell},\end{gather}

which are the zeroth- and first-order terms in (A9) and (A10).

The time evolution of the angles is given by Hamilton's equations

(3.15)\begin{gather} \bar{\theta}= \bar{\theta}_0+\tau\left[\omega_-+\partial_{\mathcal{D}}\left(\mathcal{V}_{00}-\frac{1}{4}\sum_{(\sigma,\ell)\neq (0,0)}\frac{\boldsymbol{\nabla}_{D,\sigma,\ell}\mathcal{V}_{\sigma,\ell}^2}{\varOmega_{\sigma,\ell}}\right)\right]_{\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}}, \end{gather}

(3.16)\begin{gather}\bar{\varphi} = \bar{\varphi}_0+\tau\left[-\omega_+ +\partial_{\mathcal{J}}\left(\mathcal{V}_{00}-\frac{1}{4}\sum_{(\sigma,\ell)\neq (0,0)}\frac{\boldsymbol{\nabla}_{D,\sigma,\ell}\mathcal{V}_{\sigma,\ell}^2}{\varOmega_{\sigma,\ell}}\right)\right]_{\bar{\mathcal{D}},\bar{\mathcal{J}},\bar{\zeta}}. \end{gather}

The oscillations around the gyrocentre position are evident in these expressions. Of special note is the axial oscillations in ${\zeta }$, which are of $O(\sqrt {\epsilon })$, whereas the oscillations in ${\mathcal {D}}, {\mathcal {J}}, {\theta }, {\varphi }$ are of $O(\epsilon )$. This degree of freedom stores the most energy, when particles interact with the multipole field. The oscillations in ${\mathcal {P}}$ are of $O(\epsilon ^{3/2})$, which is still $\epsilon$ times smaller than ${\mathcal {P}}$.

A comparison between the first order, second order and a numerical solution for the motion in these fields is presented in figures 2 and 3. In this simulation, we used the simplest ramp-up function

(3.17)\begin{equation} g(\zeta)=\begin{cases} 0, & \zeta<-\dfrac{L}{2R},\\ \dfrac{\zeta R}{L}+\dfrac{1}{2}, & \zeta\in \left[-\dfrac{L}{2R},\dfrac{L}{2R}\right],\\ 1, & \zeta>\dfrac{L}{2R}, \end{cases} \end{equation}

with $L/R = 5000$, in order to satisfy (3.2).

Figure 2. Particle trajectory in the $x$–$y$ plane, near the reflection point. Here, $\omega /\varOmega _c = -0.012$, $n=2$, $\epsilon g = 0.0045$, $\bar {\mathcal {P}}=0$, $\bar {\mathcal {D}}=0.65$, $\bar {\mathcal {J}}=0.00005$. Black: unperturbed trajectory ($\epsilon =0$). Blue: first-order correction. Red: second-order correction. Green: numerical solution. A thin black line shows $r/R=1$.

Figure 3. (a) Particle trajectory in the $z$–$x$ plane. (b) Particle trajectory in the $z$–$y$ plane. Both up to the reflection point. Here, $\omega /\varOmega _c = -0.012$, $n=2$, $\epsilon = 0.01$, $\bar {\mathcal {P}}_0=0.064$, $\bar {\mathcal {D}}=0.65$, $\bar {\mathcal {J}}=0.00005$. Dark line: second-order solution for the envelope of the motion. Green: numerical solution.

Figure 2 shows the particle motion in the $x$–$y$ plane, where the gyrocentre axial momentum is zero, near the reflection point. The unperturbed trajectory is marked in a black line, while the numerical solution is marked in green. The expressions in ((3.9)–(3.14)) generate the blue curve, and the expressions in (A9), (A10), (A5) and (A11) generate the red curve. The red curve approximates the numerical solution fairly well.

Figure 3 shows projections of the particle trajectory on the $z$–$x$ plane, and the $z$–$y$ planes. The motion presented in this figure is the motion starting outside the region of the multipole field, and up to the reflection point. In green is he numerical solution, and the full dark line is $\sqrt {\mathcal {D}}+\sqrt {\mathcal {J}}$, using the second-order expressions for $\mathcal {D}(\bar {\boldsymbol {P}}), \mathcal {J}(\bar {\boldsymbol {P}})$ evaluated at $(\theta =0, \varphi ={\rm \pi} )$ for the $z$–$x$ plot and $(\theta ={\rm \pi} /2, \varphi ={\rm \pi} /2)$ for the $z$–$y$ plot.

3.2 Radial excursions

Having transformed the Hamiltonian to second order in $\mathcal {V}$, we see that careful choice of electric and magnetic fields can cause the phase space volume of confined particles to be extended in $\mathcal {P}$ by the second term of (3.6) over and above the leading-order effect of $\mathcal {V}_{00}$. This allows for particles of higher energy to be reflected. However, the phase space volume is also reduced in $\mathcal {D}, \mathcal {J}$ due to the condition in (3.23), which causes particles that start out at a larger radius to have increased radial excursions and hit the machine wall.

Particles are confined radially if $(r/R)^2= \mathcal {D}+\mathcal {J}-2\sqrt {\mathcal {DJ}}\cos ( \theta +\varphi )<1$. Since both $\mathcal {D}, \mathcal {J}$ are phase dependent, as illustrated in (3.9) and (3.10), one would have to evaluate them at the appropriate angles to determine the radial excursion.

The leading-order expressions contain many terms. A useful limit is the small-gyroradius limit, which reduces the number of terms significantly. We investigate this case in the following subsection.

3.3 Small gyroradius

The zero-gyroradius limit ($\bar {\mathcal {J}}\ll \bar {\mathcal {D}}$) is a useful one, as it reduces the number of non-zero terms in $\mathcal {H}_1$ to the $(\sigma,\ell ) = (0,0),(1,0),(2,0)$ terms. In the gyro-averaged Hamiltonian, the non-zero terms generate the following mass shift and potentials:

(3.18)\begin{gather} \Delta m^{-1}=\epsilon n^2g^2\bar{\mathcal{D}}^{n-1}\left(\frac{1}{\varOmega_{11}}-\frac{1}{\varOmega_{10}}\right), \end{gather}

(3.19)\begin{gather}V= \frac{1}{2}\epsilon g^2\bar{\mathcal{D}}^n +\frac{1}{4}\epsilon^2 n^2g^4\bar{\mathcal{D}}^{2n-1}\left(\frac{1}{\varOmega_{20}}-\frac{1}{\varOmega_{21}}-\frac{1}{\varOmega_{01}}\right). \end{gather}

Since $\epsilon$ is proportional to $n^{-2}$, the mass shift term depends on $n$ only through the frequencies appearing in the brackets of (3.18), and at the power of $\bar {\mathcal {D}}$. At this order, it appears as though the mass term ($\propto (1-\Delta m^{-1})$) can be pushed into the negatives through zero. This would not cause a reflection due to terms at the next order, which are proportional to $\bar {\mathcal {P}}^2$ and $\bar {\mathcal {P}}$. Bringing this term to be of $O(1)$ is a classic example of the problem of small denominators, discussed extensively in Lichtenberg & Lieberman (1983).

The potential, (3.19) determines the reflection boundary. Both terms have a pre-factor of $n^{-2}$ through $\epsilon$ and $\epsilon ^2 n^2$, for constant $\varOmega _w/\varOmega _c$. The denominators in the brackets, however, can be either negative or positive. The sum of inverse frequencies appearing in parenthesis is plotted in figure 4. In the case of $\omega /\varOmega _c<0$, corresponding to a positively charged mirror, there is no opportunity for a resonance crossing by varying $\varOmega _c$, in the manner used in Dodin et al. (2004). However, a rotation of the same sign as the gyromotion, $\omega /\varOmega _c >0$ does allow for resonant interaction.

Figure 4. Evaluation of the term in parenthesis in (3.19). The zero crossing and the resonances are marked by solid and dashed lines.

Figure 4 also shows that, for $\omega /\varOmega _c<0$, the Miller potentials are of larger magnitude for low $n$. Combined with the $n$ dependence of $\epsilon$, low $n$ is preferable for potentials used for particle confinement.

The denominators in the mass term are plotted in figure 5. For $\omega /\varOmega _c<0$, the mass shift is positive, increasing the effective mass for motion along the axis. Between the $\omega /\varOmega _c=0$ resonance and the second resonance, the mass shift is negative, hastening the axial dynamics and making the particle less susceptible to added external forces. The sign changes again past the second resonance.

Figure 5. Evaluation of the term in parenthesis in (3.18). The resonances are marked by dashed lines.

In the case of $\omega /\varOmega _c<0$, $\bar {\mathcal {J}}=0$, $\bar {\mathcal {P}}>0$, the expression for the maximal $\mathcal {D}$ is given by

(3.20)\begin{align} \mathcal{D}& \approx \bar{\mathcal{D}} -\frac{ \sqrt{\epsilon} g}{\omega_-}\sqrt{\varOmega_b}\bar{\mathcal{P}} \bar{\mathcal{D}}^{n/2} -\frac{1}{2}\frac{\epsilon g^2}{\omega_-}\bar{\mathcal{D}}^{n} \nonumber\\ & \quad +\frac{1}{2}n\frac{\epsilon^2 g^4}{\omega_-^2}\bar{\mathcal{D}}^{2n-1} +\frac{1}{4}n\frac{\epsilon g^2}{\omega_-^2}\varOmega_b\bar{\mathcal{P}}^2\bar{\mathcal{D}}^{n-1}+\frac{7}{8}n\frac{\epsilon^{3/2}g^3}{\omega_-^2}\sqrt{\varOmega_b}\bar{\mathcal{P}}\bar{\mathcal{D}}^{3n/2-1}, \end{align}

up to second order in $\epsilon$. All the terms in this expression are positive contributions. Remembering $\bar {\mathcal {P}}\approx \bar {\mathcal {P}}_0\sqrt {1-{4 V(\bar {\zeta })}/{\varOmega _b\bar {\mathcal {P}}_0^2}}$, the axial momentum contribution to the radial excursion is large, taking into account the $n/2$ power of $\bar {\mathcal {D}}$ in the second term.

In order for the radial excursions to remain small, we require

(3.21)\begin{equation} \epsilon g < \left|\frac{\omega}{\varOmega_c}\right|,\end{equation}

along the path, up to the reflection point.

The axial reflection condition is

(3.22)\begin{equation} \frac{4 V(\bar{\zeta})}{\varOmega_b\bar{\mathcal{P}}_0^2}<1,\end{equation}

and the radial excursion limit is

(3.23)\begin{equation} \sqrt{\bar{\mathcal{D}}}+\sqrt{\bar{\mathcal{J}}} -\frac{1}{2}\frac{ \sqrt{\epsilon} g}{\omega_-}\sqrt{\varOmega_b}\bar{\mathcal{P}}_0\left(\bar{\mathcal{D}}^{n/2}-\frac{\epsilon g^2 }{\varOmega_b\bar{\mathcal{P}}_0^2}\bar{\mathcal{D}}^{3n/2}\right) -\frac{1}{4}\frac{\epsilon g^2}{\omega_-}\bar{\mathcal{D}}^{n}<1.\end{equation}

Condition (3.21) effectively sets a lower limit for $\omega$, the column rotation, and conditions (3.22) (3.23) limit the phase space of confined particles.

In figure 6, we present the confined, passing and radially lost particle populations, as a function of initial gyrocentre position and axial momentum. The trapped population (in blue), is characterized by having an initial axial energy that is smaller than the ponderomotive pseudopotential, and also small enough such that the radial component of the Lorentz force (using the multipole field and the axial velocity) would not cause large radial excursions. In the case presented here, the additional confinement due to the Miller potentials does not compensate for the large volume of phase space of particles hitting the wall. Larger rotation would have more particles reflected instead of hitting the wall. In a sense, for small rotation, the radial excursions grow faster than the confining potential.

Figure 6. Phase space of confined, passing and radially lost particles, as a function of initial axial momentum and gyrocentre position. Full-colour background: numerical solution. Blue: trapped (reflected) particles. Green: passing particles.Yellow: radially lost particles. Red curve: ponderomotive potential up to second order. Blue dashed curve: leading-order ponderomotive potential. Green curve: approximate trapped – radial loss boundary, using (3.20). Here, $n=2$, $\epsilon = 0.01$, $\bar {\mathcal {J}}=0.00005$. Panels show (a) $\omega /\varOmega _c = -0.012$, (b) $\omega /\varOmega _c = -0.06$.

3.4 Mass separation

The boundary between reflection and radial loss depends strongly on the frequency ratio $\omega /\varOmega _c$, near $\omega =0$. This allows exploration of this configuration for mass-separation purposes. The ratio $\omega /\varOmega _c \propto Z^2/m$, with $Z$ being the ionization number and $m$ the mass number. Particles with high axial energy hit the wall at all initial radial position, for small enough $\omega /\varOmega _c$, whereas for a more negative frequency ratio particles are either reflected or continue along axially.

4 Resonance

4.1 Resonance structure

The perturbation $\mathcal {H}_1$ consists of $4n+3$ terms, all but one of which are proportional to $\cos \varTheta _{\sigma,\ell }$, such that $(\sigma,\ell )\neq (0,0)$. The leading-order time evolution of $\varTheta _{\sigma,\ell }$

(4.1)\begin{equation} \dot{\varTheta}_{\sigma,\ell}=\ell(\dot{\theta}+\dot{\varphi})-2 n\dot{\theta}\approx \{\varTheta_{\sigma,\ell},\mathcal{H}_0\}=\ell\varOmega_b-\sigma n\omega_-, \end{equation}

depends only on $\omega /\varOmega _c$.

Resonance $\varOmega _{\sigma,\ell }=0$, would occur for the $(\sigma,\ell )=(2,\ell _\star )$, $0\leq \ell _\star < n$ term in the perturbation if

(4.2)\begin{equation} \frac{\omega}{\varOmega_c} = \frac{(2 n-\ell_\star)\ell_\star}{(2 n-2\ell_\star)^2}. \end{equation}

All such $\omega /\varOmega _c>0$, necessitating a negatively charged plasma column.

If such $\ell _\star$ is even, the term $(\sigma,\ell )=(1,\ell _\star /2)$ also becomes resonant.

These cases correspond to a periodic motion of the particle, with either the same periodicity as the multipole field (even $\ell _\star$) or twice the periodicity (odd $\ell _\star$).

In this work we do not allow the time evolution of the electromagnetic fields, so we would not see the back reaction of the particle reflection off the fields, or the RF waves generated by the periodic motion.

4.2 Particle reflection at resonance

In resonance, condition (3.2) is not satisfied for the $(2,\ell _\star )$, and possibly also $(1,\ell _\star /2)$, terms.

In this case, the averaging procedure is not applicable, due to the particle not sampling the entire $2{\rm \pi}$ range of $\varTheta _{\sigma,\ell }$ while entering the perturbation.

The approximate Hamiltonian for motion near an odd resonance is

(4.3)\begin{equation} \mathcal{K}_{\mathrm{axial}} =\frac{1}{4}\varOmega_b\bar{\mathcal{P}}^2+ \mathcal{V}_{0,0}+\mathcal{V}_{2,\ell_\star}\cos \bar{\varTheta}_{2,\ell_\star}-\frac{1}{4}\sum_{\substack{(\sigma,\ell)\neq (0,0)\\(\sigma,\ell)\neq (2,\ell_\star)}}\frac{\boldsymbol{\nabla}_{D,\sigma,\ell}\mathcal{V}_{\sigma,\ell}^2}{\varOmega_{\sigma,\ell}},\end{equation}

and the expressions in (3.9), (3.10), (3.11), (3.13) and (3.14) all have the sums skip the $(2,\ell _\star )$ term.

The angle $\varTheta _{2,\ell _\star }$ is a slow-changing or constant angle compared with the rate of motion on axis, corresponding to the opposite limit of (3.2) for that angel

(4.4)\begin{equation} \frac{R}{L}\frac{\varOmega_b\mathcal{P}}{2}\frac{1}{\varOmega_{2,\ell_\star}} \gg 1. \end{equation}

Because this angle depends on the particle initial conditions before interacting with the perturbation, and we assume a thermal particle distribution away from the perturbation, it is customary to take the random angle approximation. In this case, the potential $V(\bar {\zeta })$ the particle experiences as it comes into the perturbation region is increased or reduced by a phase-dependent term, which at most can be $\pm V_{2,\ell _\star }$.

The implication here can be a leaky end plug, where particles might arrive at the right phase $\varTheta _{2,\ell _\star }$ to experience a greatly reduced potential barrier. Particles that are reflected by the barrier might thermalize by collisions within the bulk of the plasma, and try again.

The approximate Hamiltonian for motion near an even resonance is

(4.5)\begin{equation} \mathcal{K}_{\mathrm{axial}} =\frac{1}{4}\varOmega_b\bar{\mathcal{P}}^2(1-\Delta m^{-1})+\sqrt{\varOmega_b} \mathcal{P} \sqrt{\epsilon} g U_{\ell_\star{/}2}\cos \bar{\varTheta}_{1,\ell_\star{/}2}+V(\bar{\zeta}).\end{equation}

This Hamiltonian has the same phase-dependent term in the effective potential as the one appearing in (4.3), and the mass term is also missing the $\ell =\ell _\star /2$ term in the sum. The new phase-dependent term, $\mathcal {V}_{1,\ell _\star /2}$, written here explicitly, is proportional to $\mathcal {P}$.

Hamilton's equations give the relation between the momentum and velocity in this case

(4.6)\begin{equation} \dot{\bar{\zeta}} =\frac{\varOmega_b}{2} \bar{\mathcal{P}}(1-\Delta m^{-1})+\sqrt{\varOmega_b}\sqrt{\epsilon} g U_{\ell_\star{/}2}\cos \bar{\varTheta}_{1,\ell_\star{/}2}. \end{equation}

The gyro-averaged velocity as a function of position is

(4.7)\begin{equation} \dot{\bar{\zeta}} ={\pm} \sqrt{\varOmega_b}\sqrt{\epsilon g^2 U_{\ell_\star{/}2}^2 \cos^2 \bar{\varTheta}_{1,\ell_\star{/}2}+\left(\bar{\mathcal{P}}_0^2-\frac{4V}{\varOmega_b}\right)(1-\Delta m^{-1})}. \end{equation}

Reflection would occur if the expression in the square root reaches zero. This resonance makes reflection less likely due to the added term which is strictly positive.

5 Conclusions

The magnetostatic end-plugging scheme generated by a ponderomotive quasipotential barrier has been further investigated. The ponderomotive barrier, which to leading order is generated by the interaction of azimuthally rotating particles with an azimuthal multipole magnetic field, is now better understood in the small rotation regime. In the small rotation case, radial and azimuthal oscillations join the axial oscillations, and are responsible for Miller-type potentials, which can take either a positive (repulsive) sign, or a negative (attractive) sign.

Plasma rotation has several qualities which makes it desirable for nuclear fusion applications; shear stabilization of instabilities, centrifugal confinement and our innovative magnetostatic ponderomotive end-plugging scheme, among others. However, a too large rotation can be undesirable. Downsides include energy investment in the rotating ions not immediately of use for the fusion reaction, magnetic pressure being spent on the confinement of ion inertia rather than thermal pressure, plasma deformation into an annulus and rotation-driven instabilities. Furthermore, in lower temperature plasma devices, the rotation is limited by the critical ionization temperature. Thus, investigating the case of small rotation is of practical use, as the aforementioned considerations may make this regime favourable.

In the small rotation regime, the Miller potentials are largest, adding the most to the confining pseudopotential. However, we find the limiting design criterion to be the radial excursion of particles. At large initial radii, where the confining potential is greatest, the radial oscillations take particles into the wall, where otherwise they would have been confined by the ponderomotive pseudopotential. The expression for the radial excursions can also be interpreted as the limit for the minimum electric field needed for this type of magnetostatic ponderomotive end plug to be useful, rather than drive much of the plasma to the wall.

One can also integrate the dynamics backwards, and use this end plug to fuel the plasma, by supplying low energy fuel next to the wall, at a position where travelling into the plasma it would end up in an interior flux surface. The particle axial dynamics is further explored, and a mass modification term is found. This effect, which is formally a first-order one, allows for determination of the particle time evolution more precisely.

The investigation performed in this paper is largely confined to the adiabatic regime, where the ramp up of the multipole field is slow compared with the azimuthal rotation. We considered two cases in which this is not the case, the resonances, where the particle rotation has either twice or once the periodicity of the multipole. We found that resonance is not conducive to confinement, either producing a ‘leaky’ potential barrier or, in addition, adding an always-positive term to the axial velocity.

Thus, this investigation ties up some of the loose ends for particle motion in this electric and magnetic fields configuration. One question remaining unanswered by the present work is the practical generation of these electric and magnetic fields. In addition to the matter of how to set up the perturbing fields, there is the matter of how the plasma responds to the imposed fields. It may be that the electric field, or equivalently the distribution of angular momentum in the plasma, would distribute itself such that $|\boldsymbol {E}\boldsymbol {\cdot } \boldsymbol {B}|$ is minimized in the presence of the multipole. In that case, the analysis of Ochs & Fisch (2023a) would pertain.