2.1 The Vlasov–Poisson system

A rigorous analysis of the problem of transverse instability at arbitrary magnetic field strength has been presented in an earlier paper (Hutchinson Reference Hutchinson2018b ) which should be consulted for general mathematical detail. We proceed more simply here by making the early approximation of high magnetic field. A linearized analytic treatment of electrostatic instability of a magnetized electron hole depends on the first-order perturbation to the distribution function $f_{1}$ caused by a potential perturbation $\unicode[STIX]{x1D719}_{1}$ . It is found by integrating Vlasov’s equation along the equilibrium (zeroth-order) helical orbits, which are the equation’s ‘characteristics’. For a collisionless situation $f$ is constant along orbits.

The integration can be expressed as an expansion in harmonics of the cyclotron frequency ( $\unicode[STIX]{x1D6FA}=eB/m_{e}$ ). However, if the magnetic field is strong enough, only the $m=0$ harmonic is important. Physically, this high-field approximation amounts to accounting only for particles’ motion along the magnetic field, and ignoring cross-field motion, taking the Larmor radius (and cross-field drifts) to be negligibly small. Vlasov’s equation is then essentially one-dimensional so we need only consider the parallel velocity distribution, denoted $f(v)$ .

Because the equilibrium is non-uniform in the (trapping) $z$ -direction, uncoupled Fourier representation of the potential variation is possible only for the transverse direction (taken as $y$ without loss of generality) orthogonal to $z$ . The $z$ -dependence of the linearized eigenmode must be expressed in a full-wave manner by writing

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{1}(\boldsymbol{x},t)=\hat{\unicode[STIX]{x1D719}}(z)\exp \text{i}(ky-\unicode[STIX]{x1D714}t).\end{eqnarray}$$

One way to derive the solution of Vlasov’s equation intuitively is to recognize that if a small time-independent potential perturbation is applied, then particles’ (parallel) energy is still conserved as they approach from far past time (and distance). Consequently the perturbation to the distribution function ( $f_{1}$ ) at fixed velocity arises purely as a result of the perturbation to the potential in the form $f_{1}(z)=q_{e}\unicode[STIX]{x1D719}_{1}(z)(\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W)$ . This equation expresses the conservation of distribution function along constant energy orbits and the fact that the potential perturbation causes an orbit at fixed velocity to correspond to an energy different (in the distant past, where the distribution is $f_{0}$ ) by $q_{e}\unicode[STIX]{x1D719}_{1}$ . This component is commonly referred to as the ‘adiabatic’ perturbation.

However when the perturbation is time dependent an additional effect of particle energization occurs. Particles no longer move with constant energy. Instead their energy has an instantaneous rate of increase equal to $q_{e}(\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}/\unicode[STIX]{x2202}t)$ , at every position in the past orbit. The energy change from the distant past can be written ${\mathcal{E}}=\int _{-\infty }^{t}q_{e}(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t)[\unicode[STIX]{x1D719}_{1}(z(\unicode[STIX]{x1D70F}),\unicode[STIX]{x1D70F})]\,\text{d}\unicode[STIX]{x1D70F}$ . The starting $f$ value (still conserved along the perturbed orbit) thus corresponds to energy smaller by ${\mathcal{E}}$ . Therefore the distribution perturbation acquires a second ‘non-adiabatic’ component $-{\mathcal{E}}(\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W)$ that, for harmonic time dependence $\propto \text{e}^{-\text{i}\unicode[STIX]{x1D714}t}$ , gives a total

(2.2) $$\begin{eqnarray}\displaystyle f_{1}=q_{e}\unicode[STIX]{x1D719}_{1}(t)\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}+q_{e}\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6F7}\text{e}^{\text{i}(ky-\unicode[STIX]{x1D714}t)}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}, & & \displaystyle\end{eqnarray}$$

where

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}(z,t)\equiv \int _{-\infty }^{t}\hat{\unicode[STIX]{x1D719}}(z(\unicode[STIX]{x1D70F}))\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F},\end{eqnarray}$$

and $z(\unicode[STIX]{x1D70F})=z(t)+\int _{t}^{\unicode[STIX]{x1D70F}}v_{z}(t^{\prime })\,\text{d}t^{\prime }$ is the position at earlier time $\unicode[STIX]{x1D70F}$ (see Hutchinson Reference Hutchinson2018b , equation (5.6)). For positive imaginary part of $\unicode[STIX]{x1D714}$ ( $\unicode[STIX]{x1D714}_{i}>0$ ) the integral converges. We denote the second term of (2.2) omitting the dependence $\text{e}^{\text{i}(ky-\unicode[STIX]{x1D714}t)}$ , as $\tilde{f}\equiv q_{e}\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6F7}\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W$ , which is the ‘non-adiabatic’ distribution perturbation.

The main formal difficulty is to find the shape of the eigenfunction $\hat{\unicode[STIX]{x1D719}}(z)$ which self-consistently satisfies the perturbed Poisson equation: an integro-differential eigenproblem. For slow time dependence relative to particle transit time, it can be argued on general grounds that the eigenmode consists of a spatial shift (by small distance $\unicode[STIX]{x1D709}$ independent of position) of the equilibrium potential profile ( $\unicode[STIX]{x1D719}_{0}(z)$ ) (see Hutchinson Reference Hutchinson2018b , §3.1) giving

(2.4) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D719}}=-\unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{0}}{\unicode[STIX]{x2202}z}.\end{eqnarray}$$

The frequencies we care about have periods not much longer than the particle transit time, at least for particles with total energy near zero; so this shift form cannot be expected to hold exactly. However, we can obtain a good approximation to the corresponding eigenvalue of our system by expressing it in terms of a ‘Rayleigh quotient’. This mathematical procedure is equivalent to requiring the conservation of total $z$ -momentum under the influence of the assumed shift eigenmode. (See Hutchinson & Zhou (Reference Hutchinson and Zhou2016), Zhou & Hutchinson (Reference Zhou and Hutchinson2017) and Hutchinson (Reference Hutchinson2018b ) §3). This amounts to a ‘kinematic’ treatment of the hole as a composite object.

The $z$ -momentum balance can be derived in an elementary way by applying the zeroth- and first-order Poisson equations $\text{d}^{2}\unicode[STIX]{x1D719}_{0}/\text{d}z^{2}=-\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x1D716}_{0}$ , and $\text{d}^{2}\unicode[STIX]{x1D719}_{1}/\text{d}z^{2}-k^{2}\unicode[STIX]{x1D719}_{1}=-\unicode[STIX]{x1D70C}_{1}/\unicode[STIX]{x1D716}_{0}$ (where $\unicode[STIX]{x1D70C}$ is the charge density) to the integral expression for the first-order hole force $\int \unicode[STIX]{x1D70C}_{0}(-\text{d}\unicode[STIX]{x1D719}_{1}/\text{d}z)\,\text{d}z$ , using judicious integrations by parts as follows:

(2.5) $$\begin{eqnarray}\displaystyle -\int \unicode[STIX]{x1D70C}_{0}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\,\text{d}z & = & \displaystyle \unicode[STIX]{x1D716}_{0}\int \frac{\text{d}^{2}\unicode[STIX]{x1D719}_{0}}{\text{d}z^{2}}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\,\text{d}z=-\unicode[STIX]{x1D716}_{0}\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\frac{\text{d}^{2}\unicode[STIX]{x1D719}_{1}}{\text{d}z^{2}}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}(\unicode[STIX]{x1D70C}_{1}-\unicode[STIX]{x1D716}_{0}k^{2}\unicode[STIX]{x1D719}_{1})\,\text{d}z.\end{eqnarray}$$

Now, since $\unicode[STIX]{x1D70C}_{0}$ is a function of $\unicode[STIX]{x1D719}_{0}$ , we can also integrate by parts the other way to find

