2.1 Adiabatic response eigenmodes

In units scaled to Debye length, electron background density and electron thermal energy, the one-dimensional Poisson equation, assuming immobile uniform ion density, is ${\rm d}^2\phi /{\rm d}z^2=n-1$. The equilibrium analysed is chosen to be of the form $\phi _0=\psi {\,{\rm sech}}^4(z/4)$. Writing for brevity $S={\,{\rm sech}}(z/4)$ and $T=\tanh (z/4)$ and denoting the $z$ derivative by prime, the linearized equilibrium (‘adiabatic’) Poisson operator for this hole potential can then be written

(2.1)\begin{equation} V_a = \left[\frac{{\rm d}^2}{{\rm d}z^2} - \frac{{\rm d}n}{{\rm d}\phi_0}\right] = \left[\frac{{\rm d}^2}{{\rm d}z^2} - \frac{\phi_0^{'''}}{\phi_0'}\right] = \left[\frac{{\rm d}^2}{{\rm d}z^2} +\frac{30}{16}S^2- 1\right]. \end{equation}

Notation $V$ reminds us that this is partly the Vlasov operator, transforming a potential into a density, and the subscript $a$ denotes adiabatic (meaning steady). The eigenmodes of this operator, satisfying $V_a|u\rangle =\lambda |u\rangle$, can be found (see Appendix A) by applying the raising operator $4 ({{\rm d}}/{{\rm d}z})-l T$ for $l=1,2,3,4,5$ to the function ${\rm e}^{uz/4}$ yielding

(2.2)$$\begin{gather} |u\rangle(z)= \exp(uz/4)\{[{-}15u^4 + (420S^2 - 225)u^2 - 945S^4 + 840S^2 - 120]T \nonumber\\ +\, u[u^4 + ({-}105S^2 + 85)u^2 + (945S^4 - 1155S^2 +274)]\}. \end{gather}$$

For real $u$, as $u z\to +\infty$ the mode is bounded (and tends to zero) only if $u$ is one of the discrete roots of the polynomial in the braces obtained by letting $S\to 0$ and $T\to {\rm sign}(u)$. These are $u=j=1,2,3,4,5$. The odd-numbered discrete modes are symmetric in $z$, and the even-numbered are antisymmetric. In contrast, imaginary $u$ values $u={\rm i}p$ give rise to the continuum modes which are formally finite at infinity. Their overlap integral over a finite domain exists; and, as the domain tends to infinity, it tends to a delta function. The continuum eigenmodes have definite parity only if $u$ reverses sign with $z$; so we write $u={\rm i} \sigma _z p$, where $\sigma _z={\rm sign}(z)$ and positive $p$ represents antisymmetric outwardly propagating waves.

The corresponding eigenvalues are

(2.3)\begin{equation} \lambda_u= u^2/16-1={-}p^2/16-1. \end{equation}

Normalized so that $\langle j |l\rangle =\delta _{jl}$, they are given in table 1.

Table 1.

Discrete eigenmodes.

In particular $\lambda _4=0$ for the shiftmode $|4\rangle \propto {\rm d}\phi _0/{\rm d}z$, and it is the predominant perturbation in essentially all linear electron hole instabilities. Moreover, it couples only to antisymmetric modes; so $|2\rangle$ is the only other discrete mode that needs to be considered.

It can be shown that the the continuum modes are ‘normalized’, in the sense that then $\langle p |q\rangle =\delta _{pq}=\delta (p-q)$, by dividing (2.2) by the factor

(2.4)\begin{equation} [8{\rm \pi}(p^2+1^2)(p^2+2^2)(p^2+3^2)(p^2+4^2)(p^2+5^2)]^{1/2}. \end{equation}

Far from the hole ($|z|\gg 1$) the normalized oscillatory continuum modes are sinusoidal with amplitude $1/\sqrt {8{\rm \pi} }=0.19947$ and parallel wavenumber $|k_\parallel |=p/4$. The eigenmodes are plotted in figure 1.

Figure 1.

Eigenmodes of the adiabatic response operator.

2.2 Including non-adiabatic response

We now must consider the full linearized Poisson equation including the non-adiabatic response arising from the solution to the time-dependent Vlasov equation, $\tilde {V}$, as well as the adiabatic response $V_a$. The form of the non-adiabatic Vlasov operator $\tilde {V}$ is discussed later. For a perturbation $\phi _1(y,z,t)=\hat \phi (z)\exp ({{\rm i}(k_\perp y-\omega t)})$ with transverse wavenumber $k_\perp$ and frequency $\omega$, Poisson's equation becomes

(2.5)\begin{equation} ({-}k_\perp^2+V_a+\tilde{V})|\hat\phi\rangle=0, \end{equation}

in which it is convenient to regard $k_\perp ^2\equiv \lambda _\perp$ as the full eigenvalue. We suppose the solution for potential can be expanded as a sum of scalar amplitudes $a_u$ times the eigenmodes $|u\rangle$ of $V_a$:

(2.6)\begin{equation} |\hat\phi\rangle=\sum_u a_u |u\rangle, \end{equation}

where the summation notation also includes an integral over the continuum modes.

Then we can invoke the orthogonal properties of the adiabatic eigenmodes and form the inner product:

(2.7)\begin{equation} \langle s|( V_a+\tilde{V})|\hat\phi\rangle = \lambda_s\langle s|s\rangle a_s+\sum_u\langle s|\tilde{V}|u\rangle a_u, \end{equation}

which must be equal to $\lambda _\perp \langle s |s\rangle a _s$ to satisfy (2.5). In particular, choosing the predominant mode $\langle 4|$ for $\langle s |$, for which $\lambda _4=0$, we get an eigenvalue equation $0=(\lambda _4-\lambda _\perp )\langle 4|4\rangle a _4+\sum _u\langle 4|\tilde {V}|u\rangle a _u$.

If all the amplitudes $a_u$ for $u\neq 4$ are negligible, then $\tilde {V}$ contributes a correction ($\lambda _{4}^{(1)}$) to the eigenvalue:

(2.8)\begin{equation} \lambda_{4}^{(1)}= \frac{\langle4|\tilde{V}|4\rangle}{\langle4|4\rangle}=\lambda_\perp{-}\lambda_4. \end{equation}

The expression $\langle 4|\tilde {V}|4\rangle /\langle 4|4\rangle$ is the ‘Rayleigh quotient’ approximation for the eigenvalue of $V_a+\tilde {V}$ (because $\lambda _4=0$). It is physically the (normalized) jetting force on particles because of the unperturbed electric field $-{\rm d}\phi _0/{\rm d}z\propto \langle 4|$, acting on the non-adiabatic density perturbation $\tilde n = \tilde {V}|4\rangle$, integrated over the entire hole. It balances the (normalized) Maxwell shear stress from the transverse kinking of the hole ($k_\perp ^2$) to make the total force zero.

Actually $\tilde {V}$ is a complicated nonlinear function of the complex frequency $\omega$ of the mode; and for specified $k_\perp$ the dispersion relation between $\omega$ and $k_\perp$ must be solved by some kind of iterative procedure searching for an $\omega$ that satisfies (2.8). The imaginary part of $\omega$ thus found determines the stability of the hole. This approximation has yielded stability results (Hutchinson Reference Hutchinson2018b, Reference Hutchinson2019a,Reference Hutchinsonb) that are in reasonable (but not perfect) agreement with simulation. The question at hand is whether analysis can determine approximately the magnitude of the other coefficients $a_u$ for $u\neq 4$, and therefore give a more accurate perturbation structure $|\hat \phi \rangle$ and $\omega$.

If instead of the approximation (2.8) we are able to evaluate all the matrix coefficients $\langle s |\tilde {V}|u\rangle$, then in principle we can regard (2.7) instead of (2.8) as a matrix eigensystem we must solve to find the $\omega$ that permits a non-zero solution for the vector $a_u$. The off-diagonal matrix entries $s\neq u$ are the coupling of the potential modes by the non-adiabatic Vlasov operator $\tilde {V}$. The condition for the existence of a solution is that the determinant of the matrix $[(-\lambda _\perp +\lambda _s)\langle s |s\rangle \delta _{su}+\langle s |\tilde {V}|u\rangle ]$ should be zero.

Such a programme faces formidable practical challenges, however, because each evaluation of $\tilde {V}|u\rangle$ requires a computation involving multiple-dimension integrations over space and velocity distribution – repeated for each mode $|u\rangle$ and each adjustment of $\omega$. Moreover, in principle, the continuum contains infinitely many modes, and the matrix contains the square of the number of modes. Obviously we require this number to be reduced. We can immediately restrict our attention just to antisymmetric modes, since by assumed symmetry no coupling to symmetric modes occurs. We later show how the entire continuum contribution can be reduced to a single amplitude. Continuum modes $|p\rangle$ extend to $|z|=\infty$, an integration range that for computation needs to be reduced, and have formally divergent inner products. We show how these difficulties are negotiated.

Formally, any mode for which $\tilde {V}|u\rangle a _u$ is not negligible must be retained in the sum of $\langle s |\tilde {V}|u\rangle$ in (2.7), giving for the dominant mode

(2.9)\begin{equation} 0=(\lambda_4-\lambda_\perp)\langle4|4\rangle a_4+\langle4|\tilde{V}|4\rangle a_4+\sum_{u\neq 4}\langle4|\tilde{V}|u\rangle a_u, \end{equation}

but also for $s\neq 4$

(2.10)\begin{equation} 0=(\lambda_s-\lambda_\perp)\langle s|s\rangle a_s+\langle s|\tilde{V}|4\rangle a_4+\sum_{u\neq 4}\langle s|\tilde{V}|u\rangle a_s. \end{equation}

The well-known approach of time-independent perturbation theory in elementary quantum mechanics (see e.g. Dirac Reference Dirac1958, § 43) regards $\tilde {V}$ as systematically small, and takes $a_s$ for $s\neq 4$ also to be first-order small relative to $a_4$. The first-order approximation of (2.10) then drops the final sum term, giving

(2.11)\begin{equation} a_s \simeq \frac{\langle s|\tilde{V}|4\rangle a_4}{(\lambda_\perp{-}\lambda_s)\langle s|s\rangle}. \end{equation}

Substituting back (with $s\to u$) into (2.9), the eigenvalue to second order is

(2.12)\begin{equation} \lambda_\perp\simeq \lambda_4+\lambda_{4}^{(1)}+\lambda_{4}^{(2)}=\lambda_4 +\frac{\langle4|\tilde{V}|4\rangle}{\langle4|4\rangle} +\sum_{u\neq 4} \frac{\langle4|\tilde{V}|u\rangle\langle u|\tilde{V}|4\rangle}{(\lambda_\perp{-}\lambda_u)\langle4|4\rangle\langle u|u\rangle}. \end{equation}

Normally the substitution $\lambda _\perp \simeq \lambda _4$ is made in the final sum; but we do not need to do that since we consider $\lambda _\perp$ to be given and $\omega$ to be changed to achieve equality in this equation.

These perturbation equations (2.11) and (2.12) are appropriate for coupled discrete modes $|j\rangle$ when the amplitudes $a_u\langle u |u\rangle$ for $u\neq 4$ are small compared with $a_4\langle 4|4\rangle$. Those equations give the first-order mode amplitude and second-order eigenvalue correction from them. However, in our case, for the continuum modes, $\tilde {V}$ is not systematically small everywhere, and an altered approach appropriate to the actual form of $\tilde {V}$ must be adopted. We shall see that in order to obtain converged integrals for the relevant continuum inner products it is necessary to use the difference between $\tilde V |q\rangle$ and a pure wave operator expression $\tilde V_w|q\rangle$.

3.1 The non-adiabatic linearized Vlasov operator

The operator $\tilde {V}$ transforms a potential perturbation $|u\rangle$ into a non-adiabatic density perturbation $\tilde n$, both of which are complex functions of $z$. It does so by solving the linearized time-dependent Vlasov equation for the non-adiabatic distribution function perturbation $\tilde f(z,v)$ and integrating it: $\tilde n =\int \tilde f \,{\rm d}v$. The solution can be found in terms of an integral over past time along the linearized Vlasov equation's ‘characteristic’, that is, the unperturbed orbit $\boldsymbol x(t)$, giving an integral expression

(3.1)\begin{equation} \varPhi(\boldsymbol x,\boldsymbol v,t)\equiv \int_{-\infty}^t \phi_1(\boldsymbol x(\tau),t-\tau ) \,{\rm d}\tau, \end{equation}

where $\phi _1(\boldsymbol x,t)$ is the potential perturbation and $\boldsymbol x(\tau )$ is the past position of the unperturbed orbit that has velocity $\boldsymbol v$ at $(\boldsymbol x,t)$. The units of time adopted are $1/\omega _{pe}$: the inverse of the electron plasma frequency. When a uniform magnetic field in the $z$ direction is present with electron cyclotron frequency $\varOmega$ and the background perpendicular velocity distribution is Maxwellian, then the transverse perpendicular velocity dependence can be expressed as a sum over cyclotron harmonics (Hutchinson Reference Hutchinson2018b). The non-adiabatic parallel distribution function perturbation is then

(3.2)\begin{equation} \tilde f(z,y,v,t) = \exp({{\rm i}(k_\perp y-\omega t)}) \sum_{m={-}\infty}^\infty \tilde f_m(z,v)= \exp({{\rm i}(k_\perp y-\omega t)}) \sum_{m={-}\infty}^\infty b_m \varPhi_m, \end{equation}

with

(3.3a,b)\begin{equation} b_m={\rm i}\left[\omega_m \frac{\partial f_{\parallel0}}{\partial W_\parallel} +m\varOmega \frac{f_{\parallel0}}{T_\perp}\right] q_e{\rm e}^{-\zeta_t^2}{\rm I}_m(\zeta_t^2),\quad \varPhi_m=\int_{-\infty}^t \hat\phi(z(\tau)){\rm e}^{-{\rm i}\omega_m(\tau-t)}\,{\rm d}\tau. \end{equation}

Here $\zeta _t^2=k^2T_\perp /\varOmega ^2m_e$, where $T_\perp$ is the perpendicular temperature, ${\rm I}_m$ is the modified Bessel function and $\omega _m=m\varOmega +\omega$. The unperturbed equilibrium parallel distribution function depends only on parallel energy $W_\parallel$ (not $z$ directly) and is written $f_{\parallel 0}$, but for the perturbed parallel distributions we omit the $\parallel$ subscript for brevity, and from now on omit the $y,t$ dependencies as implicitly $\exp ({{\rm i}(k_\perp y-\omega t)})$ (so in $\varPhi _m$, the upper limit is $t=0$). The prior integral $\varPhi _m$ is a function of parallel position $z$ and velocity $v$. The weight $b_m$ is independent of $z$ but depends on $W_\parallel$ through the $f_{\parallel 0}$ distribution. We can regard each harmonic $m$ as giving a perturbed density contribution $\tilde n_m(z)=\int \tilde f_m\, {\rm d}v$. When $\hat \phi$ is a sum over modes $|u\rangle$,

(3.4)\begin{equation} \tilde n_m=\tilde{V}_{m}|\hat\phi\rangle=\tilde{V}_{m}\sum_u a_u|u\rangle=\int \sum_u a_u \tilde f_{um} \,{\rm d}v=\sum_u a_u \tilde n_{um}, \end{equation}

where $\tilde f_{um}$ denotes $\tilde f_m$ with $|\hat \phi \rangle =|u\rangle$ substituted in (3.2), and similarly $\tilde n_{um}=\int \tilde f_{um}\,{\rm d}v$. Then $\tilde V |\hat \phi \rangle =\sum _m\tilde V_m|\hat \phi \rangle =\sum _{u,m}a_u\tilde n_{um}$.

Solving analytically for $\varPhi _m$ and $n_m$ in the non-uniform equilibrium of an electron hole seems too difficult. The approach taken here is to perform the required integrations numerically. Evaluation of $\tilde n_{um}$ is performed using the same code for each $m$, but using the different $\omega _m$ in (3.3a,b). The inner products we need are $\langle s |\tilde V|u\rangle =\int _{-\infty }^\infty s^*(z)\tilde n_{u}(z)\,{\rm d}z$.

3.2 The external continuum wave contribution

In the region far outside the hole, $|z/4|\gg 1$ where $T=\tanh (z/4)\to \pm 1$ and $S={\,{\rm sech}}(z/4)\to 0$ in (2.2), the normalized continuum modes (taking them to be antisymmetrically outward-propagating) are

(3.5)\begin{align} |p\rangle=A{\rm e}^{{\rm i}\sigma_zpz/4}=\sigma_z\frac{({-}15p^4+225p^2-120) +{\rm i}p (p^{4}-85p^{2}+274)}{[8{\rm \pi}(p^2+1^2)(p^2+2^2)(p^2+3^2)(p^2+4^2)(p^2+5^2)]^{1/2}}{\rm e}^{{\rm i}\sigma_zpz/4}, \end{align}

where $\sigma _z={\rm sign}(z)$. Thus they are purely sinusoidal waves. For large enough $|z|$ the influence of the local hole potential becomes negligible, and such waves should there be identified with the normal modes of the uniform background plasma. The applicable normal mode for the present electrostatic approximation, in the (usually well-justified) cold electrostatic limit, ignoring ion response, has dispersion relation (including the upper hybrid waves)

(3.6)\begin{equation} k_\parallel^2/k_\perp^2=\frac{\omega^2(\varOmega^2+\omega_{pe}^2-\omega^2)} {(\varOmega^2-\omega^2)(\omega_{pe}^2-\omega^2)}. \end{equation}

(In our normalized units $\omega _{pe}=1$.) Assuming that $\omega$ and $k_\perp$ are prescribed, there is only one $|k_\parallel |$ that satisfies this dispersion relation; and corresponding to it, the continuum eigenmode of $V_a$ has $k_\parallel =\sigma _zp/4$. The mode number $p$ is taken positive, but $k_\parallel$ is signed. In this subsection we analyse just this single mode. The approach is mathematically justified in the following subsection.

For a single continuum mode, one can immediately calculate the value of $\varPhi _m(z)$ for external positions and inward orbit velocity (e.g. negative $v$ on the positive-$z$ side) using $z'=z(\tau )$ and $\tau =-(z-z')/v$, as

(3.7)\begin{align} \varPhi_{m}^{{\rm wave}}(z) = \int_{\sigma_z\infty}^zA\exp({{\rm i}[\sigma_zpz'/4+\omega_m (z-z')/v]})\,{\rm d}z'/v =\frac{A{\rm e}^{{\rm i}\sigma_zpz/4} }{{\rm i}(\sigma_zpv/4-\omega_m)}=\frac{A{\rm e}^{{\rm i}k_\parallel z} }{{\rm i}(k_\parallel v-\omega_m)}, \end{align}

where dropping the infinity limit is justified by positive imaginary part of $\omega$. This $\varPhi _{m}^{{\rm wave}}$ for the external region applies for all $z$ down to where the eigenmode begins to deviate significantly, because of non-zero hole potential, from the external wave; obviously that $z$ is of the order of a few times 4.

The non-adiabatic distribution perturbation for outward (positive $v$ at positive $z$) orbit velocity, by contrast, is strongly affected by the rapid variation of potential across the hole at $z\sim 0$ whose effect is carried into the external region by the particle orbits. However, for large enough $|z|$ the effects of the hole and earlier parts of the orbit become negligible because of dephasing between oscillatory contributions from different velocities. The dephasing effect attenuates the influence on density in a distance of several times $z_l\sim 1/\omega _m$, which is finite but (for $m=0$ at least) typically exceeds the hole extent itself (${\sim }4$) because $\omega$ is small. Thus the very distant ($|z|\gg z_l$) $\tilde n$ perturbations for inward- and outward-going velocities are given approximately by integrating ${\rm d}v$ the same $\varPhi _m^{{\rm wave}}$ expression for $z\gg z_l$, but with $v$ of opposite sign.

The total wave density perturbation can then be considered to be a constant $\tilde V_m^{{\rm wave}}$ times $|p\rangle$, where

(3.8)\begin{equation} \tilde V_m^{{\rm wave}}= \int_{-\infty}^\infty \frac{ b_m(v) }{{\rm i}(k_\parallel v-\omega_m)} \,{\rm d}v. \end{equation}

If $f_{\parallel 0}$ is an unshifted Maxwellian, and only $m=0$ is included, the total wave operator $\tilde V_w\equiv \sum _m\tilde V_m^{{\rm wave}}$ is proportional to the plasma dispersion function $Z$, and in the small $k_\parallel /\omega$ limit becomes $\tilde V_w=1+(k_\parallel /\omega )^2$. However, that approximation effectively implements the limit $\varOmega \gg \omega _{pe}$, and gives a dispersion relation $k_\parallel ^2/k_\perp ^2 = \omega ^2/(\omega _{pe}^2-\omega ^2)$, rather than the full cold plasma expression for finite $\varOmega /\omega$: (3.6). The full expression arises when the $|m|=1$ terms in the harmonic sum for the kinetic electrostatic dispersion relation (see e.g. Swanson Reference Swanson1989, § 4.4) are also included, to lowest order in $\zeta _t^2$. The resulting analytic form is

(3.9)\begin{equation} \tilde V_w= 1 + \left(\frac{k_\parallel}{\omega}\right)^2 + k_\perp^2 \frac{T_\perp{/}T_\parallel}{\omega^2-\varOmega^2}. \end{equation}

The crucial point is that $(\tilde {V}-\tilde {V}_{w})|q\rangle$ is localized, tending to zero in the region beyond $z_l$, which means that overlap integrals of the form $\langle p |\tilde {V}-\tilde {V}_{w}|q\rangle$ exist finite over an infinite domain. An important secondary feature is that, for external positions $|z|\ge z_d$ (where the effective hole edge $z_d$ is great enough that $\phi _0(z_d)=0$), the inward velocity part of $\tilde V_m^{{\rm wave}}$, consisting of

(3.10)\begin{equation} \tilde V_m^{{\rm in}}= \int_{0}^\infty \frac{ b_m(-\sigma_x|v|) }{{\rm i}(-|k_\parallel v|-\omega_m)}\sigma_x \,{\rm d}|v|, \end{equation}

is exactly equal to the actual non-adiabatic density perturbation that accounts for the hole's presence and the full eigenmode structure. So the cancellation of the part of $\tilde V-\tilde V_w$ for inward orbits only is exact $(\tilde V^{{\rm in}}-\tilde V_m^{{\rm in}})|q\rangle =0$. Therefore the non-zero contribution to external region $z$ integrals of $(\tilde V-\tilde V_w)|q\rangle$ arises only from outward orbits $(\tilde V^{{\rm out}}-\tilde V_m^{{\rm out}})|q\rangle \neq 0$. The necessary integral of outward orbits, however, must be carried out to an upper limit for which $|z|\gg z_l\gg z_d$. We consider explicitly only parallel distributions $f_{\parallel 0}$ that are symmetric in $v$. In that case $b_m$ is symmetric, and the overlap integrals can be calculated for a single sign of $v$ and then doubled to give the total.

3.3 Reducing the continuum to give the dispersion matrix

Now we discuss mathematically how the previous considerations allow us to reduce the continuum contribution to effectively a single mode. Using $s,u\to q,p$ etc. to observe our notation for continuum modes, and supposing them to be normalized ($\langle q |p\rangle =\delta _{qp}$), we write the continuum part of (2.10) explicitly as an integral $a(p)\,{\rm d}p$ over an amplitude distribution function $a(p)$ rather than a sum. Proceeding with no assumption about the size of cross-coupling of secondary modes, the eigenmode equation inner product with mode $\langle s |$ gives the re-expressed (2.10) as

(3.11)\begin{equation} 0=(\lambda_s-\lambda_\perp)\langle s|s\rangle a_s+\sum_j\langle s|\tilde{V}|j\rangle a_j +\int\langle s|\tilde{V}|p\rangle a(p)\,{\rm d}p. \end{equation}

We apply this equation for the three modes $s=4,2,q$. First for the continuum mode ($\langle s |=\langle q |$) with the adiabatic terms moved to the left-hand side:

(3.12)\begin{equation} (-\lambda(q)+\lambda_\perp)a(q)=\sum_j\langle q|\tilde{V}|j\rangle a_j +\int\langle q|\tilde V|p\rangle a(p)\,{\rm d}p. \end{equation}

We write the continuum integral using $\langle q |\tilde V_w|p\rangle = \tilde V_w \delta _{qp}$ as

(3.13)\begin{align} \int \langle q|\tilde{V}|p\rangle a(p)\,{\rm d}p=\int\langle q|\tilde{V}_{w}+(\tilde{V}-\tilde{V}_{w})|p\rangle a(p)\,{\rm d}p =\tilde{V}_{w}a(q) +\int\langle q|(\tilde{V}-\tilde{V}_{w})|p\rangle a(p)\,{\rm d}p. \end{align}

The final term can be approximated by observing that, by construction, $\tilde V -\tilde V_w$ is significant only in the inner $z$ region. As can be seen in figure 1, the form of the continuum modes in the inner region is almost independent of $p$. (Actually for small $p$, which is our interest here, the differences are even smaller than that figure shows.) Therefore, to that degree of approximation, we can replace $|p\rangle$ with $|q\rangle$ in the final term, giving $\int \langle q |(\tilde {V}-\tilde {V}_{w})|p\rangle a (p)\,{\rm d}p\simeq \langle q |\tilde {V}-\tilde {V}_{w}|q\rangle \int a(p)\,{\rm d}p$, and obtain the equation

(3.14)\begin{equation} (-\lambda(q)+\lambda_\perp{-}\tilde V_w)a(q)=\sum_j\langle q|\tilde{V}|j\rangle a_j +\langle q|\tilde V-\tilde V_w|q\rangle\int a(p)\,{\rm d}p. \end{equation}

Since (by the approximate invariance of the continuum modes in the inner region) the right-hand side is a weak function of $q$, this equation implies that $a(q)$ is a resonant function of $q$ centred on the $q$ value for which its coefficient on the left-hand side is approximately zero. We can write the coefficient of $a(q)$ using (2.3) $-\lambda (q)=k_\parallel ^2+1$ and (3.9) as

(3.15)\begin{align} k_\parallel^2+k_\perp^2+1-\tilde V_w\simeq k_\perp^2\left(1-\frac{T_\perp{/}T_\parallel}{\omega^2-\varOmega^2}\right) +k_\parallel^2\left(1-\frac{1}{\omega^2}\right) = \frac{1}{4^2}\left(\frac{1}{\omega^2}-1\right)(q_0^2-q^2), \end{align}

where

(3.16)\begin{equation} \left. q_0=4k_\perp\sqrt{1+\frac{T_\perp{/}T_\parallel}{\varOmega^2-\omega^2}} \right/\sqrt{\frac{1}{\omega^2}-1}. \end{equation}

The coefficient $(-\lambda (q)+\lambda _\perp -V_w)$ is zero when $q^2=q_0^2$, which is the dispersion relation of the wave in a uniform plasma. It is equal to the cold plasma expression (3.6) when $T_\perp /T_\parallel =1$.

We integrate (3.14) ${\rm d}q/(\lambda _\perp -\lambda (q)-V_w)$, recognizing again the approximate independence of $q$ of the right-hand side to find

(3.17)\begin{equation} \int a(q) \,{\rm d}q =\int \frac{{\rm d}q}{\lambda_\perp{-}\lambda(q)-V_w} \left(\sum_j\langle q|\tilde{V}|j\rangle a_j +\langle q|\tilde V-\tilde V_w|q\rangle\int a(p)\,{\rm d}p \right). \end{equation}

In view of the resonant form of the expression for $a(q)$, we can regard the integral over this resonance $\int a(q) \,{\rm d}q$ as quantifying the total continuum perturbation. The integral encounters the two poles of the integrand at $q=\pm q_0$. Near them the real part of $(\lambda _\perp -\lambda (q)-V_w)$ becomes small. Its imaginary part comes from the small, necessarily positive, imaginary part of $\omega =\omega _r+{\rm i}\omega _i$. It can then quickly be shown that $q_0$ has a small positive imaginary part: $q_0=q_{0r}+{\rm i}\epsilon$.

The infinite integral along the real $q$ axis can be closed by returning in the upper part of the complex plane ($\mathrm {Im}(q)$ positive) where the eigenmode tends to zero. The resonant integral is then just $2{\rm \pi} {\rm i}$ times the residue at the positive $q_0$ pole, which is the mathematical justification for our physical assumption of outwardly propagating waves. Writing it $\int ({{\rm d}q}/{\lambda _\perp -\lambda (q)-V_w})=1/K$, we find

(3.18)\begin{equation} K\simeq\frac{(1/\omega^2-1)q_0}{16{\rm \pi} {\rm i}} = \frac{k_\perp}{4{\rm \pi} {\rm i}}\sqrt{\frac{1}{\omega^2}-1} \sqrt{1+\frac{T_\perp{/}T_\parallel}{\Omega^2-\omega^2}}. \end{equation}

Substituting elsewhere the value $q=q_0$ at the integrand's pole, and introducing the shorthand notation $\int a(p) \,{\rm d}p= a_q$ (omitting the subscript 0 on $q$), (3.17) becomes

(3.19)\begin{equation} (K-\langle q_0|\tilde V-\tilde V_w|q_0\rangle)a_q = \sum_j\langle q_0|\tilde{V}|j\rangle a_j. \end{equation}

Treatment of the discrete modes uses again the approximation that over the resonance $\langle j |\tilde V|p\rangle$ can be regarded as independent of $p$; then for $l=4,2$:

(3.20)\begin{align} (\lambda_\perp{-}\lambda_l)a_l=\sum_j\langle l|\tilde{V}|j\rangle a_j +\int\langle l|\tilde{V}|p\rangle a(p)\,{\rm d}p =\sum_j\langle l|\tilde{V}|j\rangle a_j +\langle l|\tilde{V}|q_0\rangle\int a(p)\,{\rm d}p. \end{align}

We regard the amplitudes of the three modes $4,2,q_0$ as composing a column 3-vector whose transpose is

(3.21)\begin{equation} \boldsymbol{\mathsf{a}}^T=(a_4,a_2,a_q)=(a_4,a_2,\int a(p)\,{\rm d}p). \end{equation}

The three relations in (3.19) and (3.20) can be considered a matrix equation $\boldsymbol{\mathsf{M}}\boldsymbol{\mathsf{a}}=0$, where

(3.22) \begin{equation} \boldsymbol{\mathsf{M}}=\left[ \begin{array}{@{}ccc@{}} \lambda_4-\lambda_\perp{+}\langle4|\tilde V|4\rangle & \langle4|\tilde V|2\rangle & \langle4|\tilde V|q_0\rangle\\ \langle2|\tilde V|4\rangle & \lambda_2-\lambda_\perp{+}\langle2|\tilde V|2\rangle & \langle2|\tilde V|q_0\rangle\\ \langle q_0|\tilde V|4\rangle & \langle q_0|\tilde V|2\rangle & -K+\langle q_0|\tilde V-\tilde V_w|q_0\rangle \end{array}\right]. \end{equation}

The complex determinant of ${\mathsf{M}}$ must be zero for a non-trivial solution. To connect more directly with the shiftmode calculation, which simply zeroes the top-left coefficient ${\mathsf{M}}_{11}$, and to avoid large values arising from ${\mathsf{M}}_{33}$, the scaled quantity (denoted $D$) we actually zero is the determinant of ${\mathsf{M}}$ divided by the co-factor of ${\mathsf{M}}_{11}$. If $\omega$ is iterated until $D=0$, the $\omega$ value found will be the corresponding mode frequency, and the mode structure will be given by the solution of $\boldsymbol{\mathsf{M}}\boldsymbol{\mathsf{a}}=0$.

Our reduction of the resonant continuum mode uses just the $m=0,\pm 1$ cyclotron harmonics to determine $q_{0}.$ That is justified by supposing that the finite Larmor radius parameter $\zeta _t=k_\perp \sqrt {T_\perp /m_e}/\varOmega$ is small, and recognizing ${\rm I}_m(\zeta _t^2)\propto \zeta _t^{2m}$. If instead $\zeta _t$ is not small (because $\varOmega \not \gg k_\perp v_{t\perp }$), the dependence ${\rm e}^{-\zeta _t^2}$ predominates in (3.2), and harmonics up to $|m|\sim 3\zeta _t$ approximate a continuous integral over the perpendicular Maxwellian, with the relevant range of $\omega _m$ approximately $\omega \pm 3k_\perp v_{t\perp }$. The resultant frequency spread exceeds ${\sim }\omega$ for the whistler mode (for which $\omega < \varOmega$) implying that high-harmonic (i.e. perpendicular Landau) damping is strong. The continuum mode contribution itself is then suppressed, and it is reasonable to ignore coupling to the continuum modes, reducing $\boldsymbol{\mathsf{M}}$ to its upper-left $2\times 2$ submatrix.

4.1 Internal numerical evaluation of inner products

The evaluation of the overlap integrals (forces) that appear in $\boldsymbol{\mathsf{M}}$ for the inner hole region is carried out numerically using methods that have been documented in detail previously (Hutchinson Reference Hutchinson2018b) for mode $|4\rangle$. In summary the process consists, for each relevant orbit energy $W$, of numerically integrating to obtain the relationship between $z$ and the prior time $\tau$ for the unperturbed orbit, and simultaneously accumulating the integral $\varPhi _m(z)$ for each mode, using a discrete (but non-uniform) $z$ mesh. Trapped orbits $W<0$ require a sum over all prior bounces in the potential well, which is represented by a multiplying bounce-resonant factor. Passing orbits use a single transit across the domain $-z_d< z< z_d$, and for the discrete modes there is no external contribution. The continuum mode at the incoming boundary ($z_{{\rm in}}=-(v/|v|)z_d$), unlike the discrete modes, has a non-zero value, which is provided by the wave expression $\varPhi _m(z_{{\rm in}})=\varPhi _m^{{\rm wave}}$, (3.7). Inside the prior time integration loop, the overlap integrals $\int s^*(z)\tilde f_m(z)\,{\rm d}z$ are simultaneously accumulated. Afterwards they are integrated ${\rm d}v$ over the relevant range of parallel energy $W$, and summed over relevant harmonics $m$ to give the hole-region contributions to $\langle s |\tilde V|u\rangle$. Except when $s$ is the continuum mode (the bottom row of $\boldsymbol{\mathsf{M}}$), these internal values provide the full evaluation of the overlap integrals. For $s=q_0$, however, extra contributions arise from the external region, which we now describe.

4.2 External analytic integration

In the external region $|z|>z_d$, the prior integral $\varPhi$ giving $\tilde f$ for discrete modes $|j\rangle$ is non-zero only for outward-moving orbits. We focus for definiteness on $z>z_d$ and positive $v$. For inward-moving orbits, by contrast, $\tilde f$ is zero externally for the operations $\tilde V|j\rangle$ and $(\tilde V -\tilde V_w)|q_0\rangle$. Consequently, the external contribution to the forces $\langle q _0|\tilde V |j\rangle$ and $\langle q _0|\tilde V -\tilde V_m|q_0\rangle$ arises from the external integration $\int _{z_d}^\infty q_0^*(z) \tilde n \,{\rm d}z$ only for outward orbits, indicated by superscript ‘out’ on $\tilde V$.

To assist with understanding, figure 2 shows a case illustrating the relevant density perturbations $\tilde n_q=\tilde V^{{\rm out}}|q\rangle$, $\tilde n_w=\tilde V_w^{{\rm out}}|q\rangle$ and $\tilde n_4=\tilde V|4\rangle$ arising from positive velocity orbits. The upper panel shows the relevant curves for the continuum mode in the inner and outer regions. The densities $\tilde n_q$ and $\tilde n_w$ are very different in the inner region $|z|\lesssim z_d$ and well outside it but they converge to each other in the wave region $|z|\gg z_l$ ($z_l\sim 20$ for this relatively high-frequency illustration). Also shown is the curve of $\tilde V_w|q\rangle /2$ which is the average of the inward and outward densities. It has a shape similar to $\tilde n_w$ but is somewhat smaller because $\tilde V^{{\rm in}}< \tilde V^{{\rm out}}$. The relevant contribution comes from the difference between $\tilde n_q$ and $\tilde n_w$.

Figure 2.

Illustrating the process of calculating the inner products involving $|q_0\rangle$ and $|4\rangle$. Passing density contributions from positive, outward, velocity orbits integrated over all relevant energies. The square symbol locates the edge of the inner hole region, chosen here as $z_d=15$. The scaled $\varPhi$ shown is for the last included energy $W=6$.

The lower panel shows $\tilde n_4$ which is substantial in the inner region and converges to zero for $|z|\gg z_l$, and one can see residual oscillations caused by discrete contributions at low velocity approximating the ${\rm d}v$ integral, which die out at large $z$. In the upper panel there are also some oscillations but they are barely visible. The ‘$\varPhi$ (scaled)’ curve illustrates (only) the shape of the final contribution to the velocity integral from high velocity, showing how long the wavelength of oscillations becomes there; its contribution to $\tilde n$ is small because $f_{\parallel 0}$ is negligible at high $v$.

Now we explain how the external integrations are performed mostly analytically. In the external region, the velocity is constant; so for $|j\rangle$, which has zero perturbed potential there, $\varPhi$ is a simple time delay factor multiplied by its value at the join between internal and external regions $\varPhi (z_d)$:

(4.1)\begin{equation} \varPhi_{|j\rangle}(z)=\exp({{\rm i}\omega_m(z-z_d)/v}) \varPhi_{|j\rangle}(z_d). \end{equation}

The corresponding expression for the continuum mode, which has non-zero external potential, may be found by substituting the wave expression (3.7) for $|q_0\rangle$ which can be integrated analytically to give an extra term, arriving at

(4.2)\begin{align} \varPhi_{|q_0\rangle}(z)=\exp({{\rm i}\omega_m(z-z_d)/v}) \varPhi_{|q_0\rangle}(z_d) +\frac{A\exp({{\rm i}k_\parallel z_d})}{{\rm i}(k_\parallel v-\omega_m)}[\exp({{\rm i}k_\parallel(z-z_d)}) -\exp({{\rm i}\omega_m(z-z_d)/v})]. \end{align}

Recalling that the wave operator is simply constant, its external contribution to $\varPhi _{|w\rangle }$ for a particular positive velocity is ${A{\rm e}^{{\rm i}k_\parallel z}}/{{\rm i}(k_\parallel v-\omega _m)}$, which exactly cancels the first term in the square bracket of $\varPhi _{|q_0\rangle }(z)$. Such cancellation is essential to produce a finite value for $\langle q _0|\tilde V -\tilde V_w|q_0\rangle$. Hence, when forming $(\tilde V -\tilde V_w)|q_0\rangle$, the required external prior integral is

(4.3)\begin{align} \varPhi_{|q_0-w\rangle}(z)\equiv \varPhi_{|q_0\rangle}(z)-\varPhi_{|w\rangle} =\exp({{\rm i}\omega_m(z-z_d)/v}) \varPhi_{|q_0\rangle}(z_d) -\frac{A\exp({{\rm i}k_\parallel z_d})}{{\rm i}(k_\parallel v-\omega_m)} \exp({{\rm i}\omega_m(z-z_d)/v}). \end{align}

To obtain the total external force we can carry out the inner product $z$-integration analytically as

(4.4)\begin{equation} \int_{z_d}^\infty q_0^*(z) \varPhi_u(z) \,{\rm d}z = \frac{A^*\varPhi_u(z_d)v{\rm e}^{-{\rm i}k_\parallel z_d}}{{\rm i}(k_\parallel v-\omega_m)} + \frac{\delta_{q_0u}|A|^2v}{(k_\parallel v-\omega_m)^2}, \end{equation}

where $\varPhi _u$ refers to the ‘ket’ mode ($j$ or $q_0-w$) and the $|A|^2$ term is present only if we are constructing $\langle q _0|\tilde V -\tilde V_w)|q_0\rangle$, as indicated by $\delta _{q_0u}$. This force quantity is multiplied by the weighting factor $b_m$ (3.3a,b) and integrated (numerically) over parallel velocity so as to produce the total inner product. For asymmetric $f_{\parallel 0}$ the process would need to be carried out also for negative velocity and $z$, but since we take $f_{\parallel 0}$ to be symmetric the integration is carried out only for positive $v$ and $z$; the result is then doubled to account for force exerted at negative $z$. The sum over relevant cyclotron harmonics $m$ is performed last.

In verification of the fairly complex coding we have compared results for mode $|4\rangle$ with prior calculations, confirmed certain analytic self-adjoint properties of the different modes and debugged the $m=0$ results by benchmark comparisons between two independent implementations of the calculation (in Fortran and Python).