3.1. Space, velocity structure of steady-state decay

At the start of each simulation, the plasma relaxes from its initial state over $\mathord {\sim } 1$ – $3 \tau _{\mathrm{bounce}}$ ; the diamagnetic field response is changed, short-wavelength electrostatic fluctuations occur at the plasma edge and plasma escapes from the central cell into the expanders. The plasma reaches a steady-state decay by $t = 6 \tau _{\mathrm{bounce}}$ for all $R_{{m}}$ simulations. At this time: (i) the particle loss time $\tau _{\mathrm{p}} = n / (\mathrm{d} n/\mathrm{d} t)$ is roughly constant and exceeds the ion bounce time ( $\tau _{\mathrm{p}} \gg \tau _{\mathrm{bounce}}$ ); (ii) the plasma beta $\beta _{{i}} = 8\pi P_{{i}}/B^2 \sim 0.1$ to within a factor of two at the origin $(r,z) = (0,0)$ , with $P_{{i}}$ the total ion pressure; (iii) the combined vacuum and diamagnetic fields attain a mirror ratio $R_{{m}} = \{ 21, 45, 69 \}$ somewhat higher than the respective vacuum values $R_{{m}} = \{ 20, 41, 64 \}$ .

Figure 1 shows the plasma’s overall structure at $t = 6 \tau _{\mathrm{bounce}}$ for each of the vacuum $R_{{m}} = 20, 41, 64$ simulations. Flute-like, electrostatic fluctuations at the plasma’s radial edge are visible in $z=0$ slices of ion density and electric fields, with the strongest and most coherent fluctuations for the $R_{{m}} = 20$ case. In figure 1(a,e,i), the axial outflow at $|z| \gtrsim 70 \;\mathrm{cm}$ is split about $r = 0$ , so more plasma escapes from the radial edge $r \gt 0$ than the core $r \sim 0$ . In figure 1(d,h,l), the radial electric field fluctuation $\delta E_r = E_r - \langle E_r \rangle _\theta$ , where $\langle \cdots \rangle _\theta$ represents an average over the azimuthal coordinate to subtract the plasma’s net radial potential. The azimuthal fluctuation $\delta E_\theta$ in figure 1(c,g,k) is defined similarly. The transverse magnetic fluctuations $\delta B_r$ and $\delta B_\theta$ have small amplitudes $\lesssim 10^{-3} B(z=0)$ , whereas the electric fluctuations $\delta E_r$ and $\delta E_\theta$ are of order $0.1 v_{\mathrm{ti0}} B(z=0)/c$ , corresponding to motional flows at thermal speeds. We therefore neglect electromagnetic fluctuations and focus solely on the azimuthal, electrostatic mode visible in figure 1.

Figure 2(a,b) shows ion density profiles along $z$ both on- and off-axis, with horizontal dashes marking the density floor imposed in Ohm’s law, (2.1). The off-axis density is measured along flux surfaces hosting the strong electric fluctuations seen in figure 1. Specifically, we pick surfaces at $r=\{9,18,22\}\;\mathrm{cm}$ and $z=0$ that have an approximate (azimuth-averaged, paraxial) flux coordinate $\psi \approx \int 2\pi \langle B_z \rangle _\theta r \mathrm{d} r = \{ 2.1, 3.8, 3.7\} \times 10^6 \;\mathrm{G\,cm^{2}}$ for the $R_{{m}}=\{20,41,64\}$ simulations, respectively. The plotted density $n_{{i}}$ is also azimuth averaged.

The density profiles peak near the turning points of $45^\circ$ pitch-angle ions, defined as the $z$ locations where the local mirror ratio $R_{{m}}(z)=2$ on axis; i.e. $\{57,44,36\} \;\mathrm{cm}$ (figure 2 c). Comparing on- and off-axis density peaks, the off-axis peak is wider and decreases more slowly towards the mirror throat and expander. This can be explained by the plasma edge’s stronger loss-cone outflow and broader pitch-angle distribution between $0^\circ$ and $45^\circ$ , compared with the plasma core at $r = 0$ (figure 3).

Figure 2(d) shows on- and off-axis electrostatic potential profiles $e\phi (z)/T_{{e}}$ . The off-axis profiles $\phi (s(z)) = -\int E_\parallel (s) \mathrm{d} s$ , with arclength $s$ in the $(r,z)$ plane, are integrated along the same flux surfaces used in figure 2(b). We notice that the $z=0$ potential well has similar depth both on- and off-axis. The density floor truncates the axial electrostatic potential at $z \approx 100$ to $110 \;\mathrm{cm}$ , so the full potential drop from the mirror throat to the domain’s $z$ boundary is not captured in our simulation. In any case, plasma outflow in the expanders is not well modelled by our electron closure, as the outflow is far from thermal equilibrium (e.g. Wetherton et al. Reference Wetherton, Le, Egedal, Forest, Daughton, Stanier and Boldyrev2021). We will restrict our attention to central-cell plasma behaviour that we suppose to be unaffected by the expanders.

Figure 4.

The 3-D rendering of ion density in $R_{{m}}=20$ simulation at $t=6\;\mathrm{\unicode{x03BC}s}$ ; colourmap is ion density in units of $\mathrm{cm}^{-3}$ . An animated movie is available in the online journal.

