2. Bounce-averaged theory in COM coordinates

Here, we briefly review the basics of BAFP theory. We move quickly, with a primary focus on magnetic mirror geometries.

In a magnetic mirror, whether simple, tandem or centrifugal, particles in the mirror live on trapped trajectories. On the shortest time scale, the gyroperiod $1/\varOmega _c$ , the particle gyrates around the magnetic field line. On the next longest time scale, the bounce period $\tau _b \sim L_\parallel / v_\parallel$ , the particle transits the mirror, completing a closed orbit.

Absent instabilities, particle deconfinement in such mirrors occurs due to collisional processes. Thus, on a collision time scale $\tau _c \gg \tau _b, 1/\varOmega _c$ , collisions push the particle onto new orbits. Eventually, the particle is pushed onto an orbit that leaves the device.

Obviously, it is inefficient (whether analytically or computationally) to resolve the gyroperiod, bounce time scales and complex trajectories when evaluating the slower diffusion processes that drive particles out of the device. Suitably averaging over these fast time scales to produce a lower-dimensional set of equations describing these slow processes is the goal of BAFP.

The first step in deriving the theory is to transform to coordinates that are constant along a trajectory. For an axisymmetric, non-relativistic plasma, a suitable set of COMs are:

(2.1) $$\begin{align}\mu &\equiv \frac {1}{2} \frac {m v_\perp ^2}{B(\boldsymbol{x})},\end{align}$$

(2.2) $$\begin{align}\epsilon &\equiv \frac {1}{2} m v_\parallel ^2 + \mu B(\boldsymbol{x}) + \psi (\boldsymbol{x}),\end{align}$$

(2.3) $$\begin{align}p_\theta &\equiv m \varOmega _r r + q \varPhi (\boldsymbol{x}).\end{align}$$

Here, $\mu$ is the magnetic moment of the particle, $\epsilon$ is the energy and $p_\theta$ is the azimuthal momentum. Additionally, $v_\parallel$ and $v_\perp$ are the velocity components parallel and perpendicular to the magnetic field, $\varOmega _r(\varPhi )$ is the rotation frequency of the magnetic surface, and $q$ and $m$ are the charge and mass of the particle. The potential energy term $\psi$ is given by

(2.4) \begin{align} \psi = q \phi _\parallel + \frac {1}{2} m r^2 \varOmega _r^2 ,\end{align}

where $\phi _\parallel$ is a potential coordinate along the flux surface (Kolmes et al. Reference Kolmes, Ochs, Rax and Fisch2024). Finally, $\varPhi$ is the flux function, given by

(2.5) \begin{align} \varPhi = r A_\theta . \end{align}

In looking at $p_\theta$ , note that the ratio of the first term to the second is given by

(2.6) \begin{align} \frac {m \varOmega _r }{q A_\theta } \sim \frac {\varOmega _r}{\varOmega _c}, \end{align}

i.e. $p_\theta$ conservation implies that particles stay on their flux surface for particles far from the Brillouin limit (Rax et al. Reference Rax, Fruchtman, Gueroult and Fisch2015). We assume in the following that we are in this regime, taking $\varPhi$ as the third COM in place of $p_\theta$ .

To reach our six dimensions of phase space, we need three other variables. These we take to be the azimuthal angle $\theta$ , the gyro-angle $\alpha$ and the distance along a field line $s$ . These variables $(\theta ,\alpha ,s)$ are the variables we will integrate out to form the bounce-averaged theory. Integrating out $\alpha$ is the basis for gyro or drift kinetics; integrating out $\theta$ makes the theory axisymmetric and integrating out $s$ results in the bounce average.

Now, in general, we will start with an equation of the form:

(2.7) \begin{align} \sqrt {g_X} \frac {\partial f}{\partial t} = \frac {\partial }{\partial X^i} \left [ \sqrt {g_X} \, \Gamma _X^i (f,\boldsymbol{X})\right ] + \sqrt {g_X} S(\boldsymbol{X}) , \end{align}

where $\Gamma ^i(f,\boldsymbol{X})$ is a flux operator, $S(\boldsymbol{X})$ is a source operator and $\sqrt {g_X}$ is the volume element in the space $\boldsymbol{X} \equiv (\boldsymbol{x},\boldsymbol{p})$ . This equation is in conservation form; integrating the left-hand side gives the total change in particle number in a given region, while the first term on the right-hand side reduces to a surface integral of the flux out of the region.

