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 (Reference Brillouin1945) and Davidson (Reference Davidson1990), 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. Reference Rubin, Rax and Fisch2023), 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 (Reference Dodin and Fisch2005).

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.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. (Reference Rubin, Rax and Fisch2023). 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 Reference Zenitani and Umeda2018), generalizing Boris’ method (Boris Reference Boris1970; Qin et al. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013), using the LOOPP code (Ochs & Fisch Reference Ochs and Fisch2021, Reference Ochs and Fisch2023b).

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 (Reference Lichtenberg and Lieberman1983).

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. (Reference Dodin, Fisch and Rax2004). 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.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.