Figure 3 shows initial ion velocity distributions, as imported into Hybrid-VPIC from CQL3D-m, at the centre of the mirror cell: $z \in (-5.9, 5.9) \;\mathrm{cm}$ for all simulations. Figure 3(a–d) sample ions from the plasma’s radial edge: $r \in [ 5.9, 11.8) \;\mathrm{cm}$ for $R_{{m}} = 20$ ; $r \in [11.8, 23.5) \;\mathrm{cm}$ for $R_{{m}} = 41$ ; $r \in [14.7, 29.4) \;\mathrm{cm}$ for $R_{{m}} = 64$ . Figure 3(e) samples ions from the plasma’s core: $r \in [0,2.9) \;\mathrm{cm}$ for $R_{{m}} = 20$ ; $r \in [0,5.9) \;\mathrm{cm}$ for $R_{{m}} = 41$ ; $r \in [0,7.4) \;\mathrm{cm}$ for $R_{{m}} = 64$ . Figure 3(f–j) shows ion distributions, selected from the same axial and radial regions as figure 3(a–e), after the simulation has reached $t = 6 \tau _{\mathrm{bounce}} \approx 6 \;\mathrm{\unicode{x03BC}s}$ . Ions diffuse mostly in $v_\perp$ ; their distribution is continuous and nearly flat across the velocity-space loss-cone boundary. The reduced distribution $F(v_\perp ) = \int f \mathrm{d} v_\parallel$ has relaxed to a monotonically decreasing shape, $\mathrm{d} F / \mathrm{d} v_\perp \lt 0$ , at the plasma edge (figure 3 f); however, the core plasma maintains $\mathrm{d} F / \mathrm{d} v_\perp \gt 0$ at low $v_\perp$ (figure 3 j). Some distribution function moments will be used in later discussion. We define $\boldsymbol{B}$ -perpendicular and parallel temperatures $T_{\mathrm{i\perp }} \equiv (1/2) \int m_{{i}} v_\perp ^2 f \mathrm{d}\boldsymbol{v}$ and $T_{\mathrm{i\parallel }} \equiv \int m_{{i}} v_\parallel ^2 f \mathrm{d}\boldsymbol{v}$ so that $T_{{i}} = (2T_{\mathrm{i\perp }} + T_{\mathrm{i\parallel }})/3$ ; temperature values for the edge ion distributions at $t=6\,\tau _{\mathrm{bounce}}$ (figure 3 f–i) are given in table 1.

Figure 4 shows a 3-D render of ion density in the $R_{{m}}=20$ simulation at $t=6\;\mathrm{\unicode{x03BC}s}$ . The flute-like ( $k_\parallel \sim 0$ ) nature of the edge fluctuations is apparent. An accompanying movie of the full time evolution from $t=0$ to $6\;\mathrm{\unicode{x03BC}s}$ is available in the online journal.

To summarize, figures 1–4 show that at the plasma’s radial edge: (i) flute-like electrostatic fluctuations appear; (ii) axial outflow and hence losses are enhanced relative to the plasma’s core at $r \sim 0$ ; and (iii) ions diffuse in $v_\perp$ to drive $\mathrm{d} F/\mathrm{d} v_\perp \lt 0$ . It is already natural to suspect that the electrostatic fluctuations diffuse ions into the loss cone and hence cause plasma to escape the mirror.

Figure 5.

Radial structure of plasma at midplane $z=0$ and at $t = 6 \,\tau _{\mathrm{bounce}} \approx 6 \;\mathrm{\unicode{x03BC}s}$ , for simulations with vacuum (a–e) $R_{{m}}=20$ , (f–j) $41$ , (k–o) $64$ . Panels (a–c), (f–h) and (k–l) show azimuth-averaged radial profiles of (a) ion density $n_{{i}}$ , (b) ion density gradient $\epsilon \rho _{\mathrm{i0}}$ , (c) azimuthal electrostatic fluctuation energy $\delta E_\theta ^2$ . Horizontal shaded bars contain the ‘edge’ ion distributions from figure 3. Vertical dashes in (a), (f) and (k) mark density floor for (2.1). Panels (d), (e), (i), (j), and (n), (o) show azimuthal Fourier spectra of density $\tilde {n}_{{i}}(r,m)$ and azimuthal electric field $\tilde {E}_\theta (r,m)$ ; Fourier transform maps $\theta \to m$ , but radius $r$ is not transformed. White rays mark azimuthal wavenumber $k \rho _{\mathrm{i0}} = 2,4,6,8,10,12$ , with $k = m/r$ . Dashed pink ray is the maximum $k = \pi /\Delta r$ resolved by the spatial grid, taking $\Delta r = \sqrt {2} \Delta x$ . Panels (f)–(j) and (k)–(o) are organized similarly.

3.2. Drift cyclotron mode identification

To establish the electrostatic mode’s nature, we need to know plasma properties at the radial edge and the mode’s wavenumber and frequency spectrum.

Figure 5(a–c,f–h,k–l) presents the radial structure of the ion density $n_{{i}}$ , and the electrostatic fluctuation energy $\delta E_\theta ^2 = \langle E_\theta ^2 \rangle _\theta - \langle E_\theta \rangle _\theta ^2$ , at the mirror midplane $z=0$ . Figure 5(d,e,i,j,n,o) also presents Fourier spectra of density $\tilde {n}_{{i}}(m,r)$ and electric component $\tilde {E}_\theta (m,r)$ as a function of azimuthal mode number $m$ and radius $r$ . Beware that Fourier spectrum normalization is arbitrary here and in all figures; Fourier amplitudes may be compared between panels within one figure, but not across distinct figures.

The density gradient $\epsilon \equiv (\mathrm{d} n_{{i}}/\mathrm{d} r)/n_{{i}}$ , in units of inverse ion Larmor radius $\rho _{\mathrm{i0}}^{-1}$ , is of order unity and increases with $R_{{m}}$ (figure 5 b,g,l); equivalently, the plasma column radius is smaller in units of $\rho _{\mathrm{i0}}$ for larger $R_{{m}}$ , despite the column’s larger physical extent.