(2.6) $$\begin{eqnarray}-\int \unicode[STIX]{x1D70C}_{0}\frac{\text{d}\unicode[STIX]{x1D719}_{1}}{\text{d}z}\,\text{d}z=\int \frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}\unicode[STIX]{x1D719}_{0}}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\unicode[STIX]{x1D719}_{1}\,\text{d}z.\end{eqnarray}$$

Combining and rearranging these expressions we get

(2.7) $$\begin{eqnarray}\displaystyle F_{E}\equiv -\unicode[STIX]{x1D716}_{0}k^{2}\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\unicode[STIX]{x1D719}_{1}\text{d}z & = & \displaystyle -\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\left(\unicode[STIX]{x1D70C}_{1}-\frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}\unicode[STIX]{x1D719}_{0}}\unicode[STIX]{x1D719}_{1}\right)\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle -\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\left(\int q_{e}\tilde{f}\,\text{d}v_{z}\right)\,\text{d}z\equiv \tilde{F}.\end{eqnarray}$$

The force $F_{E}$ consists of transfer by Maxwell stress in the $y$ -direction of $z$ -momentum ( $\text{d}/\text{d}y(E_{y}E_{x})$ ). It acts in a direction so as to increase the kink amplitude, in a manner analogous to a compressed spring. The force $\tilde{F}$ is exerted by the equilibrium potential on the non-adiabatic part of the charge density perturbation which is the jetting. Its real part is proportional to kink acceleration, and acts like a negative inertia. The eigenvalue equation (2.7) is that they must balance; we substitute into it the shift-mode form for $\unicode[STIX]{x1D719}_{1}$ , equations (2.4) and (2.1). To lowest order, satisfying the real part of the force equation, the result is an oscillation at a real frequency whose square is proportional to the ratio of the negative tension effect and the negative inertia.

Both the real and imaginary parts of the complex momentum balance equation $F_{E}=\tilde{F}$ must be zero. But $F_{E}$ is real and positive for real $k$ , and in this high-B approximation $k$ appears only in $F_{E}$ and not in $\tilde{F}$ . If we regard $k$ as a free choice, then provided the sign of $\text{Re}(\tilde{F})$ is positive, we can always satisfy $\text{Re}(\tilde{F}-F_{E})=0$ by simply choosing the appropriate value for $k$ . Consequently it is only the imaginary part $\text{Im}(\tilde{F})$ (which is independent of $k$ ) that determines whether there exists a solution of the dispersion relation with a frequency $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ in the upper half-plane ( $\unicode[STIX]{x1D714}_{i}>0$ ), implying instability. We denote contributions from trapped (negative energy, $W<0$ ) and passing ( $W>0$ ) particles with subscripts ‘ $t$ ’ and ‘ $p$ ’ respectively, and write $\tilde{F}=\tilde{F}_{t}+\tilde{F}_{p}$ . At frequencies low compared with the transit time of thermal electrons across the hole, $\text{Re}(\tilde{F}_{t})>0$ , $\text{Re}(\tilde{F}_{p})<0$ , and $|\text{Re}(\tilde{F}_{t})|>|\text{Re}(\tilde{F}_{p})|$ so $\text{Re}(\tilde{F})$ is indeed positive.

2.2 Numerical evaluation observations

A numerical implementation of the required integrations to find $\tilde{F}$ has previously been developed (Hutchinson Reference Hutchinson2018b ) for the specific hole equilibrium

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{0}(z)=\unicode[STIX]{x1D713}\,\text{sech}^{4}(z/4),\end{eqnarray}$$

where the constant $\unicode[STIX]{x1D713}$ is the maximum hole potential: the ‘depth’ of the hole. This code actually performs the integrations for arbitrary magnetic field strength but works well for the present high-B case, with some modifications, described later, newly implemented to allow accurate evaluations at $\unicode[STIX]{x1D714}_{i}\rightarrow 0$ . It shows there to be unstable solutions at low frequency as illustrated by the contours of $\tilde{F}$ plotted in figure 2. An unstable mode occurs at the intersection of the zero contours of the real and imaginary parts of $\tilde{F}$ . Notice that the real part of the frequency is small even for the deep hole $\unicode[STIX]{x1D713}=0.64$ , but the imaginary part is far smaller $\unicode[STIX]{x1D714}_{i}\lesssim 0.02\unicode[STIX]{x1D713}^{0.75}\unicode[STIX]{x1D714}_{r}$ . The shape of the contours is approximately similar for different hole depths, when the frequencies are scaled to $\unicode[STIX]{x1D714}_{r}/\unicode[STIX]{x1D713}^{0.75}$ , $\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D713}^{1.5}$ and the wavenumber to $k/\unicode[STIX]{x1D713}^{0.25}$ . The vertical contours of real part show that there is negligible influence of $\unicode[STIX]{x1D714}_{i}$ on $\text{Re}(\tilde{F})$ in this low-frequency region. The plots are for a specific finite field $\unicode[STIX]{x1D6FA}=10$ , but are essentially unchanged for any $\unicode[STIX]{x1D6FA}\gtrsim 5$ , indicating that we are well into the high-B, one-dimensional motion, regime. We wish to derive analytically the shape of these contours in order to identify the controlling physics of this instability.

Figure 2. Contours of real and imaginary parts of $\tilde{F}-F_{E}$ . Instability occurs for different $k$ at different places along the zero contour $\text{Im}(\tilde{F}-F_{E})=0$ . For the particular $k$ chosen, the eigenfrequency is show as a point.

We concentrate on the decisive imaginary part of

(2.9) $$\begin{eqnarray}\displaystyle \tilde{F}=-(\text{i}\unicode[STIX]{x1D714})\unicode[STIX]{x1D709}\int q_{e}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}q_{e}\int _{-\infty }^{t}\hat{\unicode[STIX]{x1D719}}(z(\unicode[STIX]{x1D70F}))\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}\,\text{d}v\,\text{d}z. & & \displaystyle\end{eqnarray}$$

The resulting lowest order in $\unicode[STIX]{x1D714}$ contribution to $\tilde{F}$ is proportional to $\unicode[STIX]{x1D714}^{2}$ times a real quantity (Hutchinson Reference Hutchinson2018b ). It therefore gives rise to an imaginary component $\text{Im}(\tilde{F})=2(\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D714}_{r})\text{Re}(\tilde{F})$ . The total contribution of this type includes both trapped and passing terms but the trapped real force is typically about five times larger than the passing real force, and will be our focus in the analytic approximation. However, figure 2 shows there are non-zero imaginary components at $\unicode[STIX]{x1D714}_{i}\rightarrow 0$ . We therefore seek, in addition, the imaginary component of the trapped particle force $\tilde{F}_{t}$ of lowest order in $\unicode[STIX]{x1D714}$ that does not depend on $\unicode[STIX]{x1D714}_{i}$ . It comes from accounting to higher order for the variation of the $\text{e}^{-i\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}$ factor, and is contributed by the resonance between the eigenfrequency $\unicode[STIX]{x1D714}$ and the bounce frequency $\unicode[STIX]{x1D714}_{b}$ of some trapped particles. There is also an important imaginary component of $\tilde{F}_{p}$ that does not depend on $\unicode[STIX]{x1D714}_{i}$ .

Figure 3. Example of contributions to the trapped particle force from different particle energies, showing the resonance between bounce frequency and eigenfrequency. Components: solid line, real; dashed line, imaginary; dotted line, $\unicode[STIX]{x1D714}_{r}/2\unicode[STIX]{x1D714}_{i}\times \text{imaginary}$ .