To get our bounce-averaged Fokker–Planck equation, we first need to convert to the six-dimensional (6-D) space $\boldsymbol{Y}$ that contains the COMs, i.e. $\boldsymbol{Y} \equiv (\epsilon ,\mu ,\varPhi ,\theta ,\alpha ,s)$ . Because the metric in phase space $\boldsymbol{X} \equiv (\boldsymbol{x},\boldsymbol{v})$ is $g_{ij} = \delta _{ij}$ , corresponding to a unit volume element, the volume element $\boldsymbol{\mathbf{d}V} =\sqrt {g_Y} \boldsymbol{\mathbf{d}Y}$ in the new coordinates is

(2.8) \begin{align} \sqrt {g_Y} = 2\sqrt {\left | \frac {\partial X^m}{\partial Y^i} \frac {\partial X^n}{\partial Y^j} \delta _{ij} \right | } = \frac {\sqrt {2}}{m^{3/2}\sqrt {\epsilon - \mu B - \psi }}. \end{align}

Here, the factor of 2 comes from the fact that positive and negative $v_\parallel$ are condensed into the same portion of COM space. We also convert the flux $\Gamma$ to $Y$ space. For instance, for a Fokker–Planck flux operator:

(2.9) \begin{align} \Gamma ^i_X = A_X^i f + D_X^{ij} \frac {\partial f}{\partial X^j}, \end{align}

we take

(2.10) \begin{align} \Gamma ^i_Y = \frac {\partial Y^i}{\partial X^m} A_X^m f + \frac {\partial Y^i}{\partial X^m} \frac {\partial Y^j}{\partial X^n}D_X^{mn} \frac {\partial f}{\partial Y^j}, \end{align}

Now, we assume that $f$ is independent of $(\theta ,\alpha ,s)$ , i.e. take $f(Z)$ , $\boldsymbol{Z} \equiv (\epsilon ,\mu ,\varPhi )$ . Then, we integrate the Fokker-Planck equation over the non-COM coordinates in $\boldsymbol{Y}$ . The integrals over $\theta$ and $\alpha$ are trivial, leading to two factors of $2\pi$ . The integral over $s$ is more non-trivial, but results in

(2.11) \begin{align} \sqrt {g_Z} \frac {\partial f}{\partial t} = \frac {\partial }{\partial Z^i} \left ( \sqrt {g_Z} \, \Gamma ^i (f,\boldsymbol{Z})\right ) + \sqrt {g_Z} S_Z(\boldsymbol{Z}) , \end{align}

where

(2.12) $$\begin{align}\sqrt {g_Z} &\equiv 4 \pi ^2 \int _{s_1}^{s_2}\, {\rm d}s \sqrt {g_Y},\end{align}$$

(2.13) $$\begin{align}\Gamma _Z^i &\equiv \frac {1}{\sqrt {g_Z}} \int _{s_1}^{s_2}\, {\rm d}s \sqrt {g_Y} \, \Gamma ^i_Y ,\end{align}$$

(2.14) $$\begin{align}S_Z &\equiv \frac {1}{\sqrt {g_Z}} \int _{s_1}^{s_2}\, {\rm d}s \sqrt {g_Y} S_Y,\end{align}$$

and where $s_1$ and $s_2$ are the two outer limits of the orbit.

Noting that $v_\parallel (s) = \sqrt {2 (\epsilon - \mu B(s) - \psi (s) )/ m }$ , we see that the volume element $\sqrt {g_Z}$ is proportional to the total bounce time $\tau _b = \int _{s_1}^{s_2} \,{\rm d}s/v_\parallel (s)$ , while the operators $\Gamma _Z^i$ and $S_Z^i$ are averaged based on the time ${\rm ds}/v_\parallel (s)$ they spend in each region. Thus, the phase-space volume average and temporal bounce average are equivalent.

3. Multiple wells and the need for multiple domains

Bounce-averaged theory provides a massive simplification to the collisional Vlasov equations. However, it relies on one major assumption: that $f$ can be written as a function of the COM coordinates alone: $f = f(\boldsymbol{Z})$ . As we will see, this condition is often violated in scenarios relevant for modern mirrors. When it is violated, a single function $f(\boldsymbol{Z})$ on COM space is no longer adequate to describe the complete system. Then, the problem must be solved on a set of bounded domains, with appropriate boundary conditions enforced on the shared boundaries between domains.

To explore these issues, we will make two simplifications to the above mentioned theory. First, we will look at equations for a single value of $\varPhi$ , looking only at the 2-D COM space of $(\epsilon ,\mu )$ . Despite this reduction in dimension, we note that the methods presented here generalise straightforwardly to 3-D COM space.