The mode spectra of $\tilde {n}$ and $\tilde {E}_\theta$ suggest a partial decoupling of density and electric fluctuations (figure 5, d, e, i, j, n, o). In all simulations, low $m \sim 2$ – $4$ density fluctuations are not accompanied by a strong $E_\theta$ signal (figure 5 d, e, i, j, n, o). The $R_{{m}}=20$ simulation shows a strong mode in both density and $E_\theta$ fluctuations at $m \approx 9$ – $10$ and equivalent angular wavenumber $k \rho _{\mathrm{i0}} \approx 2$ – $4$ (figure 5 d,e). We identify this Fourier signal with phase-coherent fluting at the same $m$ visible to the eye in figure 1(b,c). In contrast, the $R_{{m}}=41,64$ simulations show a decoupling of density and $E_\theta$ fluctuations. The strongest density fluctuations reside at $r \sim 1$ – $2\rho _{\mathrm{i0}}$ , $m \sim 7$ – $8$ and $k\rho _{\mathrm{i0}} \approx 2$ – $6$ (figure 5 i,n), whereas the electrostatic fluctuations reside at larger $r \sim 2$ – $4\rho _{\mathrm{i0}}$ , $m \sim 15$ – $30$ and $k \rho _{\mathrm{i0}} \sim 4$ – $12$ (figure 5 j,o).

The fluctuations have $k_\parallel \ll k_\perp$ and are thus flute-like, which we checked by plotting $\boldsymbol{E}$ in approximate flux-surface coordinates (not shown). Electric-field fluctuations terminate at the mirror throats and do not extend into the expanders; fluctuations may be artificially truncated by the density floor in (2.1).

Figure 6.

Time-azimuth Fourier spectra of density $\tilde {n}(\omega ,m)^2$ (a–c) and electric field $\tilde {E}_\theta (\omega ,m)^2$ (d–f) for simulations with $R_{{m}}=\{20, 41, 64\}$ (left to right). Panels (g)–(l) show corresponding $(\omega ,k)$ of unstable DCLC modes predicted by (3.4) for edge $F(v_\perp )$ at $t \approx 6 \;\mathrm{\unicode{x03BC}s}$ (g–i) or $t = 0$ (j–l). In panels (a)–(f), the full $\omega$ range within Nyquist-sampling limits is shown; signals with $\omega \gtrsim 2 \Omega _{\mathrm{i0}}$ alias in frequency. White dotted lines plot ion diamagnetic drift velocity $\omega /k = v_{\mathrm{Di}}$ . Shaded vertical bar in (a), (d) marks grid resolution limit $k \gt \pi /\Delta r$ with $\Delta r = \sqrt {2}\Delta x$ . In (g)–(l), we plot both stable- and unstable-mode frequencies $\textrm {Re}(\omega )$ (black, blue), and also the corresponding unstable-mode growth rates $\textrm {Im}(\omega )$ (green). In (l) only, red curves plot $\textrm {Im}(\omega )$ for higher- $\omega /k$ modes with $\textrm {Re}(\omega ) \in [4\Omega _{\mathrm{i0}},14\Omega _{\mathrm{i0}}]$ beyond the plot extent. Black dotted lines plot $\omega /k = v_{\mathrm{Di}}$ .

Joint time-frequency and azimuthal-mode spectra of density and electric field fluctuations, $\tilde {n}(\omega ,m)$ and $\tilde {E}_\theta (\omega ,m)$ , are presented in figure 6(a–f). Fluctuations are sampled at radii $r = \{3.34, 2.98, 2.69 \} \rho _{\mathrm{i0}}$ , respectively, over $t = 3$ to $6 \,\tau _{\mathrm{bounce}}$ ; $\omega$ is angular frequency. Positive $\omega /k$ corresponds to the ion diamagnetic drift direction. We interpret Fourier power at $\omega \lt 0$ as high- $\omega$ signal that is aliased in frequency space and would otherwise be contiguous in physical $(\omega ,k)$ . Assuming so, both $\tilde {n}$ and $\tilde {E}_\theta$ show a mode spectrum with a phase speed $\omega /k \gt 0$ comparable to the ion diamagnetic drift speed $v_{\mathrm{Di}}$ (white dotted lines, figure 6 a–f). We compute

(3.1) \begin{equation} v_{\mathrm{Di}} \approx \frac {T_{\mathrm{i\perp }}/m_{{i}}}{\Omega _{\mathrm{i0}}} \left (-\frac {1}{n_{{i}}} \frac {\mathrm{d} n_{{i}}}{\mathrm{d} r}\right ) = v_{\mathrm{ti0}} \frac {T_{\mathrm{i\perp }}}{T_{\mathrm{i0}}} |\epsilon | \rho _{\mathrm{i0}} \end{equation}

using $\epsilon \rho _{\mathrm{i0}} = \{ -1, -1, -1.5 \}$ and $T_{\mathrm{i\perp }}$ measured at $t=6\tau _{\mathrm{bounce}}$ (values reported in § 3.1). The spectra align with $\omega /k = v_{\mathrm{Di}}$ within a factor of $2$ .

A fundamental mode appears at $\omega \in [\Omega _{\mathrm{i0}}, 2 \Omega _{\mathrm{i0}}]$ in all simulations. Fluctuation power extends to $\omega \gtrsim 3 \Omega _{\mathrm{i0}}$ in all simulations, perhaps up to $\omega \gtrsim 7 \Omega _{\mathrm{i0}}$ in the $R_{{m}}=64$ simulation, but by eye we do not discern discrete harmonics above $2 \Omega _{\mathrm{i0}}$ . A low-frequency $\omega \ll \Omega _{\mathrm{i0}}$ mode with non-zero $m$ appears chiefly in $\tilde {n}$ and weakly in $\tilde {E}_\theta$ ; we identify this slower motion as fluid interchange and discuss it further in § 4.4.