One can exchange the order of integration in (2.9) so as to perform the velocity integration last. The notation $\text{d}\tilde{F}_{t}/\text{d}u$ is used to denote the quantity that when integrated $\text{d}u$ over any velocity-dependent variable $u$ gives $\tilde{F}_{t}$ . Figure 3 shows an example of the real and imaginary parts of $\text{d}\tilde{F}_{t}/\text{d}(-W)^{1/2}$ . The area under the curves gives the total force. The parameter $(-W)^{1/2}$ runs from zero for the highest energy (marginally) trapped particles to $\sqrt{\unicode[STIX]{x1D713}}$ for the deepest trapped particles. As will be shown shortly, it is approximately proportional to the bounce frequency, because of the anharmonic shape of the electron hole potential $\unicode[STIX]{x1D719}_{0}$ , which falls exponentially to zero in the wings. In figure 3 we observe for this low-frequency ( $\unicode[STIX]{x1D714}_{r}=0.02$ , $\unicode[STIX]{x1D714}_{i}=2\times 10^{-4}$ ) case a resonant response at a low value of $(-W)^{1/2}$ (weakly trapped particles) corresponding to $\unicode[STIX]{x1D714}_{r}\simeq \unicode[STIX]{x1D714}_{b}$ . It is obvious that in addition to the substantial real force, arising from the area under the solid curve, there is a smaller imaginary force that is dominated by the contribution from the resonance. Resonances also occur at odd harmonics of the fundamental bounce frequency of particles having lower $\unicode[STIX]{x1D714}_{b}$ that are closer to zero energy (the trapped/passing boundary). However their contribution is much smaller and can be ignored.

The low-level non-resonant imaginary contribution that extends across the entire $(-W)^{1/2}$ range in figure 3 arises from the component $\text{Im}(\tilde{F})=2(\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D714}_{r})\text{Re}(\tilde{F})$ (as is verified by the dotted line, which equals $\text{Re}(\tilde{F})$ far from resonance). It tends to zero as $\unicode[STIX]{x1D714}_{i}\rightarrow 0$ , as expected from analysis; but the narrow resonant contribution does not. Instead it becomes narrower and higher with a convergent total area. This fact poses a numerical challenge at small $\unicode[STIX]{x1D714}_{i}$ for the code that calculates $\tilde{F}$ by discrete integration. The challenge is overcome by adopting the approach made familiar in plasma physics in the context of Landau damping: namely displacing the contour of integration in the $\unicode[STIX]{x1D714}_{b}\leftrightarrow v$ plane so that it remains below the pole at $\unicode[STIX]{x1D714}_{b}=\unicode[STIX]{x1D714}$ as $\unicode[STIX]{x1D714}_{i}$ decreases toward (or even beyond) zero. Keeping the integration contour sufficiently below the pole limits the narrowness and height of the resonance, allowing it to be numerically resolved. We now obtain approximate analytic expressions for the important contributions to the imaginary force at low frequency.

3.1 Bounce orbit integration

First we obtain the relationship between the bounce frequency of trapped particles, $\unicode[STIX]{x1D714}_{b}$ , and (unperturbed) orbit energy $W$ . To approximate analytically for low frequency we observe that orbits near the separatrix ( $W\rightarrow 0$ ) have low $\unicode[STIX]{x1D714}_{b}$ because they dwell most of their time near the turning points. The potential and its gradient are small there, and for the $\,\text{sech}^{4}(z/4)$ hole they can be approximated as $-(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}|z|)\simeq \unicode[STIX]{x1D719}_{0}\simeq \unicode[STIX]{x1D713}16\text{e}^{-|z|}$ . Actually, the mechanism of Debye shielding requires, in the far wing, that $\unicode[STIX]{x1D719}_{0}\propto \text{e}^{-|z|}$ for any hole shape that does not have infinite gradients of $f$ at the separatrix (Hutchinson Reference Hutchinson2017). Therefore the treatment has much wider application than the specific $\,\text{sech}^{4}$ hole shape. In the wing (where $\unicode[STIX]{x1D719}_{0}\sim -W$ ) we may thus take $\text{d}\unicode[STIX]{x1D719}_{0}/\unicode[STIX]{x1D719}_{0}\simeq -\text{d}|z|$ , and can write the relationship between time $t$ to pass from the turning point to a smaller $|z|$ , corresponding to a larger $\unicode[STIX]{x1D719}_{0}$ , as

(3.1) $$\begin{eqnarray}t(z)=\int \frac{-\text{d}|z|}{|v|}\simeq \frac{1}{\sqrt{2}}\int \frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\unicode[STIX]{x1D719}_{0}\sqrt{\unicode[STIX]{x1D719}_{0}+W}}=\frac{1}{\sqrt{2}}\frac{2}{\sqrt{-W}}\tan ^{-1}\sqrt{\frac{\unicode[STIX]{x1D719}_{0}}{(-W)}-1},\end{eqnarray}$$

whose inverse is

(3.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D719}_{0}}{-W}=1+\tan \left(t\sqrt{-W/2}\right).\end{eqnarray}$$

A quarter period of the orbit has been reached when $-\unicode[STIX]{x1D719}_{0}/W$ is large, which occurs when the value of $t\sqrt{-W/2}$ is $\simeq \unicode[STIX]{x03C0}/2$ . Consequently the orbit period is approximately $t_{b}=2\unicode[STIX]{x03C0}/(\sqrt{-W/2})$ , and $\unicode[STIX]{x1D714}_{b}=\sqrt{-W/2}$ .

We shall shortly also need the following average over the orbit, which can be evaluated for the first quarter period using the relation (3.2), initially neglecting the distinction between $\sqrt{-W/2}$ and $\unicode[STIX]{x1D714}_{b}$ ,

(3.3) $$\begin{eqnarray}\displaystyle \left\langle \frac{\unicode[STIX]{x1D719}_{0}}{-W}\cos \unicode[STIX]{x1D714}_{b}t\right\rangle \simeq \frac{1}{t}\int _{0}^{t}\cos \unicode[STIX]{x1D714}_{b}t^{\prime }+\sin \unicode[STIX]{x1D714}_{b}t^{\prime }\,\text{d}t^{\prime }=[\sin \unicode[STIX]{x1D714}_{b}t-\cos \unicode[STIX]{x1D714}_{b}t+1]/\unicode[STIX]{x1D714}_{b}t. & & \displaystyle\end{eqnarray}$$

Taking $\unicode[STIX]{x1D714}_{b}t=\unicode[STIX]{x03C0}/2$ , we get an average equal, by symmetry, to that for a full orbit,

(3.4) $$\begin{eqnarray}\left\langle \unicode[STIX]{x1D719}_{0}\sin \unicode[STIX]{x1D714}_{b}\unicode[STIX]{x1D70F}^{\prime }\right\rangle =\left\langle \unicode[STIX]{x1D719}_{0}\cos \unicode[STIX]{x1D714}_{b}t\right\rangle \simeq \frac{4}{\unicode[STIX]{x03C0}}(-W)\simeq \frac{8}{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D714}_{b}^{2};\end{eqnarray}$$

( $\unicode[STIX]{x1D70F}^{\prime }$ is a shifted time, measured from the centre of the orbit: $\unicode[STIX]{x1D714}_{b}\unicode[STIX]{x1D70F}^{\prime }=\unicode[STIX]{x03C0}/2+\unicode[STIX]{x1D714}_{b}t$ ).

Although these results are exact as $W\rightarrow 0$ , and over-estimate $\unicode[STIX]{x1D714}_{b}$ by only a factor $\sqrt{2}$ at the other extreme $-W\rightarrow \unicode[STIX]{x1D713}$ , their inaccuracy will turn out to be numerically significant. It arises because the estimate of $t_{b}$ has neglected the extra time it takes the orbit to pass from the place where the relationship $\text{d}\unicode[STIX]{x1D719}_{0}/\unicode[STIX]{x1D719}_{0}\simeq -\text{d}|z|$ is broken by the rounded potential peak, to the exact centre of the hole $z=0$ . We track the correction by writing $\unicode[STIX]{x1D714}_{b}=A\sqrt{-W/2}$ where $A$ is a correcting factor close to but slightly below unity. Approximately the same factor applies to the mapping of $W$ to $\unicode[STIX]{x1D714}_{b}^{2}$ in (3.4) which will be written $\langle \unicode[STIX]{x1D719}_{0}\sin \unicode[STIX]{x1D714}_{b}\unicode[STIX]{x1D70F}^{\prime }\rangle \simeq (8/\unicode[STIX]{x03C0})\unicode[STIX]{x1D714}_{b}^{2}/A^{2}$ .