Second, we will consider a piecewise-constant series of potentials $\psi (n(s))$ and magnetic fields $B(n(s))$ , with the $n$ th piece of width $\Delta s_n$ , so that the integrals in the bounce-averaged theory become sums. Of course, this is equivalent in the limit of $\Delta s \rightarrow 0$ and the number of partitions going to infinity, but working with a small discrete number of areas will allow us to develop the discussion much more clearly.

Consider, then, particle dynamics along a field line. As is often conventional, we will take the mirror dynamics to be symmetric about the mirror midplane, and so consider a sequence of $\psi (n)$ , $B(n)$ , with $n \in \{0,\ldots ,N-1\}$ , and with $n=0$ corresponding to the midplane and $n = N-1$ corresponding to the mirror throat. We normalise the potentials and fields to their midplane values, with $\psi (0) = 0$ and $B(0) = 1$ , and arbitrary units for $\psi$ and $\epsilon$ .

Given sequences $\psi (n)$ and $B(n)$ , we can determine the particles allowed in any given field line segment $n$ . For a particle to be allowed in a region, it must have positive parallel kinetic energy, $v_\parallel ^2 \gt 0$ . From (2.2), we see that this requires

(3.1) \begin{align} \epsilon \geqslant \mu B(n) + \psi (n) . \end{align}

Each pair $(\psi (n)$ , $B(n))$ thus determines an allowed region of $(\epsilon ,\mu )$ space, i.e. the region satisfying (3.1).

Figure 2. (a) Discretised magnetic field $B$ and potential energy $\psi$ as a function of axial segment $n$ and (b) COM-space acessibility plot for electrons in a magnetic mirror with a (typical) outward-pointing electric field. In the COM-space accessibility plot, each line $n\in \{0,1,2\}$ represents the boundary below which particles do not have enough kinetic energy to enter that axial segment.

The simplest case is provided by electrons in a simple mirror (which also applies to both electrons and ions in a centrifugal mirror). An ambipolar potential generally forms (Pastukhov Reference Pastukhov1974; Najmabadi, Conn & Cohen Reference Najmabadi, Conn and Cohen1984; Ochs, Munirov & Fisch Reference Ochs, Munirov and Fisch2023), confining electrons and repelling ions. Thus, for electrons, both $\psi$ and $B$ increase as a function of line segment index $n$ , and so particles are best trapped at the midplane (figure 2). This is the ideal situation for BAFP theory: trajectories that can access segment $n=1$ also access segment $n=0$ , so the distribution function $f(\epsilon ,\mu )$ is clearly single-valued. To apply a BAFP theory, operator averages are performed over segment $n = 0$ for the particles above line 0 and below line 1, and over segments $n\in \{0,1\}$ for the particles that are above line 1 and below line 2. A similar situation will occur whenever $\partial \psi /\partial s \geqslant 0$ and $\partial B/\partial s \geqslant 0$ for the entire field line from the midplane to throat. Such a plasma can even be modelled using midplane momentum coordinates, since every allowed trapped trajectory reaches the midplane ( $n=0$ ).

Figure 3. (a) Discretised magnetic field $B$ and potential energy $\psi$ as a function of axial segment $n$ and (b) COM-space acessibility plot for ions in a magnetic mirror with a (typical) outward-pointing electric field. Because of the decreasing electric potential towards the edge of the device, some ions get trapped between the mirror throat and the midplane, i.e. at $n=1$ . These ions are referred to as ‘Yushmanov-trapped’.

A somewhat more complex case is provided by ions in a simple mirror. Because the potential is repulsive, $\psi$ now decreases as a function of $n$ (figure 3). Thus, a class of ions can develop which are trapped off the midplane in the potential well, despite the higher magnetic field there. These are known as ‘Yushmanov-trapped’ particles (Yushmanov Reference Yushmanov1966; Post Reference Post1987). In a plasma with Yushmanov-trapped particles, a midplane-momentum coordinate theory becomes insufficient, because these particles exist in the plasma but do not reach the midplane. However, the distribution function $f(\boldsymbol{Z})$ is still clearly single-valued, so BAFP theory can still be used, solved on a single domain in COM space. The operator averages are performed over $n=1$ for the Yushmanov particles, over $n=0$ for the particles above line 0 and below line 1, and over $n\in \{0,1\}$ for the passing particles that traverse both segments.