To help interpret figure 6(a–f), we compute the linearly unstable $(\omega ,k)$ for DCLC in a planar-slab plasma with a spatial density gradient $\epsilon$ and uniform background magnetic field ( ${\nabla } B = 0$ ). In such a plasma, a dispersion relation for exactly perpendicular electrostatic waves can be obtained by integrating over unperturbed orbits and Taylor expanding $f$ in particle guiding-centre coordinate, following Stix (Reference Stix1992, § 14–3, (9)). The dispersion relation is then

(3.2) \begin{equation} D = 1 + \sum _s \chi _s = 0,\end{equation}

where the perpendicular ( $k=k_\perp$ ) susceptibility of species $s$ reads

(3.3) \begin{align} \chi _{\mathrm{s}} = \left (\frac {\omega _{\mathrm{p}s}}{\Omega _{\mathrm{s}}}\right )^2 \bigg [ \left ( 1 - \frac {\epsilon \omega }{k} \right ) \frac {1}{k^2} &\sum _{n=-\infty }^{\infty } \frac {n}{\omega - n} \int \mathrm{d}^3\boldsymbol{v} \left ( \frac {1}{v_\perp } \frac {\partial f}{\partial v_\perp } \right ) J_n^2 \nonumber \\ - \frac {\epsilon }{k} &\sum _{n=-\infty }^{\infty } \frac {1}{\omega - n} \int \mathrm{d}^3\boldsymbol{v} f J_n^2 \bigg ] . \end{align}

In (3.3), variables are written in a species-specific dimensionless form: $\omega /\Omega _{\mathrm{s}} \to \omega$ , $k\rho _s \to k$ , $\epsilon \rho _s \to \epsilon$ and $v_\perp /v_{\mathrm{t}s} \to v_\perp$ , where $\Omega _{\mathrm{s}}$ is signed (i.e. $\Omega _{{e}} \lt 0$ ) and $\rho _s \equiv v_{\mathrm{t}s}/\Omega _{\mathrm{s}}$ . The decision of how to define $v_{\mathrm{t}s}$ (with or without $\sqrt {2}$ ) is given to the user. The plasma frequency $\omega _{\mathrm{p}s} = \sqrt {4\pi n_{\mathrm{s}} q_{\mathrm{s}}^2/m_{\mathrm{s}}}$ for each species. The Bessel functions $J_n = J_n(k_\perp v_\perp )$ as usual, with $k_\perp =k$ . Equations (3.2)–(3.3) simplify for cold fluid electrons to yield

(3.4) \begin{equation} D = 1 + \chi _{{i}} + \frac {\omega _{\mathrm{pe}}^2}{\Omega _{{e}}^2} + \frac {\omega _{\mathrm{pe}}^2}{|\Omega _{{e}}|} \frac {\epsilon }{k \omega } = 0 , \end{equation}

where the variables $k$ , $\epsilon$ and $\omega$ are now in dimensional units. Equation (3.4) is the slab DCLC dispersion relation also used by Lindgren et al. (Reference Lindgren, Birdsall and Langdon1976, (2)), Ferraro et al. (Reference Ferraro, Littlejohn, Sanuki and Fried1987, (19)), and Kotelnikov et al. (Reference Kotelnikov, Chernoshtanov and Prikhodko2017, (17) and (A14)). In our sign convention, $\epsilon \lt 0$ obtains DCLC with $\omega /k \gt 0$ in the ion diamagnetic drift direction. Equation (3.4) also hosts normal modes with $k \lt 0$ and high phase velocity in the electron diamagnetic drift direction (Lindgren et al. Reference Lindgren, Birdsall and Langdon1976, § 2.A.1.b), which do not appear in our simulations and so are omitted from our discussion.

The unstable- and normal-mode solutions to (3.4), presented in figure 6(g–l), are computed as follows. First, we take $\Omega _{\mathrm{i0}}$ , $\rho _{\mathrm{i0}}$ and $v_{\mathrm{ti0}}$ as defined in § 2.3 to normalize all variables in (3.4). Plasma parameters used for the $R_{{m}} = \{ 20,41,64 \}$ simulations, respectively, are $\epsilon \rho _{\mathrm{i0}} = \{ -1, -1, -1.5 \}$ ; $n_{{i}} = \{ 4, 1.2, 0.5 \} \times 10^{12} \;\mathrm{cm^{-3}}$ ; $B = \{ 8.6, 4.1, 2.7 \} \times 10^{3} \;\mathrm{G}$ . Both $\epsilon$ and $n_{{i}}$ describe the plasma edge at the midplane $z=0$ (figure 5). We take $B$ at $(r,z)=(0,0)$ to match the variable normalization throughout this manuscript; $B$ at the plasma edge differs by $\lesssim 10\,\%$ . Reduced ion distributions $F(v_\perp ) = \int f \,\mathrm{d} v_\parallel$ are measured directly from the plasma edge (figure 3). Bessel function sums are computed using all terms with index $|n| \leqslant 40$ . The waves and particles at hand have $k_\perp \rho _{\mathrm{i0}} \lesssim 20$ and $v_\perp /v_{\mathrm{ti0}} \lesssim 2$ , so the Bessel function argument $(k_\perp \rho _{\mathrm{i0}}) (v_\perp /v_{\mathrm{ti0}}) \lesssim 40$ . Terms with $n \gt 40$ contribute little to $\chi _{{i}}$ because the first positive oscillation of $J_n(\xi )$ peaks at $\xi = j_n' \gt n$ , where $j_n'$ is the smallest positive zero of $J_n'$ (Watson Reference Watson1922, §15.3), and $J_n(\xi ) \to 0$ quickly as $\xi \to 0$ for $\xi \lesssim n$ .