3.2 Resonant force

To evaluate the resonant contribution of trapped particles to the force, it is convenient to integrate (2.3) twice by parts using the fact that $\hat{\unicode[STIX]{x1D719}}=-\unicode[STIX]{x1D709}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)=\unicode[STIX]{x1D709}(\text{d}v/\text{d}\unicode[STIX]{x1D70F})/q_{e}$ ,

(3.5) $$\begin{eqnarray}\frac{q_{e}\unicode[STIX]{x1D6F7}(t)}{\unicode[STIX]{x1D709}}=\int _{-\infty }^{t}\frac{\text{d}v}{\text{d}\unicode[STIX]{x1D70F}}\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=v(t)+\text{i}\unicode[STIX]{x1D714}z(t)+(\text{i}\unicode[STIX]{x1D714})^{2}\int _{-\infty }^{t}z(\unicode[STIX]{x1D70F})\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}.\end{eqnarray}$$

This is exact for an exact shift mode. If the eigenmode deviates from a shift mode (for example by having a shift $\unicode[STIX]{x1D709}$ that is a non-uniform function of $|z|$ ) but retains the properties of the shift mode in that it is antisymmetric and localized to the hole, then the treatment remains valid provided we reinterpret in $\unicode[STIX]{x1D6F7}$ the meaning of $v$ to be $\int q_{e}\hat{\unicode[STIX]{x1D719}}\,\text{d}\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D709}$ and $z=\int v\,\text{d}\unicode[STIX]{x1D70F}$ .

The benefit of this re-expression is that, measuring $\unicode[STIX]{x1D70F}$ from where $z=0$ , $z(\unicode[STIX]{x1D70F})$ is an approximate square wave in time on low bounce frequency orbits. Also, when we integrate the resulting $\tilde{f}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)$ over the hole extent, to obtain $\tilde{F}_{t}$ , the antisymmetry of $(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)$ annihilates any symmetric part of $\tilde{f}$ and hence of $\unicode[STIX]{x1D6F7}$ . Therefore the term $v(t)$ can be dropped. Choose the zero of $\unicode[STIX]{x1D70F}$ and $t$ to be where $v$ is positive, during the bounce period chosen to minimize $|t|$ . Then the integral in (3.5) can be separated into two parts: $\int _{0}^{t}$ with $|t|\leqslant t_{b}/2$ plus $\int _{-\infty }^{0}$ . Only $\int _{-\infty }^{0}$ contributes to the resonant term. The other (non-resonant) terms will be dealt with later.

Represent the square wave during the final orbit as $z(t)=\text{sign}(t)z_{A}$ . The resonant part of the integral may then be evaluated as

(3.6) $$\begin{eqnarray}\displaystyle \frac{q_{e}\unicode[STIX]{x1D6F7}_{R}}{(\text{i}\unicode[STIX]{x1D714})^{2}\unicode[STIX]{x1D709}} & = & \displaystyle \int _{-\infty }^{0}z(\unicode[STIX]{x1D70F})\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=\frac{1}{1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}}}\int _{-t_{b}}^{0}z(\unicode[STIX]{x1D70F})\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}\nonumber\\ \displaystyle & = & \displaystyle \frac{z_{A}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}}{1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}}}\frac{1}{\text{i}\unicode[STIX]{x1D714}}([\text{e}^{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}}]_{-t_{b}/2}^{0}-[\text{e}^{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}}]_{-t_{b}}^{-t_{b}/2})\nonumber\\ \displaystyle & = & \displaystyle \frac{z_{A}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}}{\text{i}\unicode[STIX]{x1D714}}\frac{(1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}/2})^{2}}{1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}}}=\frac{z_{A}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}}{\text{i}\unicode[STIX]{x1D714}}\frac{1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}/2}}{1+\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}/2}}\left[=\frac{z_{A}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}}{\unicode[STIX]{x1D714}}\tan (\unicode[STIX]{x1D714}t_{b}/4)\right].\end{eqnarray}$$

The denominator $1+\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}/2}$ becomes zero, giving resonances, when $\unicode[STIX]{x1D714}t_{b}=2\unicode[STIX]{x03C0}\ell$ with $\ell$ an odd integer. In the vicinity of the $\ell$ th (odd) resonance

(3.7) $$\begin{eqnarray}\frac{1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}/2}}{1+\text{e}^{\text{i}\unicode[STIX]{x1D714}t_{b}/2}}\simeq \frac{2}{1-\text{e}^{\text{i}\unicode[STIX]{x03C0}(\unicode[STIX]{x1D714}-\ell \unicode[STIX]{x1D714}_{b})/\unicode[STIX]{x1D714}_{b}}}\simeq \frac{-2\unicode[STIX]{x1D714}_{b}}{\text{i}\unicode[STIX]{x03C0}(\unicode[STIX]{x1D714}-\ell \unicode[STIX]{x1D714}_{b})}.\end{eqnarray}$$

Thus at resonance

(3.8) $$\begin{eqnarray}\tilde{f}_{R}=q_{e}\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6F7}_{R}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\simeq -\unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}z_{A}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}\frac{(\text{i}\unicode[STIX]{x1D714})^{2}2\unicode[STIX]{x1D714}_{b}}{\text{i}\unicode[STIX]{x03C0}(\unicode[STIX]{x1D714}-\ell \unicode[STIX]{x1D714}_{b})},\end{eqnarray}$$

in which $t$ represents the endpoint $z(t)$ (zero when $t=0$ ) of the orbit: the position where $\tilde{f}$ is being evaluated. Also, in general, for small positive imaginary part $\unicode[STIX]{x1D714}_{i}$ , the integral through (strictly close under) a resonance is

(3.9) $$\begin{eqnarray}\int _{-}^{+}\frac{g(\unicode[STIX]{x1D714}_{b})}{(\ell \unicode[STIX]{x1D714}_{b}-\unicode[STIX]{x1D714}_{r}-\text{i}\unicode[STIX]{x1D714}_{i})}\,\text{d}\unicode[STIX]{x1D714}_{b}\simeq \frac{\text{i}\unicode[STIX]{x03C0}}{\ell }g(\unicode[STIX]{x1D714}_{r}/\ell )\end{eqnarray}$$

for any slowly varying function $g$ .

We now need to get the total non-adiabatic force by integrating $-q_{e}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)\tilde{f}_{R}$ over the relevant phase space $\text{d}z\,\text{d}v$ . Our current interest is the trapped particles, which occupy a finite area of phase space bounded by the separatrix. The best way to carry out this area integral is not along fixed values of $v$ or $z$ , but along fixed values of energy, following the orbits along which particles move as a function of time (denoted $\unicode[STIX]{x1D70F}^{\prime }$ ). It is convenient to represent the orbit energy by the value of the velocity at the hole centre $z=0$ where the potential has its single maximum, $\unicode[STIX]{x1D713}$ , because every orbit does in fact pass through this position. We write this velocity $v_{\unicode[STIX]{x1D713}}$ (taken positive) so that (the negative quantity) $W=q_{e}\unicode[STIX]{x1D713}+\frac{1}{2}v_{\unicode[STIX]{x1D713}}^{2}=q_{e}\unicode[STIX]{x1D719}(z)+\frac{1}{2}v(z)^{2}$ . Then, since at constant energy $v\,\text{d}v=v_{\unicode[STIX]{x1D713}}\,\text{d}v_{\unicode[STIX]{x1D713}}$ , we have $\text{d}z\,\text{d}v\rightarrow v\,\text{d}\unicode[STIX]{x1D70F}^{\prime }\,\text{d}v=\text{d}\unicode[STIX]{x1D70F}^{\prime }v_{\unicode[STIX]{x1D713}}\,\text{d}v_{\unicode[STIX]{x1D713}}=\text{d}\unicode[STIX]{x1D70F}^{\prime }\,\text{d}W\simeq -\text{d}\unicode[STIX]{x1D70F}^{\prime }4\unicode[STIX]{x1D714}_{b}\,\text{d}\unicode[STIX]{x1D714}_{b}/A^{2}$ (using $\unicode[STIX]{x1D714}_{b}^{2}=-A^{2}W/2$ ). Incidentally, the numerical integration code is implemented in the same way, but it uses none of the present approximations.