Figure 4. (a) Discretised magnetic field $B$ and potential energy $\psi$ as a function of axial segment $n$ and (b) COM-space acessibility plot for a scenario with a constant magnetic field and an internal potential maximum. This scenario exhibits a bifurcation of trajectories around $n=1$ , and thus the bounce-averaged theory for this space requires three solution domains.

As the magnetic geometry and potentials become more complex, the single-valuedness can quickly break down. For instance, consider a case where the magnetic field is constant, but the potential has an intermediate peak (figure 4). Now, as in the 1-D double well in figure 1, $\epsilon$ and $\mu$ are not sufficient to define a particle trajectory, since a particle with insufficient $\epsilon$ to traverse the peak can be trapped on one side or the other. For instance, in figure 4, in region $r_1$ where $(\epsilon ,\mu )$ is below line 1 and above line 2, the particle can either be on a trajectory that is trapped in segment $n=0$ or a trajectory trapped in segment $n=2$ . There is no reason for particles on these different trajectories to share a value for the distribution function; thus, $f(\epsilon ,\mu )$ fails to be single-valued in this region. (A similar bifurcation occurs for a local maximum of the magnetic field; the only difference is that changing $B$ would change the slope of the lines in $(\epsilon ,\mu )$ space, rather than their $\epsilon$ -intercept.)

When such a bifurcation of the trajectory happens, the domain on which the BAFP equations are solved must be split, to keep $f(\epsilon ,\mu )$ single-valued on each domain. Each domain consists of (a) an area in $(\epsilon ,\mu )$ space and (b) a continuous subset of the field line segments $n$ over which the bounce-average is taken. For the case in figure 4, three domains are required: the domain $d_1$ consisting of region $r_2$ , bounce-averaged over field line segments $n \in \{0,1,2\}$ ; the domain $d_2$ consisting of regions $r_0$ and $r_1$ , bounce-averaged over field line segments $n \in \{0\}$ ; and the domain $d_3$ consisting of region $r_1$ , bounce-averaged over field line segments $n \in \{2\}$ .

To solve the problem, the domains must be stitched together. This requires two boundary conditions at the boundary (line 1). First, $f(\epsilon ,\mu )$ must be continuous at the boundary. Second, the phase-space flux must be continuous, i.e.

(3.2) \begin{align} n_i \Gamma _Z^i \bigr |_{d_1} = n_i \Gamma _Z^i \bigr |_{d_2} + n_i \Gamma _Z^i \bigr |_{d_3}, \end{align}

where $n_i$ is a shared normal vector to the surface in $(\epsilon ,\mu )$ space.

Figure 5. Directed graph structure of boundary conditions for the field configuration in figure 4.

We can see from the flux conservation equation that the two domains $d_2$ and $d_3$ merge into $d_1$ . This suggests a representation of the relationship between the domains as a directed graph, with the directionality going from the domain representing a single passing trajectory to the domains representing bifurcated trajectories (figure 5).

3.1. Ambipolar fields and the need for an automatic method

Historically, especially in the study of tandem mirrors with thermal barriers, these domain relations have been imposed at the beginning of the problem and left fixed for the computation (Matsuda & Stewart Reference Matsuda and Stewart1986). However, this is likely to be insufficient for accurate simulations of modern mirrors, which often have potential profiles determined by sloshing kinetic ion distributions (Egedal et al. Reference Egedal, Endrizzi, Forest and Fowler2022; Endrizzi et al. Reference Endrizzi2023; Frank et al. Reference Frank2024). Correctly determining the resulting potentials requires solving for the electric potential in each region self-consistently over the course of the simulation so as to enforce quasineutrality, i.e. to enforce

(3.3) \begin{align} \sum _s Z_s \int _s f_s\, {\rm d}\boldsymbol{Z} = 0 \end{align}

for each region (Frank et al. Reference Frank2024; Kolmes et al. Reference Kolmes, Ochs and Fisch2025). As a result, it might not be known at the beginning of the simulation where the potential peak might occur, and thus over which field line segments and regions of $(\epsilon ,\mu )$ to define the domains. Thus, dynamic detection of such domains is a necessity.

Figure 6. (a) Discretised magnetic field $B$ and potential energy $\psi$ as a function of axial segment $n$ and (b) COM-space acessibility plot for a scenario exhibiting trajectory bifurcation without a local maximum in either potential or magnetic field. This ‘Yushmanov trajectory bifurcation’ process is discussed in more detail in Appendix A.