We then compute $D$ on a discrete mesh of $(k, \textrm {Re}(\omega ), \textrm {Im}(\omega ))$ ; for each $k$ , we identify normal modes (whether stable, damped or growing) by seeking local minima of $D$ with respect to the complex $\omega$ mesh. Our $D \approx 0$ solutions are not exact. To test our solution scheme, we refined our solutions to $D=0$ by applying a manual root-finder to each $(k,\omega )$ normal mode for one set of plasma parameters, and we saw no significant difference.

Figure 6(g–i) uses $F(v_\perp )$ measured from the plasma edge at $t = 6 \tau _{\mathrm{bounce}} \approx 6 \;\mathrm{\unicode{x03BC}s}$ , showing DCLC modes at marginal instability (more precisely, drift-cyclotron modes since the loss cone is filled).

Figure 6(j–l) uses $F(v_\perp )$ measured at $t=0$ instead to show that initial distributions with empty loss-cones and spatial gradient $\epsilon \rho _{\mathrm{i0}} \sim \mathcal{O}(1)$ drive strongly unstable, broad-band electrostatic modes with fastest growth towards high $k\rho _{\mathrm{i0}} \gg 1$ and $\omega \gg \Omega _{\mathrm{i0}}$ . The $R_{{m}}=64$ simulation (figure 6 l) is an exception, because its CQL3D-m model predicts a larger population of trapped cool ions that helps stabilize DCLC. Figure 6(l) also reveals three branches of unstable modes, each with distinct $\omega /k$ , that we speculate may be drift waves associated with distinct hot and cool plasma populations (figure 3 a,d). The slowest branch is visible with $\textrm {Re}(\omega )$ between $2$ to $4\Omega _{\mathrm{i0}}$ ; the corresponding $\textrm {Im}(\omega )$ are plotted in green. The faster phase speed branches have unstable $\textrm {Re}(\omega ) \gt 4 \Omega _{\mathrm{i0}}$ extending to at least $14\Omega _{\mathrm{i0}}$ ; the corresponding $\textrm {Im}(\omega )$ are plotted in light red.

What is learned from comparing the simulation spectra versus linear theory in figure 6? First, marginally stable DCLC mode growth may explain high $k\rho _{\mathrm{i0}} \gtrsim 5$ fluctuations residing in the device during steady-state decay. How do we explain the fundamental mode between $\omega =\Omega _{\mathrm{i0}}$ and $2\Omega _{\mathrm{i0}}$ for simulations with $R_{{m}}=41,64$ , since that mode is predicted to be linearly stable at late times? It may be an initially excited mode that did not damp and so persists to late times; this appears possible for the $R_{{m}}=41$ simulation, where the fundamental is unstable at $t=0$ . Or, it may be excited by nonlinear flow of wave energy from unstable to stable modes; such an explanation may be needed for the $R_{{m}}=64$ simulation, in which cool plasma at the radial edge should quench DCLC growth of the fundamental mode at both $t=0$ and $t=6\;\mathrm{\unicode{x03BC}s}$ . We have interpreted the $t=0$ and $t=6\;\mathrm{\unicode{x03BC}s}$ as most- and least-unstable scenarios for DCLC growth, but the plasma may also transition through other states that destabilize the fundamental mode.

3.3. Ion scattering

To establish a causal link between $\delta E$ fluctuations and axial ion losses, we quantify ion scattering in the $R_{{m}}=20$ simulation as follows. We measure velocity jumps over a short time interval $\delta t \equiv t_1 - t_0 = 0.25 \Omega _{\mathrm{i0}}^{-1}$ for $\mathcal{O}(10^7)$ PIC macroparticles sampled from $|z| \in [0, 5.9] \;\mathrm{cm}$ . Our approach is similar to many other PIC simulation studies; see Yerger et al. (Reference Yerger, Kunz, Bott and Spitkovsky2025) for a recent discussion of nuances in constructing and interpreting such velocity jump moments. Figure 7(a) shows the probability distribution of the $\boldsymbol{B}$ -perpendicular and parallel velocity jumps, $\delta v_\perp = v_\perp (t_1) - v_\perp (t_0)$ and $\delta v_\parallel = v_\parallel (t_1) - v_\parallel (t_0)$ . The distributions are not Gaussian and have long tails. The perpendicular jumps $\delta v_\perp$ are much larger than $\delta v_\parallel$ , as expected for flute-like ( $k_\perp \gg k_\parallel$ ) electrostatic modes and as evident in figure 3.

Figure 7.

Ion scattering measured in $R_{{m}}=20$ simulation, at midplane $z\in [-5.9,5.9] \;\mathrm{cm}$ unless said otherwise. All diffusion coefficients are normalized to $v_{\mathrm{ti0}}^2\Omega _{\mathrm{i0}}$ . (a) Probability distribution of ion velocity jumps, normalized to $v_\perp (t_1)$ and $v_\parallel (t_1)$ , for particles at all radii. (b) Radial profile of ion diffusion $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle / \delta t$ (solid black) compared with $\langle \delta v_\perp \delta v_\perp \rangle / \delta t$ (dotted blue). (c) Predicted radial profile of ion diffusion due to fluctuating fields $\delta E_\theta ^2$ (dotted blue), $\delta E_r^2$ (thin solid blue), and $\delta E_\perp ^2 = \delta E_\theta ^2 + \delta E_r^2$ (thick solid blue), compared with $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle / \delta t$ (black). (d) Numerical convergence in particles per cell for radial profile of $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle / \delta t$ . (e) Effect of measurement time $\delta t$ upon radial profile of $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle / \delta t$ . (f) Effect of measurement time $\delta t$ upon diffusion measured at the midplane $(r,z) \approx (8.2, 0) \;\mathrm{cm}$ (blue curve), and near the beam-ion turning point at $(r,z) \approx (6.8, 50) \;\mathrm{cm}$ (orange curve). (g) The 2-D map of diffusion $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle / \delta t$ computed in discrete $(r,z)$ bins (pixels); only bins with $\gt 100$ particles are shown. Light blue and orange boxes mark measurement locations used in (f).