The resonant force as $\unicode[STIX]{x1D714}_{i}\rightarrow 0$ can then be written, using first the parity relations, and then (3.8), (3.9) and (3.4), as

(3.10) $$\begin{eqnarray}\displaystyle \tilde{F}_{R} & = & \displaystyle -q_{e}\int _{0}^{\sqrt{-2q_{e}\unicode[STIX]{x1D713}}}\int _{-t_{b}/2}^{t_{b}/2}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\tilde{f}_{R}(\unicode[STIX]{x1D70F}^{\prime })\text{d}\unicode[STIX]{x1D70F}^{\prime }v_{\unicode[STIX]{x1D713}}\,\text{d}v_{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & \simeq & \displaystyle -q_{e}\int _{0}^{\sqrt{-2q_{e}\unicode[STIX]{x1D713}}}\int _{0}^{t_{b}/2}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}|z|}[\tilde{f}_{R}(\unicode[STIX]{x1D70F}^{\prime })-\tilde{f}_{R}(-\unicode[STIX]{x1D70F}^{\prime })]\,\text{d}\unicode[STIX]{x1D70F}^{\prime }v_{\unicode[STIX]{x1D713}}\,\text{d}v_{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & \simeq & \displaystyle -q_{e}\unicode[STIX]{x1D709}\int _{\sqrt{\unicode[STIX]{x1D713}}/2}^{0}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\frac{(\text{i}\unicode[STIX]{x1D714})^{2}z_{A}2\unicode[STIX]{x1D714}_{b}}{\text{i}\unicode[STIX]{x03C0}(\unicode[STIX]{x1D714}-\ell \unicode[STIX]{x1D714}_{b})}\int _{0}^{t_{b}/2}\unicode[STIX]{x1D719}_{0}2\text{i}\sin (\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime })\,\text{d}\unicode[STIX]{x1D70F}^{\prime }4\unicode[STIX]{x1D714}_{b}\,\text{d}\unicode[STIX]{x1D714}_{b}/A^{2}\nonumber\\ \displaystyle & \simeq & \displaystyle q_{e}\unicode[STIX]{x1D709}\int _{0}^{\sqrt{\unicode[STIX]{x1D713}}/2}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\frac{(\text{i}\unicode[STIX]{x1D714})^{2}z_{A}2\unicode[STIX]{x1D714}_{b}}{\text{i}\unicode[STIX]{x03C0}(\unicode[STIX]{x1D714}-\ell \unicode[STIX]{x1D714}_{b})}\text{i}\left\langle \unicode[STIX]{x1D719}_{0}\sin (\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime })\right\rangle t_{b}4\unicode[STIX]{x1D714}_{b}\,\text{d}\unicode[STIX]{x1D714}_{b}/A^{2}\nonumber\\ \displaystyle & \simeq & \displaystyle q_{e}\unicode[STIX]{x1D709}\left.\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\right|_{R}(\text{i}\unicode[STIX]{x1D714}_{r})^{2}z_{A}\frac{2\unicode[STIX]{x1D714}_{R}}{\ell }\text{i}\left\langle \unicode[STIX]{x1D719}_{0}\cos (\unicode[STIX]{x1D714}_{r}t)\right\rangle 8\unicode[STIX]{x03C0}/A^{2}\nonumber\\ \displaystyle & \simeq & \displaystyle -\text{i}q_{e}\unicode[STIX]{x1D709}z_{A}\left.\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\right|_{R}128\unicode[STIX]{x1D714}_{R}^{5}/A^{4},\quad \text{for}~\ell =1,\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{R}=\unicode[STIX]{x1D714}_{r}/\ell$ is the resonant bounce frequency. This surprisingly simple expression has not to the author’s knowledge been previously discovered.

It should be remarked that (3.10) has no explicit dependence on the hole depth $\unicode[STIX]{x1D713}$ . However, the (positive) value of $\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W$ varies generally as ${\sim}\unicode[STIX]{x1D713}^{-1/2}$ (Hutchinson, Haakonsen & Zhou Reference Hutchinson, Haakonsen and Zhou2015; Hutchinson Reference Hutchinson2017); more specifically, for a $\,\text{sech}^{4}(z/4)$ hole it can be shown analytically that

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle f_{0}=f_{0s}\left[\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\sqrt{-W}+\frac{15}{16}\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D713}}}W+\exp (-W)\text{erfc}(\sqrt{-W})\right], & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}=f_{0s}\left[\frac{15}{16}\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D713}}}-\exp (-W)\text{erfc}(\sqrt{-W})\right], & \displaystyle\end{eqnarray}$$

where $f_{0s}$ is the distribution function at the separatrix ( $W=0$ , $\exp (0)\text{erfc}(0)=1$ ), and for unit density $f_{0s}=1/\sqrt{2\unicode[STIX]{x03C0}}$ . We will write this $\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W\simeq D(15/16)/\sqrt{2\unicode[STIX]{x1D713}}$ , where $D$ is a correction factor approximately equal to $1-(16/15)\sqrt{\unicode[STIX]{x1D713}/\unicode[STIX]{x03C0}}$ , somewhat smaller than unity. For a general (non-shift) antisymmetric potential perturbation $\hat{\unicode[STIX]{x1D719}}$ one should take $\unicode[STIX]{x1D709}z_{A}=-\int \int \hat{\unicode[STIX]{x1D719}}(\unicode[STIX]{x1D70F}^{\prime })\,\text{d}\unicode[STIX]{x1D70F}^{\prime }\,\text{d}\unicode[STIX]{x1D70F}$ , in accordance with the remarks following (3.5). The fifth power dependence upon $\unicode[STIX]{x1D714}_{R}$ is one guarantee that $\ell =3$ and higher harmonic resonances can be ignored. So, henceforth we consider only $\ell =1$ . For small frequencies $\unicode[STIX]{x1D714}_{R}$ , the turning point position $z_{A}$ is approximately equal to the half-width of the hole, and varies only logarithmically at low frequency, because it occurs at $2\unicode[STIX]{x1D714}_{R}^{2}=-W=\unicode[STIX]{x1D719}_{0}(z)\simeq 16\unicode[STIX]{x1D713}\text{e}^{-|z|}$ in the wings, so for an exact shift mode

(3.13) $$\begin{eqnarray}z_{A}\simeq \ln \left(\frac{16\unicode[STIX]{x1D713}}{2\unicode[STIX]{x1D714}_{R}^{2}}\right)=\ln \left(\frac{8\unicode[STIX]{x1D713}^{-1/2}}{(\unicode[STIX]{x1D714}_{R}^{2}/\unicode[STIX]{x1D713}^{3/2})}\right)\simeq \ln \left(\frac{8\unicode[STIX]{x1D713}^{-1/2}}{(0.05)^{2}}\right)=8.-\ln (\unicode[STIX]{x1D713})/2,\end{eqnarray}$$

using the observed numerical scaling $\unicode[STIX]{x1D714}_{r}/\unicode[STIX]{x1D713}^{3/4}\sim 0.05$ . The value $z_{A}=9$ will be used as a generic estimate. The dominant dependency is $\tilde{F}_{R}\propto \unicode[STIX]{x1D709}\unicode[STIX]{x1D714}_{r}^{5}/\unicode[STIX]{x1D713}^{1/2}$ , and we will ignore the weak dependence on $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D714}_{r}$ of $z_{A}$ and of the correction factors $D$ and $A$ , in estimating the constant of proportionality. For $\unicode[STIX]{x1D713}=0.16$ , $D=0.76$ and $A$ can be as small as $0.9$ , in which case $D/A^{4}\simeq 1.15$ and the combined correction factors are moderate giving $\tilde{F}_{R}\simeq 1000i\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}_{r}^{5}/\sqrt{\unicode[STIX]{x1D713}}$ .