Making matters worse, the bifurcation of trajectories does not generally require such an easy-to-spot feature as a local peak in the electric or magnetic field (which is the basis for such calculations in stellarator neoclassical transport codes (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999; Velasco et al. Reference Velasco, Calvo, Parra and García-Regaña2020)). Consider the field arrangement in figure 6. In segment $n=1$ , this has neither a local maximum in the magnetic or electric field, and thus does not seem like it would represent a barrier between local wells. Nevertheless, as can be seen in figure 6, there is a region of $(\epsilon ,\mu )$ space that can be accessed by particles with $n \in \{0,2\}$ , but not for $n=1$ . Thus, this field set-up also has a trajectory bifurcation and must be solved using multiple grids.

The conditions required for the particular situation in figure 6, which we term a ‘Yushmanov trajectory bifurcation’, can actually be solved for (Appendix A); it turns out that such a bifurcation will occur when

(3.4) \begin{align} \frac {\partial \psi }{\partial B} \lt 0 \quad\text{and}\quad \frac {\partial ^2\psi }{\partial B^2} \gt 0. \end{align}

Nevertheless, in anticipating arbitrary field arrangements involving multiple wells and, particularly, in looking forward to non-relativistic plasmas in three COM dimensions, it is clearly useful to have a generalisable method that can take the allowed regions (in COM space) for each part (in $\boldsymbol{x}$ ) of the plasma, and return the requisite domains and boundary conditions. The development of such a method is the subject of the next sections.

4. Formalising the requirements for cohesive domains

As a starting point, we assume that we have a set of segments $\mathcal{N} = \{0,\ldots , N-1\}$ , with allowed areas of $(\epsilon ,\mu )$ space defined for each segment $n \in \mathcal{N}$ , of the form $\epsilon \geqslant b_n(\mu )$ . This will in general chop the $(\epsilon ,\mu )$ space up into a number of different regions $r_i$ , with edges defined by the boundary functions $b_n(\mu )$ .

We note that in formulating the problem in this way, the generalisation to 3-D $(\epsilon ,\mu ,\varPhi )$ space is straightforward, since each allowed region will now be given by $\epsilon \geqslant b_n(\mu ,\varPhi )$ , with the boundaries between regions given by 2-D surfaces rather than 1-D curves. The following discussion thus applies immediately to more complicated and higher-dimensional problems, though to keep things clear, we will restrict the discussion to non-relativistic $(\epsilon ,\mu )$ space.

Now, consider a specific region $r_i$ . This region will have a set $\mathcal{N}_i \subseteq \mathcal{N}$ of the segments for which $r_i$ is an allowed region. Ordering the set $\mathcal{N}_i$ in terms of increasing $n$ , we can in general then split $\mathcal{N}_i$ into maximal continuous subsequences $\mathcal{C}_i^m$ , defined such that each $\mathcal{C}_i^m$ has no breaks in the integer sequence of $n$ values. The set of all $\mathcal{C}_i^m$ for a given region $r_i$ we denote $\boldsymbol{\mathcal{C}}_i$ . For a given $r_i$ , the size (cardinality) of the set $\boldsymbol{\mathcal{C}}_i$ represents the number of distinct trajectories in that region. For instance, returning to figure 6, the indicated region in $(\epsilon ,\mu )$ , which we can label as region $r_0$ , has $\mathcal{N}_0 = \{0,2\}$ , $\mathcal{C}_0^0 = 0$ , $\mathcal{C}_0^1 = 2$ and $\boldsymbol{\mathcal{C}}_0 = \{\mathcal{C}_0^0,\mathcal{C}_0^1\} = \{\{0\},\{2\}\}$ . The fact that $|\boldsymbol{\mathcal{C}}_0| = 2$ means that there are two distinct trajectories; namely, those trapped in axial segment $n=0$ and those trapped in axial segment $n=2$ .

At this point, it is helpful to define a population, denoted $p_i^m$ , as the combination of the region $r_i$ (a volume in COM space) and a maximal continuous subsequence of allowed axial segments $C_i^m$ in that region. The population defines everything other than the boundary conditions about the bounce-averaged Fokker–Planck problem in COM space: namely, it defines the region over which the problem is to be solved and the axial segments over which the bounce average is to be taken. Boundary conditions then come from stitching together connecting populations.

However, solving the Fokker–Planck problem for each population so defined would be extremely inefficient, since there can be many regions even for problems (such as those in figures 2 and 3) that can be solved on a single domain. Thus, our goal will be to build up collections of compatible populations into maximal simulation domains.