Ion velocities may jump due to both adiabatic and non-adiabatic motion. To separate these motions, introduce an energy $\mathcal{E} = m_{{i}} v^2/2 + e \langle \phi \rangle _{\theta ,t}$ , where $\langle \cdots \rangle _{\theta ,t}$ is an average over both azimuth angle $\theta$ and time from $t_0$ to $t_1$ . We expect $\mathcal{E}$ to be conserved by particles gyrating in slowly varying $\boldsymbol{E}$ and $\boldsymbol{B}$ fields, at lowest order in a Larmor-radius expansion. Therefore, we attribute jumps in $\mathcal{E}$ to a non-adiabatic kick in perpendicular velocity that we call $\delta v_{\perp \mathcal{E}}$ . We use $\delta \mathcal{E} = \mathcal{E}(t_1) - \mathcal{E}(t_0) = m v_\perp \delta v_{\perp \mathcal{E}} - m (\delta v_{\perp \mathcal{E}})^2/2$ to compute

(3.5) \begin{equation} \frac {\delta v_{\perp \mathcal{E}}}{v_\perp (t_1)} = 1 - \sqrt {1 - \frac {2 \delta \mathcal{E}}{m v_\perp ^2(t_1)}} .\end{equation}

Equation (3.5) requires $\delta \mathcal{E}$ to not exceed the particle’s final perpendicular energy:

(3.6) \begin{equation} \delta \mathcal{E} \lt m v_\perp ^2(t_1)/2 .\end{equation}

Figure 7(a) shows that the probability distribution of $\delta v_{\perp \mathcal{E}}$ , computed only for those particles satisfying (3.6), is marginally narrower than that of $\delta v_\perp$ , as expected if non-adiabatic kicks are the main contribution to $\delta v_\perp$ .Footnote ⁵

The ion diffusion $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle /\delta t$ as a function of radius is shown in figure 7(b); its value is normalized to $v_{\mathrm{ti0}}^2 \Omega _{\mathrm{i0}}$ in all of figure 7(b–g). Here, $\langle \cdots \rangle$ is a velocity-distribution moment computed in radial bins. The use of $\delta v_{\perp \mathcal{E}}$ decreases the measured diffusion as compared with $\langle \delta v_\perp \delta v_\perp \rangle /\delta t$ , as expected.

The ion diffusion due to fluctuating fields $\delta E_\perp (\boldsymbol{r})$ can be described by a diffusion coefficient similar to those used in quasilinear models,

(3.7) \begin{equation} D_{\perp \perp } = \frac {1}{2} \left ( \frac {e}{m_{{i}}} \delta E_\perp \right )^2 \tau _{\mathrm{c}} ,\end{equation}

where $\tau _{\mathrm{c}}$ is a yet-unknown wave–particle correlation time. Equation (3.7) assumes (i) weak but coherent kicks $\delta v_\perp \approx (e/m_{{i}}) \delta E_\perp \tau _{\mathrm{c}}$ ; (ii) a uniform random distribution of angles between $\boldsymbol{v}_\perp$ and $\delta \boldsymbol{E}_\perp$ to obtain a factor of $1/2$ accounting for kicks in gyrophase instead of $v_\perp$ magnitude. For a scattering-measurement time $\delta t \lt \tau _{\mathrm{c}}$ , we expect

(3.8) \begin{equation} D_{\perp \perp } \approx \frac { \langle \delta v_{\perp \mathcal{E}} \, \delta v_{\perp \mathcal{E}} \rangle }{ \delta t } ,\end{equation}

also replacing $\tau _{\mathrm{c}} \to \delta t$ in $D_{\perp \perp }$ .

Choosing $\delta t \lt \tau _{\mathrm{c}}$ is unusual for studies of particle diffusion, as the resulting (3.8) describes a more ‘ballistic’ than diffusive process. But, a short $\delta t$ helps us. When using a longer $\delta t \gg \tau _{\mathrm{c}}$ , at least two issues arise. First, ions gyrate in and out of the scattering zone, as the zone’s radial width is similar to an ion Larmor radius. A typical ion may get one or a few kicks, gyrate out of the scattering zone and drift adiabatically, re-enter the scattering zone to be kicked again, and so on, resulting in a random walk with intermittent large time gaps. The scattering zone’s finite radial width may also introduce bias in the correlation time $\tau _{\mathrm{c}}$ , because a typical inboard (small $r$ ) ion gyrating in and out of the scattering zone sees a redshift $\omega - k_\perp v_\perp$ , whereas a typical outboard (large $r$ ) ion instead sees a blueshift $\omega + k_\perp v_\perp$ . Second, a longer $\delta t$ needed to sample multiple gyration periods $2\pi \Omega _{\mathrm{i0}}^{-1}$ will introduce axial bounce effects. In the $R_{{m}}=20$ simulation, $\tau _{\mathrm{bounce}} \approx 40 \Omega _{\mathrm{i0}}^{-1}$ , and even fewer ion gyrations are executed within $\tau _{\mathrm{bounce}}$ for the higher $R_{{m}}$ cases.