3.3 Non-resonant force

The non-resonant contribution to (3.5) is the sum of the non-resonant ( $\int _{0}^{t}$ ) integral plus the term $\text{i}\unicode[STIX]{x1D714}z(t)$ . If treated using the square-wave approximation $z(t)=\text{sign}(t)z_{A}$ it is

(3.14) $$\begin{eqnarray}\frac{q_{e}\unicode[STIX]{x1D6F7}_{NR}(t)}{\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D709}}=z(t)+(\text{i}\unicode[STIX]{x1D714})\int _{0}^{t}z(\unicode[STIX]{x1D70F})\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=\text{sign}(t)z_{A}\text{e}^{\text{i}\unicode[STIX]{x1D714}t},\end{eqnarray}$$

taking into account that the signs of $z$ and $t$ are the same. When multiplied by the antisymmetric quantity $-(\text{i}\unicode[STIX]{x1D714})^{2}q_{e}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)$ and integrated for positive and negative $z(t)$ and $v$ this gives rise to a force

(3.15) $$\begin{eqnarray}\tilde{F}_{NR}=-(\text{i}\unicode[STIX]{x1D714})^{2}\unicode[STIX]{x1D709}\int _{0}^{\sqrt{-2q_{e}\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\int _{0}^{t_{b}/2}z_{A}2\cos \unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime }q_{e}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\,\text{d}\unicode[STIX]{x1D70F}^{\prime }v_{\unicode[STIX]{x1D713}}\,\text{d}v_{\unicode[STIX]{x1D713}}.\end{eqnarray}$$

This expression is manifestly real when $\unicode[STIX]{x1D714}_{i}=0$ , and provides the main contribution to $\text{Re}(\tilde{F}_{t})$ . But it also makes an important contribution, first order in $\unicode[STIX]{x1D714}_{i}$ , to $\text{Im}(\tilde{F}_{t})$ : approximately $(2\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D714}_{r})\tilde{F}_{NR}(\unicode[STIX]{x1D714}_{r})$ .

It is less clear for the non-resonant contribution that the square-wave approximation to $z(t)$ is accurate, since contributions to this integral also arise from deeply trapped particles, which have almost sinusoidal bounce motion. Nevertheless, proceeding as before to replace $-q_{e}(\text{d}\unicode[STIX]{x1D719}_{0}/\text{d}z)=(\text{d}v/\text{d}\unicode[STIX]{x1D70F}^{\prime })$ , we get $-\int _{0}^{t_{b}/2}\cos \unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime }q_{e}(\text{d}\unicode[STIX]{x1D719}/\text{d}z)\,\text{d}\unicode[STIX]{x1D70F}^{\prime }=\int _{0}^{t_{b}/2}\cos \unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime }(\text{d}v/\text{d}\unicode[STIX]{x1D70F}^{\prime })\,\text{d}\unicode[STIX]{x1D70F}^{\prime }=-v_{\unicode[STIX]{x1D713}}+v(t_{b}/2)\cos (\unicode[STIX]{x1D714}t_{b}/2)+\int \unicode[STIX]{x1D714}\sin \unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime }v(\unicode[STIX]{x1D70F}^{\prime })\,\text{d}\unicode[STIX]{x1D70F}^{\prime }$ . The final integral term can be ignored by $\unicode[STIX]{x1D714}$ ordering and the second because $v(t_{b}/2)=0$ everywhere except on a negligible measure at $\unicode[STIX]{x1D714}t_{b}/2\simeq n\unicode[STIX]{x03C0}$ ; so we can substitute $-v_{\unicode[STIX]{x1D713}}$ for $-\int _{0}^{t_{b}/2}\cos \unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}^{\prime }q_{e}(\text{d}\unicode[STIX]{x1D719}/\text{d}z)\,\text{d}\unicode[STIX]{x1D70F}^{\prime }$ , giving

(3.16) $$\begin{eqnarray}\displaystyle \tilde{F}_{NR} & \simeq & \displaystyle -(\text{i}\unicode[STIX]{x1D714})^{2}\unicode[STIX]{x1D709}\int _{0}^{\sqrt{-2q_{e}\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}z_{A}2v_{\unicode[STIX]{x1D713}}v_{\unicode[STIX]{x1D713}}dv_{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & \simeq & \displaystyle -(\text{i}\unicode[STIX]{x1D714})^{2}\unicode[STIX]{x1D709}\left\langle \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}2z_{A}\right\rangle \frac{1}{3}(2\unicode[STIX]{x1D713})^{3/2}=\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D709}w_{t}\,f_{0s}\sqrt{2\unicode[STIX]{x03C0}}\frac{5}{8}\unicode[STIX]{x1D713},\end{eqnarray}$$

where $\langle (\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W)2z_{A}\rangle$ denotes an average over energy whose weighting $\propto v_{\unicode[STIX]{x1D713}}^{2}$ (emphasizing shallowly trapped orbits) justifies writing it as $w_{t}f_{0s}(15/16)\sqrt{\unicode[STIX]{x03C0}/\unicode[STIX]{x1D713}}$ , with $w_{t}\sim 2z_{A}$ a hole width somewhat greater than unity. (But the $z_{A}$ here is not just for resonant orbits; it is an average over all trapped particles.) This expression for unit density is $\frac{5}{8}w_{t}\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D709}\unicode[STIX]{x1D713}=5\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D709}\unicode[STIX]{x1D713}$ when $w_{t}=8$ which is a plausible estimate.

The non-resonant $\unicode[STIX]{x1D714}_{i}$ -term involves in addition a small contribution from passing particles. The total value $\text{Re}(\tilde{F}/\unicode[STIX]{x1D714}^{2})$ has been evaluated more precisely elsewhere (Hutchinson & Zhou Reference Hutchinson and Zhou2016; Hutchinson Reference Hutchinson2018b ), and for unit density can be written in terms of the function $h(\unicode[STIX]{x1D712})=-(2/\sqrt{\unicode[STIX]{x03C0}})\unicode[STIX]{x1D712}+(2\unicode[STIX]{x1D712}^{2}-1)\text{e}^{\unicode[STIX]{x1D712}^{2}}\text{erfc}(\unicode[STIX]{x1D712})+1=\unicode[STIX]{x1D712}^{2}-2/\sqrt{\unicode[STIX]{x03C0}}\frac{4}{3}\unicode[STIX]{x1D712}^{3}+O(\unicode[STIX]{x1D712}^{4}),$ where $\unicode[STIX]{x1D712}^{2}=\unicode[STIX]{x1D719}_{0}(z)$ , as

(3.17) $$\begin{eqnarray}\tilde{F}_{NR}=\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}^{2}\int h\left(\sqrt{\unicode[STIX]{x1D719}_{0}(z)}\right)\,\text{d}z=\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D713}w_{i}\simeq \unicode[STIX]{x1D709}\unicode[STIX]{x1D714}^{2}\int \unicode[STIX]{x1D719}_{0}(z)\,\text{d}z,\end{eqnarray}$$

which can be considered to define a width $w_{i}\simeq \int \unicode[STIX]{x1D719}_{0}\,\text{d}z/\unicode[STIX]{x1D713}$ . For shallow holes of the chosen shape, it is $w_{i}\simeq \int \,\text{sech}^{4}(z/4)\,\text{d}z=16/3$ , giving $\tilde{F}\simeq 5.3\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D713}$ , in excellent agreement with the estimate developed from (3.16), confirming that the passing particle contribution is relatively unimportant here.

3.4 Passing-particle imaginary force

To obtain the lowest-order imaginary part of the passing-particle force that is independent of $\unicode[STIX]{x1D714}_{i}$ we proceed in a similar manner using two integrations by parts, except with different constants of integration than (3.5).