Thus, define a domain $d$ as a collection of connected and pairwise-compatible populations $\{p_i^m\}$ . To define ‘connected’ and ‘pairwise-compatible’, it is helpful to consider what happens as we try to extend the domain to a new population. Thus, consider a population $p_i^m$ in domain $d$ , with region $r_i$ and continuous sequence $\mathcal{C}_i^m$ . Then, try to extend the domain $d$ by adding a population with region $r_j$ bordering $r_i$ . What are the conditions that allow this to happen?

First, it only makes sense to extend the domain if there is a population on $r_j$ that contacts trajectories from $p_i^m$ ; i.e. if

(4.1) \begin{align} \exists \mathcal{C}_j^n \in \boldsymbol{\mathcal{C}}_j : (r_i \text{ borders } r_j) \quad \mathrm{and}\quad (\mathcal{C}_j^n \cap \mathcal{C}_i^m \neq \{\}). \end{align}

If this condition is not satisfied, then trajectories in region $r_j$ are located in a spatially distinct part of the plasma from those trajectories described by $p_i^m$ and thus do not interact in a Fokker–Planck manner with the population $p_i^m$ . Thus, populations $p_i^m$ and $p_j^n$ are said to be connected if

(4.2) \begin{align} (r_i \text{ borders } r_j) \quad \mathrm{and} \quad (\mathcal{C}_j^n \cap \mathcal{C}_i^m \neq \{\}). \end{align}

However, connectedness is not sufficient to ensure that the domain can be extended. Consider, for instance, the regions $r_1$ and $r_2$ on either side of line 1 in figure 4. Then, consider the population $p_2^0 \equiv (r_2,\{0,1,2\})$ . Going to region $r_1$ , we see that there are two populations that are connected to $p_2^0$ : $p_1^0 \equiv (r_1,\{0\})$ and $p_1^1 \equiv (r_1,\{2\})$ . Precisely because there are multiple connected populations in region $r_1$ , the domain cannot be extended, because this multiplicity means that trajectories bifurcate at this point and a boundary condition (3.2) must be enforced at line 1. Thus, for $p_i^m$ in region $r_i$ to be compatible with a population $p_j^n$ in region $r_j$ , we must have

(4.3) \begin{align} \nexists \; l : (l\neq n)\quad \mathrm{and}\quad (\mathcal{C}_j^l \cap \mathcal{C}_i^m \neq \{\}). \end{align}

However, this is still not sufficient, since this compatibility relation must be reciprocal; condition (4.3) shows that $p_2^0$ is incompatible with the $p_1^0$ and $p_1^1$ in figure 5, but not that $p_1^0$ and $p_1^1$ are incompatible with $p_2^0$ . Thus, for $p_i^m$ in region $r_i$ to be compatible with a population $p_j^n$ in region $r_j$ , we must also have

(4.4) \begin{align} \nexists \; l : (l\neq m)\quad \mathrm{and} \quad (\mathcal{C}_i^l \cap \mathcal{C}_j^n \neq \{\}). \end{align}

If, for populations $p_i^m$ and $p_j^n$ , conditions (4.3)–(4.4) hold, then the populations are compatible and can exist in the same domain.

Unfortunately, it can be shown (Appendix B) that the compatibility property is not transitive. Thus, to add a population to the domain, connectedness must only be checked for one region, but compatibility must be checked with all populations already in the domain. This non-transitivity is why the domain was defined to be a connected and pairwise-compatible set of populations.

To completely partition the space for solution of the Fokker–Planck problem, we must assign every population to a domain. These domains are then connected by boundary conditions along shared boundaries. As shown in Appendix B, this decomposition is not, in general, unique. Nevertheless, an example algorithm for performing a valid decomposition is described in Appendix C.

At the end, each domain $d_a$ (which is a set of populations) will be associated with a set of connected regions $\mathbb{R}_a$ and a single continuous set of axial segments $\mathbb{C}_a$ . Thus, for each domain $d_a$ , the BAFP problem is to be performed over the area in $\boldsymbol{Z}$ that is spanned by $\mathbb{R}_a$ , with the bounce average taken over the axial segments contained in $\mathbb{C}_a$ .

5. Enforcing boundary conditions

Now that the space is partitioned into domains, we must define boundary conditions for each domain. As in §3, we need to enforce continuity and particle conservation at each boundary; the key aspect will be to identify which populations must be matched.