In figure 7(c) we compare (3.8) with the ion diffusion measured from individual particles. The fluctuating energy density is azimuth averaged as $\delta E_r^2 = \langle E_r^2 \rangle _\theta - {\langle E_r \rangle _\theta }^2$ , and similarly for $\delta E_\theta ^2$ ; the sum $\delta E_\perp ^2 = \delta E_\theta ^2 + \delta E_r^2$ . The diffusion due to $\delta E_\theta ^2$ agrees especially well with the particle measurement, whereas the diffusion due to $\delta E_r^2$ agrees less well.

Numerical noise might drive axial losses from the plasma edge in the same way that we are attributing to DCLC, because PIC particle count decreases at the plasma edge. To check this possibility, figure 7(d) shows that the measured ion scattering is converged in the number of particles per cell used. We are confident that ion scattering is not due to numerical noise because (i) the DCLC electric fields have much larger energy density than numerical noise at the grid scale, and figure 7(c) shows good agreement in radial profiles of electric fields and scattering, (ii) we see weak to no $N_{\mathrm{ppc}}$ dependence of scattering rates, whereas if scattering were due to noise, we might expect either an outwards shift in $r$ as $N_{\mathrm{ppc}}$ increases (for fixed DCLC amplitude), or a decrease in scattering rate if noise suppresses DCLC amplitude; (iii) ion scattering is clearly anisotropic (figure 7 a), whereas numerical scattering should be insensitive to $v_\parallel$ versus $v_\perp$ because the grid scale is much smaller than the ion Larmor radius.

In figure 7(e) we show the effect of $\delta t$ upon the radial profiles of measured ion diffusion. Figure 7(f) then samples the ion scattering at its radial peak $r = 8.2 \;\mathrm{cm}$ (blue curve) and shows its dependence upon many more values of $\delta t$ . We see that the diffusion moment scales linearly with small $\delta t$ as expected from (3.8); for comparison, the black dotted line shows an exactly linear correlation with $\delta t$ . As $\delta t$ becomes $\gtrsim \Omega _{\mathrm{i0}}^{-1}$ , waves and particles decorrelate and the diffusion rate begins to fall. We perform a similar calculation at $(r,z) \approx (6.8,50) \;\mathrm{cm}$ (figure 7 f, orange curve) to conclude that $\tau _{\mathrm{c}}$ is shorter near the fast ion turning point. If $\tau _{\mathrm{c}} \sim 1/\Omega _{{i}}$ (where $\Omega _{{i}}$ varies with $z$ , unlike $\Omega _{\mathrm{i0}}$ ), the lower $\tau _{\mathrm{c}}$ can be easily explained by the $2\times$ increase in $B$ magnitude.

Finally, figure 7(g) shows $\langle \delta v_{\perp \mathcal{E}} \delta v_{\perp \mathcal{E}} \rangle /\delta t$ as a function of $(r,z)$ in the mirror’s central cell. The ion scattering at all $z$ is well localized to the same flux surfaces between beam-ion turning points. Scattering is strongest towards $z=0$ , where the central-cell field is relatively uniform.

We conclude from figure 7(e,f) that particle scattering has a longer correlation time $\tau _{\mathrm{c}}$ and reaches a larger amplitude at the mirror midplane $z=0$ , as compared with near the beam-ion turning points. Ions at $z=0$ , and throughout the central cell where $B \approx B(z=0)$ , should be more important for regulating DCLC growth and saturation than ions at the turning points.

3.4. Particle confinement time

Because the loss cone is full – i.e. $F(v_\perp )$ is roughly constant within the loss cone (figure 3) – our simulated mirrors are a collisionless analogue of the GDT at the Budker Institute (Ivanov & Prikhodko Reference Ivanov and Prikhodko2017). Ions scatter across the loss-cone boundary as fast as (or faster than) untrapped ions can stream out of the mirror, implying an effective mean free path shorter than the device’s length. The particle confinement time $\tau _{\mathrm{p}} \equiv N / |\mathrm{d} N/ \mathrm{d} t|$ , where $N$ is the total number of ions, then scales like the eponymous ‘gas dynamic’ time,

(3.9) \begin{equation} \tau _{\mathrm{GD}} = \frac { 2 L_{\mathrm{p}} R_{{m}} } { v_{\mathrm{ti}\parallel } } ,\end{equation}

adapted from Endrizzi et al. (Reference Endrizzi2023, § 3) with $v_{\mathrm{ti}\parallel }$ a characteristic parallel thermal velocity.

Figure 8.

(a) Particle confinement time measured between $t=5$ to $6\tau _{\mathrm{bounce}}$ , for mirrors of varying $R_{{m}}$ (blue, orange, green) and device length $L_{\mathrm{p}}$ (circle, triangle, star markers) as a function of $\tau _{\mathrm{GD}}$ (3.9). Small blue markers vary $T_{{e}}$ for $R_{{m}}=20$ ; large blue marker is fiducial $T_{{e}}=1.25\;\mathrm{keV}$ . (b) Diffusion time scale $1/\nu _{\perp \perp }$ (3.10) modelled from $\delta E_\theta ^2$ (solid markers) and $\delta E_\perp ^2$ (hollow markers), as a function of $\tau _{\mathrm{GD}}$ . In both panels, diagonal dotted line is $\tau _{\mathrm{p}} = \tau _{\mathrm{GD}}$ .