(3.18) $$\begin{eqnarray}\frac{q_{e}\unicode[STIX]{x1D6F7}(t)}{\unicode[STIX]{x1D709}}=[v(t)-v_{\infty }]+\text{i}\unicode[STIX]{x1D714}[z(t)-z_{\infty }(t)]+(\text{i}\unicode[STIX]{x1D714})^{2}\int _{-\infty }^{t}[z(\unicode[STIX]{x1D70F})-z_{\infty }(\unicode[STIX]{x1D70F})]\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}.\end{eqnarray}$$

Here $v_{\infty }$ denotes the distant velocity (outside the hole) of the orbit under consideration; and $z_{\infty }(t)=v_{\infty }t+\text{const}.$ denotes the position as a function of time for motion at a constant velocity $v_{\infty }$ which extrapolates the distant past orbit ignoring $\unicode[STIX]{x1D719}(z)$ . Thus $[z(\unicode[STIX]{x1D70F})-z_{\infty }(\unicode[STIX]{x1D70F})]$ is zero in the distant past, rises as it passes through the hole and then remains constant past the hole. As before, it is the final integral term that gives the imaginary force (real part of $\unicode[STIX]{x1D6F7}$ ). We approximate the integral by treating $[z(\unicode[STIX]{x1D70F})-z_{\infty }(\unicode[STIX]{x1D70F})]$ as the unit step function $H(\unicode[STIX]{x1D70F})\equiv [1+\text{sign}(\unicode[STIX]{x1D70F})]/2$ times an amplitude $z_{p}=\int _{-\infty }^{\infty }[v(t)-v_{\infty }]\,\text{d}t$ , and note that $\int _{-\infty }^{t}H(\unicode[STIX]{x1D70F})\text{e}^{-\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D70F}-t)}\,\text{d}\unicode[STIX]{x1D70F}=H(t)[1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t}]/(-\text{i}\unicode[STIX]{x1D714})$ . Then, ignoring other terms that are subsequently annihilated by symmetry, we have for small $\unicode[STIX]{x1D714}t$ a term independent of $\unicode[STIX]{x1D714}_{i}$

(3.19) $$\begin{eqnarray}\text{Re}\{q_{e}\unicode[STIX]{x1D6F7}(t)\}\simeq \text{Re}\{(-\text{i}\unicode[STIX]{x1D714})\unicode[STIX]{x1D709}z_{p}[1-\text{e}^{\text{i}\unicode[STIX]{x1D714}t}]H(t)\}\simeq (\text{i}\unicode[STIX]{x1D714}_{r})^{2}\unicode[STIX]{x1D709}z_{p}tH(t).\end{eqnarray}$$

From the second term of (3.18) we also have a component of $q_{e}\unicode[STIX]{x1D6F7}(t)$ equal to $\text{i}\unicode[STIX]{x1D714}_{r}z_{p}H(t)$ , which could be pursued, but it contributes imaginary force only $\propto \unicode[STIX]{x1D714}_{i}$ and is subdominant to the similar trapped component, so we do not bother. The required imaginary contribution to passing force at $\unicode[STIX]{x1D714}_{i}=0$ is thus

(3.20) $$\begin{eqnarray}\displaystyle \text{i}\text{Im}(\tilde{F}_{p}) & \simeq & \displaystyle -\int \int q_{e}\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\unicode[STIX]{x1D709}(\text{i}\unicode[STIX]{x1D714})^{3}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}z_{p}\unicode[STIX]{x1D70F}^{\prime }H(\unicode[STIX]{x1D70F}^{\prime })\,\text{d}\unicode[STIX]{x1D70F}^{\prime }v_{\infty }\,\text{d}v_{\infty }\nonumber\\ \displaystyle & = & \displaystyle -q_{e}\unicode[STIX]{x1D709}(\text{i}\unicode[STIX]{x1D714})^{3}\int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}z_{p}\int _{0}^{\infty }-\unicode[STIX]{x1D70F}^{\prime }\frac{\text{d}v}{\text{d}\unicode[STIX]{x1D70F}^{\prime }}\,\text{d}\unicode[STIX]{x1D70F}^{\prime }v_{\infty }\,\text{d}v_{\infty }\nonumber\\ \displaystyle & = & \displaystyle -q_{e}\unicode[STIX]{x1D709}(\text{i}\unicode[STIX]{x1D714})^{3}\int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}z_{p}\int _{0}^{\infty }[v(\unicode[STIX]{x1D70F}^{\prime })-v_{\infty }]\,\text{d}\unicode[STIX]{x1D70F}^{\prime }\,v_{\infty }\,\text{d}v_{\infty }\nonumber\\ \displaystyle & = & \displaystyle -q_{e}\unicode[STIX]{x1D709}(\text{i}\unicode[STIX]{x1D714})^{3}\int \frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\frac{z_{p}^{2}}{2}v_{\infty }\,\text{d}v_{\infty }.\end{eqnarray}$$

The integral quantity $z_{p}(v_{\infty })=\int (1-|v_{\infty }/v|)\,\text{d}z$ is approximately the value of its integrand at the hole centre times the hole width: $z_{p}(v_{\infty })\simeq (1-v_{\infty }/\sqrt{v_{\infty }^{2}+2\unicode[STIX]{x1D713}})w_{p}$ . Its dependence for $v_{\infty }^{2}\gtrsim \unicode[STIX]{x1D713}$ is $z_{p}\simeq w_{p}\unicode[STIX]{x1D713}/v_{\infty }^{2}$ ; so $\int _{\unicode[STIX]{x1D713}}^{\infty }z_{p}^{2}\,\text{d}v_{\infty }^{2}\simeq w_{p}^{2}\unicode[STIX]{x1D713}$ , and for small $\unicode[STIX]{x1D713}$ we can multiply this by the value of $\text{d}f_{0}/\text{d}W$ at small energy. The remainder of the integral can be approximated as having constant $z_{p}\simeq w_{p}$ so $\int _{0}^{\unicode[STIX]{x1D713}}z_{p}^{2}\,\text{d}v_{\infty }^{2}\simeq w_{p}^{2}\unicode[STIX]{x1D713}$ , giving

(3.21) $$\begin{eqnarray}\int _{0}^{\infty }\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\frac{z_{p}^{2}}{4}\text{d}v_{\infty }^{2}\simeq \left.\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\right|_{v_{\infty }=0}w_{p}^{2}\frac{\unicode[STIX]{x1D713}}{2}\end{eqnarray}$$

(which can be considered the definition of $w_{p}$ ). Thus finally

(3.22) $$\begin{eqnarray}\text{Im}(\tilde{F}_{p})\simeq q_{e}\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}^{3}\left.\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}W}\right|_{0}\frac{w_{p}^{2}}{2}\unicode[STIX]{x1D713}.\end{eqnarray}$$

This is negative because $(\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W)|_{0}=-f_{0s}=-1/\sqrt{2\unicode[STIX]{x03C0}}$ for unit density. Also notice that despite depending on $\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W$ this is not a resonance effect but is contributed by all energies from zero to a few times $\unicode[STIX]{x1D713}$ .

4.1 Complex frequency

We now have analytic approximations of the three terms that at low frequency predominate the relation that the total imaginary part of the hole force must vanish,

(4.1) $$\begin{eqnarray}2\frac{\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x1D714}_{r}}\text{Re}(\tilde{F})+\text{Im}(\tilde{F}_{t})+\text{Im}(\tilde{F}_{p})=\unicode[STIX]{x1D709}[C_{i}\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{r}\unicode[STIX]{x1D713}+C_{t}\unicode[STIX]{x1D713}^{-1/2}\unicode[STIX]{x1D714}_{r}^{5}+C_{p}\unicode[STIX]{x1D714}_{r}^{3}\unicode[STIX]{x1D713}]=0.\end{eqnarray}$$