There are two types of boundaries. External boundaries are defined by the edge of accessible $\boldsymbol{Z}$ for the total modelled space. These boundaries are either reflecting, if they represent the boundary of 0 kinetic energy that defines accessibility, or approximately absorbing (Dirichlet with a value of 0), if they represent the loss cone. (Though it may be noted that in more collisional mirrors (Mirnov & Riutov Reference Mirnov and Riutov1979; Rognlien & Cutler Reference Rognlien and Cutler1980) or mirrors with significant secondary electron emission from the walls (Konkashbaev, Landman & Ulinich Reference Konkashbaev, Landman and Ulinich1978; Skovorodin Reference Skovorodin2019), a non-zero loss cone boundary may be appropriate.)

The remaining boundaries are internal, representing matching conditions between domains. In general, if we look at a domain $d_a$ , each internal boundary can be defined by two regions; the region $r_i$ (with population $p_i^m$ in domain $d_a$ ) and the region $r_j$ . If $r_j$ does not represent a region within the domain $d_a$ , then this represents an outer boundary of the domain $d_a$ , where the boundary conditions must be applied. We can denote this boundary by $b_{ij}^a$ .

Now, note that of the two regions $r_i$ and $r_j$ , one of these regions will always have a set of accessible axial segments that is a strict superset of the other region; i.e.:

(5.1) \begin{align} \mathcal{N}_{i} \subset \mathcal{N}_j \quad\text{ or }\quad \mathcal{N}_j \subset \mathcal{N}_i. \end{align}