To test the relation $\tau _{\mathrm{p}} \propto \tau _{\mathrm{GD}}$ , we measure $\tau _{\mathrm{p}}$ between $t = 5$ to $6\, \tau _{\mathrm{bounce}}$ , and $v_{\mathrm{ti}\parallel } = \langle v_\parallel ^2 \rangle ^{1/2}$ at $t=6\,\tau _{\mathrm{bounce}}$ , in each of the $R_{{m}}=\{20,41,64\}$ simulations with $L_{\mathrm{p}} = 98 \;\mathrm{cm}$ on axis. We also measure $\tau _{\mathrm{p}}$ and $v_{\mathrm{ti}\parallel }$ in additional $R_{{m}} = 20$ simulations with varying $T_{{e}} = {0,2.5,5,10}$ keV and longer central cells (larger $L_{\mathrm{p}}$ ); the latter are constructed as follows. Split the ‘original’ mirror device in half at $z=0$ . Between the mirror halves, insert a cylindrical plasma of length $98$ or $168$ cm, thereby increasing the entire mirror’s half-length $L_{\mathrm{p}}$ by $1.5$ or $2\times$ . The cylinder has, at all $z$ , the same velocity distribution and magnetic field $\boldsymbol{B}$ as in the original mirror at $z=0$ . The simulation domain is made larger; mesh voxel dimensions ( $\Delta x$ , $\Delta y$ , $\Delta z$ ) are the same as in § 2. The cylinder’s magnetic field is unphysical because it has $\mathrm{d} B_z/\mathrm{d} r \ne 0$ and $B_r = 0$ , implying non-zero current $c {\nabla }\times \boldsymbol{B}/(4\pi )$ , so we exclude this current from the $\boldsymbol{j} \times \boldsymbol{B}$ term in Ohm’s law(2.1).

The confinement time $\tau _{\mathrm{p}} \sim \mathcal{O}(10^2) \;\mathrm{\unicode{x03BC}s}$ , and $\tau _{\mathrm{p}}$ scales linearly with $\tau _{\mathrm{GD}}$ as expected (figure 8 a). Gas-dynamic confinement explains losses from the $R_{{m}} = 20$ and $41$ simulations very well. Raising electron temperature $T_{{e}}$ from $0$ to $10 \;\mathrm{keV}$ lowers $\tau _{\mathrm{p}}$ from $57$ to $35 \;\mathrm{\unicode{x03BC}s}$ for the $R_{{m}} = 20$ simulations. For comparison, the collisional (aka ‘classical’) confinement time is $0.1$ – $0.2 \;\mathrm{s}$ , using (3.4) of Endrizzi et al. (Reference Endrizzi2023) with $n=3\times 10^{13}\;\mathrm{cm^{-3}}$ , $R_{{m}} = 20$ to $64$ , and beam energy $25\;\mathrm{keV}$ .

The $R_{{m}} = 64$ shows 20 % better particle confinement than predicted by (3.9). Why? The larger plasma radius and hence longer flux-tube length $\gt 2 L_{\mathrm{p}}$ between mirror throats only explains $\mathord {\sim } 5\,\%$ of the disagreement. We speculate that electrostatic potential effects may explain the remaining disagreement. In the $R_{{m}} =20,41$ cases, beam ions diffuse in $v_\perp$ and escape with high $v_\parallel$ ; electrostatic effects are weak since $T_{{e}} \ll T_{{i}}$ , so (3.9) accurately describes the beam ion confinement. The $R_{{m}} = 64$ case has more cool, low-temperature ions (figure 3 i) that can be trapped by the sloshing-ions’ potential peak at $z \sim 30 \;\mathrm{cm}$ (Kesner Reference Kesner1973, Reference Kesner1980); confinement is thus improved.

Instabilities in many settings are self-regulating; i.e. unstable waves drive phase-space flow that quenches the waves’ own energy source, driving the system to equilibrium (e.g. Kennel & Petschek Reference Kennel and Petschek1966). If DCLC self-regulates, then we may expect its amplitude to grow in time until the diffusion rate into the loss cone balances the axial outflow rate: $\nu _{\perp \perp } \propto \tau _{\mathrm{GD}}^{-1}$ . We test this by computing a diffusion rate into the loss cone as

(3.10) \begin{equation} \nu _{\perp \perp } \equiv \frac {1}{N_{\mathrm{\ell }}} \int \frac { D_{\perp \perp }(r) }{ {v_{\perp \mathrm{LC}}}^2 } n_{{i}}(r) 2\pi r \ell \,\mathrm{d} r , \end{equation}

which is a density-weighted average of $D_{\perp \perp }$ (3.7) over a cylindrical kernel of axial length $\ell$ and radial profile $n_{{i}}(r)$ , normalized to $N_{\mathrm{\ell }} = \int n_{{i}}(r) 2\pi r \ell \mathrm{d} r$ . We take $v_{\perp \mathrm{LC}} = v_{\mathrm{ti}\parallel }/\sqrt {R_{{m}} - 1}$ to approximate the ions’ $v_\perp$ at the loss cone boundary, we take $\ell = 12 \;\mathrm{cm}$ centred at $z=0$ , and we take $\tau _{\mathrm{c}} = \Omega _{\mathrm{i0}}^{-1}$ . Given these assumptions, and given that (3.7) is not from a self-consistent quasilinear theory, we interpret (3.10) as no more accurate than an order-of-magnitude scaling. In figure 8(b) we compute the diffusion time scale $1/\nu _{\perp \perp }$ using either $\delta E_\theta ^2$ or $\delta E_\perp ^2$ as defined as in figure 7(c). We observe that $1/\nu _{\perp \perp }$ has similar magnitude as $\tau _{\mathrm{GD}}$ , as expected. But, no trend is obvious from the scatter and few data points.