Here the primary $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D713}$ dependences have been written explicitly, leaving coefficients $C_{i}\simeq 10$ , $C_{t}\simeq 1000$ and $C_{p}\simeq -w_{p}^{2}/2\sqrt{2\unicode[STIX]{x03C0}}\simeq -7$ (for $w_{p}=6$ , a plausible value chosen with hindsight) for unit density $f_{0s}=1/\sqrt{2\unicode[STIX]{x03C0}}$ that are approximately constant. Writing $\hat{\unicode[STIX]{x1D714}}_{r}=\unicode[STIX]{x1D714}_{r}/\unicode[STIX]{x1D713}^{3/4}$ and $\hat{\unicode[STIX]{x1D714}}_{i}=\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D713}^{3/2}$ the relation becomes

(4.2) $$\begin{eqnarray}0=\text{Im}(\tilde{F})=\unicode[STIX]{x1D709}\unicode[STIX]{x1D713}^{13/4}\hat{\unicode[STIX]{x1D714}}_{r}[C_{i}\hat{\unicode[STIX]{x1D714}}_{i}+C_{t}\hat{\unicode[STIX]{x1D714}}_{r}^{4}+C_{p}\hat{\unicode[STIX]{x1D714}}_{r}^{2}],\end{eqnarray}$$

which demonstrates the universality (observed approximately in figure 2) when expressed in the scaled quantities $\hat{\unicode[STIX]{x1D714}}_{r}$ and $\hat{\unicode[STIX]{x1D714}}_{i}$ . It also reproduces the quantitative results of figure 2, and specifically the zero contour of $\text{Im}(\tilde{F})$ . First however, figure 4(a) demonstrates the agreement of the analytic formulas for the three imaginary components of the force with the corresponding parts of the force calculated by the numerical integration code. All the forces in that plot are normalized in the form $\hat{F}\equiv \tilde{F}/\unicode[STIX]{x1D709}\unicode[STIX]{x1D713}^{13/4}$ , in accordance with (4.2), giving an essentially universal plot, independent of $\unicode[STIX]{x1D713}$ (when small) expressed in terms of $\hat{\unicode[STIX]{x1D714}}$ . The agreement of the analytic curves with the appropriate range of the full numerical integration is very good, confirming that the dominant terms (of (4.1)) have been correctly identified and quantitatively estimated.Footnote ¹

Figure 4. (a) Numerical integration (solid lines) agrees asymptotically with the analytic approximations (dashed lines) of the three terms in the imaginary force. The curve is universal when expressed in scaled quantities $\hat{\unicode[STIX]{x1D714}}_{r}$ , $\hat{\unicode[STIX]{x1D714}}_{i}$ and $\hat{F}$ . (b) The contours of force resulting from the analytic expressions, showing good shape agreement with figure 2.

Figure 4(b) shows the shape of the force contours for comparison with figure 2. The agreement is within the level of variation expected. The position of the bottom right-hand end of the zero contour, at $\unicode[STIX]{x1D714}_{i}=0$ , is decided by setting equal and opposite the second and third terms, whose solution can be rendered as

(4.3) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D714}}_{r}=\sqrt{-C_{p}/C_{t}}\simeq \sqrt{7/1000}=0.084.\end{eqnarray}$$

Near the bottom left-hand end we can neglect the trapped term and find

(4.4) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D714}}_{i}\simeq (-C_{p}/C_{i})\hat{\unicode[STIX]{x1D714}}_{r}^{2}\simeq 0.7\hat{\unicode[STIX]{x1D714}}_{r}^{2}.\end{eqnarray}$$

This parabolic $\hat{\unicode[STIX]{x1D714}}_{r}$ -dependence appears to be present in figure 2 (although numerical rounding and other uncertainties make the zero contour shape imprecise) and is obvious in figure 4(b). The top of the contour is determined by maximizing $-C_{t}\hat{\unicode[STIX]{x1D714}}_{r}^{4}-C_{p}\hat{\unicode[STIX]{x1D714}}_{r}^{2}$ which gives

(4.5) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D714}}_{r}=\sqrt{\frac{-C_{p}}{2C_{t}}}\simeq 0.06\end{eqnarray}$$

and then

(4.6) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D714}}_{i}=\frac{1}{2}\left(\frac{-C_{p}}{C_{i}}\right)\unicode[STIX]{x1D714}_{r}^{2}=\frac{1}{4}\left(\frac{-C_{p}}{C_{i}}\right)\left(\frac{-C_{p}}{C_{t}}\right)\simeq 1.2\times 10^{-3}.\end{eqnarray}$$

These expressions agree with the corresponding positions on the zero contour in figure 4(b) (as they must).

It is clear from these observations that the effect of the resonant force is to limit the maximum $\unicode[STIX]{x1D714}_{r}$ at which instability can occur. In other words, it is a stabilizing term (contradicting Newman et al. (Reference Newman, Goldman, Spector and Perez2001), Vetoulis & Oppenheim (Reference Vetoulis and Oppenheim2001) and Berthomier et al. (Reference Berthomier, Muschietti, Bonnell, Roth and Carlson2002)), which as $\unicode[STIX]{x1D714}_{r}$ increases eventually overcomes the others because it varies as a higher power. If its coefficient $C_{t}$ were reduced, for example by decreasing the distribution gradient $\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W$ for shallowly trapped particles, or for a perturbation whose shift was smaller in the wings of the hole, the effect would be to permit instability to higher $\unicode[STIX]{x1D714}_{r}$ and with higher $\unicode[STIX]{x1D714}_{i}$ , following the $\hat{\unicode[STIX]{x1D714}}_{i}\simeq 0.35\hat{\unicode[STIX]{x1D714}}_{r}^{2}$ trajectory further up. The stabilizing nature of the resonance might appear counter-intuitive since $\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}W$ is positive causing resonant particles to give energy to the eigenmode. What makes this effect stabilizing, however, is that the eigenmode is a ‘negative energy’ mode (Lashmore-Davies Reference Lashmore-Davies2005), reflecting the negative inertia of the hole. Therefore adding energy to it tends to reduce its amplitude.

4.2 Wavenumber

As explained earlier, the wavenumber must be chosen to satisfy the real part of (2.7) $\text{Re}(\tilde{F})=F_{E}$ . The real frequency $\unicode[STIX]{x1D714}_{r}$ gives the predominant real part of $\tilde{F}$ at low frequency, namely $\text{Re}(\tilde{F})=(C_{i}/2)\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}_{r}^{2}\unicode[STIX]{x1D713}$ , see (3.17) and (4.1). The Maxwell stress force can be evaluated in dimensionless units as

(4.7) $$\begin{eqnarray}F_{E}=\unicode[STIX]{x1D709}k^{2}\int \left(\frac{\text{d}\unicode[STIX]{x1D719}_{0}}{\text{d}z}\right)^{2}\,\text{d}z=\frac{128}{315}\unicode[STIX]{x1D709}k^{2}\unicode[STIX]{x1D713}^{2},\end{eqnarray}$$

from which we conclude

(4.8) $$\begin{eqnarray}k^{2}=\frac{315}{128}\frac{C_{i}}{2}\frac{\unicode[STIX]{x1D714}_{r}^{2}}{\unicode[STIX]{x1D713}},\quad \text{giving}~\hat{k}=\frac{k}{\unicode[STIX]{x1D713}^{1/4}}=\sqrt{\frac{315C_{i}}{256}}\hat{\unicode[STIX]{x1D714}}_{r}=3.5\hat{\unicode[STIX]{x1D714}}_{r}.\end{eqnarray}$$

At the maximum unstable $\unicode[STIX]{x1D714}_{i}$ , $\hat{\unicode[STIX]{x1D714}}_{r}\simeq 0.06$ , we then require $\hat{k}\simeq 0.2$ , which is in reasonable agreement with figure 2. We have recovered the observed scaling of the universal solution: $k=\hat{k}\unicode[STIX]{x1D713}^{1/4}$ , and approximately the absolute value of $\hat{k}$ required for maximum instability growth.