Thus, we can define region $r_i$ to be `higher access’ than $r_j$ if $\mathcal{N}_i \supset \mathcal{N}_j$ , and ‘lower access’ if the reverse holds. In general, trajectories can bifurcate in going from the higher-access region to the lower-access region, but not vice versa.

Now let us establish the boundary conditions for boundary $b_{ij}^a$ between region $r_i$ (in domain $d_a$ ) and $r_j$ (not in domain $d_a$ ). We start by identifying the higher-access region. Then, taking the population for region $r_i$ in domain $d_a$ to be $p_i^m$ , we start by finding the connected populations in region $r_j$ . This set of populations will define a set $\mathcal{D}_{ij}^a$ of domains connected to domain $d_a$ at that boundary, i.e. the domains which contain those connected populations. If $r_i$ is the higher-access region, then the boundary conditions now take the form:

(5.2) $$\begin{align}f(\boldsymbol{Z})\bigr |_{b_{ij}^a,d_a} &\underbrace {= f(\boldsymbol{Z})\bigr |_{b_{ji}^b,d_b} = \ldots }_{d_b \in \mathcal{D}_{ij}^a}\end{align}$$

(5.3) $$\begin{align}n_i \Gamma _Z^i \bigr |_{b_{ij}^a,d_a} &= \sum _{d_b \in \mathcal{D}_{ij}^a} n_i \Gamma _Z^i \bigr |_{b_{ji}^b,d_b}.\end{align}$$

Note that if $r_j$ has no populations that are connected to $p_i^m$ , then $\mathcal{D}_{ij}^a = \{\}$ and the boundary conditions reduce to the reflecting boundary conditions of an external boundary.

If, instead, $r_i$ is the lower-access region, then there will be only one population $p_j^n$ that is connected to $p_i^m$ in domain $d_a$ . Denote the domain associated with $p_j^n$ as $d_b$ . Then, we can establish the boundary conditions as in the previous paragraph, but starting with $p_j^n$ in $d_b$ . Explicitly, we first find all the populations $p_i^l$ in region $r_i$ connected to $p_j^n$ and thus construct the set of domains $\mathcal{D}_{ji}^b$ containing those connected populations. Note that $\mathcal{D}_{ji}^b$ will, by definition, include $d_a$ . Then, the boundary conditions are

(5.4) $$\begin{align}f(\boldsymbol{Z})\bigr |_{b_{ji}^b,d_b} &\underbrace {= f(\boldsymbol{Z})\bigr |_{b_{ij}^c,d_c} = \ldots }_{d_c \in \mathcal{D}_{ji}^b}\end{align}$$

(5.5) $$\begin{align}n_i \Gamma _Z^i \bigr |_{b_{ji}^b,d_b} &= \sum _{d_c \in \mathcal{D}_{ji}^b} n_i \Gamma _Z^i \bigr |_{b_{ji}^b,d_b}.\end{align}$$

The above mentioned procedure defines all internal boundary conditions for the BAFP problem.

6. Graph structure of domain connections

The rules outlined in §4 and §5 are sufficient to properly set up a BAFP problem, to solve it either analytically or (more likely if application is necessary) computationally. However, it is also useful when thinking about a problem to be able to quickly visualise how the various domains fit together. The above mentioned rules, which are focused on conditions at each boundary (of which each domain can have quite a number), do not make developing this intuition particularly straightforward. The goal, then, is to quickly encode the boundary condition connections into the connections between domains, which can then be visualised.

The condition for two domains to be connected is simple: domains $d_a$ and $d_b$ will be connected if

(6.1) \begin{align} \exists p_i^m \in d_a, p_j^n \in d_b: (p_i^m \text{ and } p_j^n \text{ are connected}). \end{align}

Furthermore, we can establish a ‘direction of connection’, based on whether $p_i^m$ or $p_j^n$ is in a higher-access region. We draw the direction from the higher-access region to the lower-access region, representing the direction in which trajectories can bifurcate. In other words, if $\mathcal{N}_i \supset \mathcal{N}_j$ , the direction of connection is from $d_a$ to $d_b$ . Note that if two domains connect at multiple boundaries, the direction of connection can, in principle, be bi-directional.

With such connections established, the boundary condition structure of a BAFP problem can then be represented as a graph. This graph can help to clarify the multi-well structure of the BAFP problem, in particular, encoding the bifurcation structure of trajectories.

For instance, figure 7 shows a ‘tiered well’, where as $\epsilon$ decreases, trajectories first bifurcate at the local potential and magnetic maximum at $n=3$ , and then bifurcate again at the local potential and magnetic maxima at $n=1$ and $n=5$ . Although there are only three resulting regions in $(\epsilon ,\mu )$ space, the separate wells mean that there are seven separate domains (figure 8). The boundaries between these domains are connected as shown in figure 9. Domain $d_6$ represents the domain of trajectories which transit the entire device. Then, this bifurcates into domains $d_4$ and $d_5$ , which represent the trajectories trapped on either side of $n=3$ , but which pass over the potential maxima at $n=1$ and $n=5$ . Finally, the trajectories bifurcate again at these final potential maxima. In this way, the graph structure succinctly encodes the domain connections that result from the well structure.

Of course, the ‘tiered well’ is a very straightforward example with a very clear logic to the domain connections. The established rules for establishing cohesive domains and determining the connections between domains really prove their value when considering more complicated scenarios. For instance, consider the field arrangement in figure 10. One might scratch their head for quite a while coming up with a set of consistent domains to describe this problem, but the rules in §4, as implemented in the algorithm in Appendix C, can quickly come up with a set of consistent domains (figure 11). Furthermore, the methods described here can also quickly reveal the connections between these domains (figure 12). Thus, it can be seen that formulating the problem in this way opens up the possibility of solving arbitrarily complicated multi-well BAFP problems.

Figure 7. (a) A ‘tiered well’ field arrangement and (b) COM accessibility plot. In this scenario, in order of decreasing $\epsilon$ , trajectories first bifurcate around axial segment $n=3$ , then bifurcate again around segments $n=1$ and $n=5$ .

Figure 8. Domain decomposition for the ‘tiered well’ in figure 7. Each domain $d_a$ consists of a list of populations, and as a result has an associated set of regions $\mathbb{R}_a$ (the shaded area in each plot) and continuous set of axial segments $\mathbb{C}_a$ , given in the title of the plot. Here, several domains (e.g. $d_0$ , $d_1$ , $d_2$ and $d_3$ ) share the same region in $(\epsilon ,\mu )$ , but represent different populations because they occur at different positions along the field line.

Figure 9. Graph structure of domain connections for the ‘tiered well’ in figure 7, with the domains defined in figure 8. The forks in the graph are associated with the trajectory bifurcations, first at $n=3$ (splitting $d_6$ into $d_4$ and $d_5$ ), then again around segments $n=1$ and $n=5$ .

Figure 10. An arbitrary complicated field arrangement, with many crossings in the accessibility boundary lines.

Figure 11. An algorithmically solved consistent domain decomposition of the scenario in figure 10, showing the shaded set of regions $\mathbb{R}_a$ and continuous set of axial segments $\mathbb{C}_a$ for each domain $d_a$ .

Figure 12. Graph structure of domain connections for the scenario in figure 10, with the domains defined in figure 11. The graph structure is much more complicated than for the tiered well scenario in figures 7, 8 and 9, but if one chooses any set of connected nodes, it is possible to identify the boundary where the conditions are enforced.