\varv\varv

Representation and modeling of charged particle distributions in tokamaks

Andreas Bierwage¹ [email protected] Michael Fitzgerald² Philipp Lauber³ Mirko Salewski⁴ Yevgen Kazakov⁵ Žiga Štancar^6,2 ¹QST, Rokkasho and Naka Fusion Institutes, Japan
²Culham Centre for Fusion Energy of UKAEA, Culham Science Centre, Abingdon, United Kingdom
³Max-Planck-Institut für Plasmaphysik, Garching, Germany
⁴Department of Physics, Technical University of Denmark, Kgs. Lyngby, Denmark
⁵Laboratory for Plasma Physics, LPP-ERM/KMS, Partner in the Trilateral Euregio Cluster (TEC), Brussels, Belgium
⁶Jožef Stefan Institute, Ljubljana, Slovenia

Abstract

Experimental diagnostics, analysis tools and simulations represent particle distributions in various forms and coordinates. Algorithms to manage these data are needed on platforms like the ITER Integrated Modelling & Analysis Suite (IMAS), performing tasks such as archiving, modeling, conversion and visualization. A method that accomplishes some of the required tasks for distributions of charged particles with arbitrarily large magnetic drifts in axisymmetric tokamak geometry is described here. Given a magnetic configuration, we first construct a database of guiding center orbits, which serves as a basis for representing particle distributions. The orbit database contains the geometric information needed to perform conversions between arbitrary coordinates, modeling tasks, and resonance analyses. Using that database, an imported or newly modeled distribution is mapped to an exact equilibrium, where the dimensionality is reduced to three constants of motion (CoM). The orbit weight is uniquely given when the input is a true distribution: one that measures the number of physical particles per unit of phase space volume. Less ideal inputs, such as distributions estimated without drifts, or models of particle sources, can also be processed. As an application example, we reconstruct the drift-induced features of a distribution of fusion-born alpha particles in a large tokamak, given only a birth profile, which is not a function of the alpha’s CoM. Repeated back-and-forth transformations between CoM space and energy-pitch-cylinder coordinates are performed for verification and as a proof of principle for IMAS.

Tokamak plasma, equilibrium, particle distribution, guiding centers, fusion products

1 Introduction and workflow

1.1 Motivation, purpose and scope

The ITER Integrated Modelling & Analysis Suite (IMAS) PinchesFEC18 offers various ways to store particle distributions for the study of magnetically confined fusion plasmas. Different representations include different choices of coordinates and different discretization methods (mesh grids or marker particles). The Energetic Particle Topical Group of the International Tokamak Physics Activity (ITPA) is currently driving an action to equip IMAS with tools to model and convert distributions of fast particles between different representations that arise in experimental and computational work.

In experimental work, the observed position space of an energetic particle diagnostic is often given by the line-of-sight of the diagnostic in cylinder coordinates $(R,z,\zeta)$ , where $R$ is the major radius, $z$ the height and $\zeta$ the toroidal angle. In tokamaks, $\zeta$ is often an ignorable symmetry coordinate. For many diagnostics, the observed velocity-space is given by weight functions that are derived from energy and momentum conservation Salewski15 ; Salewski16 ; Salewski18a ; Salewski19 . The diagnostic weights are often expressed as functions of kinetic energy $E=Mv^{2}/2$ and velocity pitch $\lambda=v_{\parallel}/v$ , where $v$ is the magnitude of the velocity vector ${\bm{v}}$ , $v_{\parallel}$ its component parallel to the magnetic field vector ${\bm{B}}$ , and $M$ the particle mass. Diagnosticians hence have a good understanding of the observation regions of diagnostics in phase space $(E,\lambda,R,z)$ . Synthetic diagnostics calculating expected measurements based on numerical simulations are usually based on such representations of local velocity distributions.

Given the uncertainties of experimental measurements, no distinction is usually made between the gyro-averaged distribution of physical particles and the distribution of their guiding centers (GC), at least in the plasma core. The modeling of distribution functions is therefore often limited to GC statistics (except near the wall where the geometry of gyro-orbits matters for the interpretation of loss diagnostics and for the estimation of heat loads on plasma-facing hardware).

If one assumes that the plasma is not only toroidally symmetric but also quiescent, the distribution functions of confined charged particles and their GCs can be taken to depend only on three constants of motion (CoM). Analyses of resonant instabilities are usually performed with respect to such an equilibrated reference state, where the orbit time $\tau$ , gyrophase $\xi$ and toroidal angle $\zeta$ are ignorable symmetry coordinates. The respective conserved actions — namely, the kinetic energy $E$ , magnetic moment $\mu$ , and canonical toroidal angular momentum $P_{\zeta}$ — combined with an index $\sigma$ specifying the sign of $v_{\parallel}$ on passing orbits, constitute a natural set of CoM coordinates $(E,\mu,P_{\zeta};\sigma)$ that is common in theoretical analyses, although other useful (and equivalent) sets of CoM exist. Recently, an experimental measurement of a fast ion distribution function in three-dimensional (3-D) CoM space was demonstrated for fast-ion D- $\alpha$ spectroscopy Stagner22 . Efforts are underway to apply this orbit tomography method StagnerPhdThesis to other fast-ion diagnostics Jaerleblad21 .

In order to make predictions of a measurement based on a stability calculation or carry out stability analyses based on measured data, we must be able to transform distribution functions between 3-D CoM space and 4-D representations given in various coordinates as illustrated in Fig. 1. Related problems have been tackled to various degrees by different research groups with different codes and methods. For instance, TRANSP/NUBEAM BreslauTRANSP17 and LIGKA/HAGIS LauberFEC20 employ binning and smoothing techniques to map distributions represented by marker particles onto a mesh and compute their gradients in CoM space. While TRANSP/NUBEAM uses Monte-Carlo sampling, LIGKA/HAGIS uses a database of GC orbits for all particle species. Related activities were reported by the ASCOT TholerusASCOT17 and MEGA BierwageVSTART12 groups. These methods have been used for many years, but detailed documentations are rarely published. The above references BreslauTRANSP17 ; LauberFEC20 ; TholerusASCOT17 ; BierwageVSTART12 all point to presentations given at technical meetings. Useful elements can be found in many papers, but it is necessary to assemble that information into a complete goal-oriented workflow.

Refer to caption — Figure 1: Examples of coordinate representations and transformations (arrows) of a toroidally symmetric distribution function. We require that every 4-D distribution $f({\bm{Z}})$ be reduced to an exact equilibrium $f_{\rm orb}({\rm CoM})$ , so that the mapping between 4-D GC phase space and 3-D CoM orbit space becomes unique and readily reversible. — Position coordinates: major radius $R$ ; height $z$ ; poloidal magnetic flux $\Psi_{\rm P}(R,z)$ ; poloidal angle $\vartheta$ . — Velocity coordinates: GC velocity components $v_{\perp}$ and $v_{\parallel}$ relative to the magnetic field ${\bm{B}}$ ; kinetic energy $E=Mv^{2}/2$ with $v^{2}=v_{\parallel}^{2}+v_{\perp}^{2}$ ; pitch $\lambda=v_{\parallel}/v$ ; sign $\sigma=v_{\parallel}J_{\parallel}/|v_{\parallel}J_{\parallel}|$ relative to the parallel current density $J_{\parallel}$ . — Mixed coordinates: magnetic moment $\mu=Mv_{\perp}^{2}/(2B)$ ; conserved pitch $\Lambda=\mu B_{0}/E$ ; canonical toroidal angular momentum $P_{\zeta}$ ; coordinates $\Psi_{\rm P,ref}$ , $\lambda_{\rm ref}$ , etc., evaluated at a certain reference point, such as a midplane crossing or a turning point.

Thus motivated, this paper documents the workflow used in our computer code VisualStart BierwageVSTART12 ; Bierwage12a . Since 2009, we have used this code to initialize hybrid simulations with modeled or numerically computed beam ion distributions Bierwage13a ; Bierwage14a ; Bierwage17a ; Bierwage18 and to analyze their resonant interactions with Alfvén waves Bierwage14b ; Bierwage16a ; Bierwage16b . VisualStart represents GC distributions using marker particles, which can be directly used in full- $f$ simulations¹¹1For instance, see Fig. 3 in Ref. Bierwage13a , and Fig. 4 in Ref. Bierwage14a for fast ion tails, and Fig. 4 in Ref. Bierwage17a for a birth distribution. and yield a quiet start Bierwage12a . The representation in CoM space with an underlying GC orbit database will also allow us to directly embed the distribution functions in general phase-space transport theories Zonca15b ; Chen16 ; Falessi19 ; Zonca21 , both formally and technically.

Here we propose the orbit-based representation and modeling technique as a solution for ITER IMAS and other platforms that need to store and process distributions of charged particles in magnetically confined plasmas. The subject of smoothing is discussed briefly in the following section, but the documentation and verification of smoothing algorithms is left for future work, as is the associated problem of initializing delta- $f$ simulations that require the evaluation of gradients. We note that a similar and in many ways complementary effort has been undertaken by S. Benjamin et al. BenjaminBachelorThesis .

1.2 Basic rules for practical situations

The theoretical foundations for transforming distribution functions are relatively straightforward, but it is also clear that compromises and creativity are needed when tackling real-world data with numerical techniques, where accuracy has to be traded for computational performance and where data may have varying degrees of quality and completeness. For instance, the singularities and topological boundaries of CoM space Rome79 require attention when cutting a mesh to define phase space volume elements. Experimentally measured distributions and numerically computed distributions are often incomplete and noisy, so some amount of modeling is needed, including smoothing, interpolation and extrapolation.

In view of these challenges, some basic rules should be followed in order to maintain the physical and numerical integrity of the data, independently of the particular workflow used to process distribution functions:

(i)

Smoothing, interpolation and extrapolation are best applied at the preparation stage, where a new distribution function is first constructed or imported. Such potentially irreversible manipulations should be minimized in all subsequent operations.
(ii)
Pay attention to the fact that, in toroidal geometry, resolution and noise levels inevitably vary in space and in time, because the mirror force and magnetic drifts cause GCs to be distributed nonuniformly along their orbits. In particular:
- •
  
  Inherent topological discontinuities and singularities, as found in CoM space, should be preserved; e.g., through mesh accumulation and by smoothing only around, not across trapped-passing boundaries and loss cones.
- •
  
  The optimal representation method (maximizing accuracy) may differ between archiving a distribution function and using it in a simulation. While archiving merely requires adequate resolution, simulations should also minimize signal-noise correlations Bierwage12a .
These issues are most pronounced in but not limited to cases with large magnetic drifts.
(iii)

Ensure that the distribution to be transformed is an exact equilibrium in the given magnetic configuration; i.e., its shape must not vary in time when evolved using unperturbed equations of motion.

The equilibrium constraint (iii) is crucial, so it deserves further explanation: Any distribution function $f$ can be mapped to CoM space, where the result $f_{\rm orb}$ measures the time-averaged particle densities on unperturbed GC orbits. This operation is indicated by a blue arrow in Fig. 1. However, when the original $f$ is not an equilibrium, this mapping operation is time-dependent and it is impractical to store the information needed to invert it. The equilibrium constraint reduces the effective dimensionality to three CoM by determining the longitudinal distribution of GCs along their orbit, which is now time-independent and with it all coordinate transformations as well. The coordinate transformations $f\leftrightarrow f_{\rm orb}$ can then be readily performed in any direction, as indicated by a red double arrow in Fig. 1.

Nonequilibrium distributions are very likely to be encountered in IMAS, where one can expect to receive particle data from a measurement for which the magnetohydrodynamic (MHD) equilibrium is not exactly known or from a simulation where waves or magnetic islands were present. Enforcing the equilibrium constraint means that the input data are transformed into a distribution that is consistent with a modeled magnetic configuration, which generally differs somewhat from the conditions where the input data originated from.

1.3 Orbit-based representation and modeling

The above considerations suggest that it is useful and advantageous to represent distribution functions in CoM space with an underlying database of unperturbed GC orbits. With that geometric information, it is straightforward to reconstruct a 4-D distribution in arbitrary coordinates by loading a suitable number of markers uniformly in time along each GC orbit Bierwage12a . Projections to various coordinates and views of user-defined sub-spaces with user-defined resolution can then be readily obtained by binning the markers on demand. This minimizes the amount of transformations and associated data corruption, as demanded in items (i) and (ii) above. Marker particles are a natural choice for representing fast particle distributions since the latter are often computed using particle codes. Moreover, by organizing the markers on orbit contours, one can identify and deal with inherent singularities and discontinuities as required in item (ii). Last but not least, by representing a distribution function in CoM orbit space, the equilibrium constraint (iii) is automatically enforced.

We note that the implementation and enforcement of the equilibrium constraint only requires a conversion to CoM coordinates: $f\rightarrow f_{\rm orb}$ . The Jacobians required for such transformations can be determined using the Monte-Carlo method as implemented in TRANSP/NUBEAM by Breslau & Liu BreslauTRANSP17 . In this way, it is possible to avoid integration along GC orbits.

However, the computational cost associated with building an orbit database is more than repaid by the information one gains about the shape of the orbits, the longitudinal particle distribution, and the transit times. The transit times are needed for resonance analyses and their variation near trapped-passing boundaries can be used for mesh accumulation. The information about the orbit geometry is also useful for reduced models, for instance to capture finite-orbit-width effects Lauber18 , and it allows to define and make use of a variety of specialized GC coordinates, like the radii $R_{\rm m}$ where an orbit crosses the midplane. Last but not least, the information contained in the orbit database is also valuable for carrying out modeling tasks, such as the following example:

Given a model function $G_{\rm mdl}$ that may not be an exact solution of a kinetic equation for the GCs of physical particles, we wish to construct an equilibrium GC orbit distribution $f_{\rm orb}$ . A common example are models of the form $G_{\rm mdl}(E,\lambda,\Psi_{\rm P})$ that contain measured radial profiles as functions of flux $\Psi_{\rm P}$ and a velocity distribution computed by a bounce-averaged Fokker-Planck code. In such a case, the orbit database allows to incorporate the effect of magnetic drifts a posteriori Bierwage17c . There are, of course, various possibilities to perform this operation; in other words, there is some freedom to interpret the values of $G_{\rm mdl}$ as weights in an orbit distribution $f_{\rm orb}$ .²²2Functions such as $G_{\rm mdl}({\bm{Z}})$ described in this example will here be called “quasi-distributions”, because they do not measure the true number of objects of interest (here GCs) in a volume element ${\rm d}{\bm{Z}}$ but require more elaborate transformations than a mere conversion of coordinates. See Section 2.An intuitive and physically meaningful choice is to perform an orbit-time average $\left<...\right>_{\rm orb}$ followed by a suitable normalization: $f_{\rm orb}={\rm const.}\times\left<G_{\rm mdl}\right>_{\rm orb}$ .

1.4 Workflow

The workflow described in this paper is shown schematically in Fig. 2 and consists of five steps:

•

Step 1: Define an axisymmetric magnetic field ${\bm{B}}$ ; e.g., using a dedicated MHD equilibrium code.
•
Step 2: Define samples in 3-D CoM orbit space, with cell indices $(i,j,k)$ . The blue and red workflow branches in Fig. 2 begin with different inputs:
1. (blue)
  
  Import marker particles from another code as in earlier works Bierwage14a ; Bierwage17a ; Bierwage18 . The weights of the imported samples are reinterpreted as orbit weights: $w_{\rm import}\rightarrow W_{\rm orb}(i,j,k)$ .
2. (red)
  
  Use a meshing algorithm to divide the CoM space into cells of size $\Delta\mathcal{V}_{\rm orb}(i,j,k)$ . This path is described in detail in this paper.
•

Step 3: Integrate along GC orbits starting from the samples defined in Step 2. Store in a database.
•

Step 4: Represent an equilibrium distribution function $f_{\rm orb}$ using weighted markers loaded on the orbits in the database. When loading $N_{\tau}$ markers uniformly in time, their weight is an $N_{\tau}$ ’s fraction of the orbit weight: $w=W_{\rm orb}/N_{\tau}$ . In the case of the red workflow in Fig. 2, $W_{\rm orb}=\Delta\mathcal{V}_{\rm orb}\times f_{\rm orb}$ is the product of the volume element $\Delta\mathcal{V}_{\rm orb}$ and the desired phase space density function $f_{\rm orb}$ . The latter may be imported or constructed from a model.
•

Step 5: Visualize and verify. Iterate if necessary; e.g., changing the resolution in CoM space (Step 2), the number of markers on the orbits, or the weight assignment (Step 4). Export when ready.

Note that MHD wave spectra can be obtained from the data available in Step 1, and resonance maps can be computed from the orbit database built in Step 3. Procedures for smoothing the distribution, computing phase space gradients and analyzing the stability of resonant modes remain to be added in or around Step 5.

1.5 Outline

Clear definitions of symbols and terminology are a prerequisite for successful conversion operations. Thus, we begin our treatise in Section 2 with a discussion of different classes of distributions functions and how they are transformed. Step 1 of our workflow is part of Section 3, where we describe the tokamak plasma that serves as a working example, define coordinates and normalizations. Step 2, the construction of a phase space mesh, is covered in Section 4, where we also discuss orbit classes and the phase space topology in our scenario. The GC equations of motion we solve to construct the orbit database in Step 3 are given in Section 5, and the formula for the orbit volume elements is derived in Section 6. Steps 4 & 5 are covered in Sections 7 and 8, where we construct a distribution function $f$ from a model, visualize and analyze the result, and verify the integrity and accuracy of our scheme by transforming $f$ back and forth between two different representations. Modeling is a rich subject on its own, and our example is only meant to highlight the utility of an orbit database and to serve as a proof-of-principle for IMAS applications. The Appendices contain supplementary information on binning procedures and numerical accuracy.

2 Classes of distributions and transformations

There are not only different coordinates and different numerical representations of distributions; the latter also come in different types whose identity must be clearly defined since it determines how they are transformed. Here, we distinguish between three types of distributions: (i) density functions, (ii) histograms, and (iii) quasi-distributions.

The first two can be characterized as follows:

(i)

The symbols $f({\bm{Z}})$ represent phase space densities, which have units of $[{\rm m}^{-6}{\rm s}^{3}]$ and whose values are independent of the coordinates ${\bm{Z}}$ appearing in the argument, so $f({\bm{Z}}({\bm{v}},{\bm{x}}))=f({\bm{v}},{\bm{x}})$ . The argument merely specifies the representation space. The same is true for velocity and volume integrals:


$\displaystyle n({\bm{x}})=\int{\rm d}^{3}{\bm{v}}\,f$	$\displaystyle:\text{density field},$	(1a)
$\displaystyle\nu({\bm{v}})=\int_{V}{\rm d}^{3}{\bm{x}}\,f$	$\displaystyle:{\parbox{82.51282pt}{velocity distribution \\ in a volume $V$.}}$	(1b)

(ii)

The symbols $h({\bm{Z}})$ represent histograms, which have the units of the inverse volume element $|{\rm d}{\bm{Z}}|^{-1}$ of the coordinates used for “binning” (counting the objects of interest). Histograms are convenient to construct and integrate since no Jacobian is required (except for coordinate conversions).

The phase space density function $f$ is related to histogram functions $h$ , $h^{\prime}$ , etc., as

f({\bm{Z}})=\mathcal{J}^{-1}_{\bm{Z}}h({\bm{Z}})=\mathcal{J}^{-1}_{{\bm{Z}}^{\prime}}h^{\prime}({\bm{Z}}^{\prime})=...

(2)

where $\mathcal{J}_{\bm{Z}}$ is the Jacobian for the transformation from Cartesian coordinates ${\bm{Z}}_{\rm c}=({\bm{v}},{\bm{x}})=(v_{x},v_{y},v_{z},x,y,z)$ to another arbitrary set of phase space coordinates ${\bm{Z}}$ .

A phase space density function $f$ is an exact solution of a kinetic equation for the objects of interest; in our case, guiding centers (GC). In this work, we assume that all GC coordinates are defined such that they preserve the form of the kinetic equation

$\displaystyle 0=$	$\displaystyle\;(\partial_{t}+\dot{\bm{Z}}_{\rm c}\cdot\partial_{{\bm{Z}}_{\rm c}})f$	(3)
$\displaystyle=$	$\displaystyle\;\partial_{t}f+(\underbrace{\dot{\mathcal{J}}_{\rm gc}}\limits_{0}{\bm{Z}}_{\rm gc}+\mathcal{J}_{\rm gc}\dot{\bm{Z}}_{\rm gc})\cdot\frac{\partial_{{\bm{Z}}_{\rm gc}}f}{\mathcal{J}_{\rm gc}}$
$\displaystyle=$	$\displaystyle\;(\partial_{t}+\dot{\bm{Z}}_{\rm gc}\cdot\partial_{{\bm{Z}}_{\rm gc}})f.$	(4)

This is the case when the Jacobian $\mathcal{J}_{\rm gc}$ for the transformation ${\bm{Z}}_{\rm c}\rightarrow{\bm{Z}}_{\rm gc}$ from Cartesian to GC coordinates satisfies the phase space conservation law

\dot{\mathcal{J}}_{\rm gc}\equiv\partial_{t}\mathcal{J}_{\rm gc}+\partial_{{\bm{Z}}_{\rm gc}}\cdot\left(\mathcal{J}_{\rm gc}\dot{\bm{Z}}_{\rm gc}\right)=0.

(5)

Equation (5) is trivially satisfied for canonical coordinates, since their volume elements differ only by a constant factor, so $\mathcal{J}_{\rm gc}={\rm const}$ . All our CoM coordinates — namely, $(E,\mu,P_{\zeta};\sigma)$ and any equivalent set of GC orbit coordinates, such as $(v,v_{\parallel,{\rm ref}},R_{\rm ref})$ — are, by definition, conserved by the GC equations of motion that we solve. This means that, even when the equations of motion are not written in Hamiltonian form, our orbit coordinates combined with the triplet of ignorable angle variables $(\tau,\xi,\zeta)$ — orbit time, gyrophase, toroidal angle — can be considered to be canonical coordinates, whose Jacobians are constants satisfying Eq. (5) down to the numerical accuracy of our particle-pushing algorithm.

A well-known set of noncanonical GC coordinates is $(\mu,v_{\parallel},{\bm{x}}_{\rm gc})$ , whose volume elements transform as

{\rm d}^{3}{\bm{v}}{\rm d}^{3}{\bm{x}}=\frac{2\pi B^{*}_{\parallel}}{M}{\rm d}\mu{\rm d}v_{\parallel}{\rm d}^{3}{\bm{x}}_{\rm gc}.

(6)

The Jacobian factor $B^{*}_{\parallel}$ satisfies Eq. (5) as shown by Littlejohn Littlejohn83 . Equation (6) is very useful, because it can be easily converted to other velocity coordinates that are defined locally in position space — such as the pitches $\lambda=v_{\parallel}/v$ and $\Lambda=\mu B_{0}/E$ — since their Jacobian factors can be derived analytically (A) Moseev19 .

Both (i) densities and (ii) histograms are true distributions in the sense that, when integrated over a phase space volume element $\Delta\mathcal{V}$ ,

\int_{\Delta\mathcal{V}}{\rm d}{\bm{Z}}\;h({\bm{Z}})=\int_{\Delta\mathcal{V}}{\rm d}^{3}{\bm{v}}{\rm d}^{3}{\bm{x}}\;f=\mathcal{N}(\Delta\mathcal{V}),

(7)

they yield the true number $\mathcal{N}$ of objects contained inside that volume element. Functions or data that do not satisfy this condition for the specified kind of objects constitute the third type in our classification scheme:

(iii)

The symbol $G({\bm{Z}})$ represents quasi-distributions, which do not yield the true number of objects by mere volume integration in the configuration at hand, but require more complex transformations.

Quasi-distributions often arise in experimental or modeling work. A typical example in experimental measurements is the line-of-sight density, where the transformation $G\rightarrow f$ takes the form of an inverse problem.

In fact, since the magnetic configuration in an experiment is rarely known accurately, any experimentally measured particle statistics should be treated as quasi-distributions when used as input for modeling distributions in a numerically constructed MHD equilibrium. Conversely, it should be kept in mind that the spatial distribution of GCs differs from that of physical particles since the latter is broadened by gyration around the GCs. This means that, strictly speaking, the GC distributions that we deal with in this work constitute quasi-distributions from the viewpoint of physical particles. Even if our GC phase space density function $f$ is normalized to give the correct number of physical particles in the entire plasma volume, $\mathcal{N}(\mathcal{V})=\mathcal{N}_{\rm phys}(\mathcal{V})$ , the values of $\mathcal{N}(\Delta\mathcal{V})$ and $\mathcal{N}_{\rm phys}(\Delta\mathcal{V})$ given by Eq. (7) for an arbitrary subvolume $\Delta\mathcal{V}$ may differ in the case of fast ions with a steep density gradient.

In general, the process of “modeling” essentially consists of choosing a physically meaningful transformation $G\rightarrow f$ . We will return to this topic in Sections 7 and 8, where we apply our methods to model a distribution of fusion-born alpha particles, based on a quasi-distribution $G(\Psi_{\rm P})$ that represents their birth profile.

Finally, we emphasize that when a GC orbit distribution function is given in the form $f_{\rm orb}(E,\mu,P_{\zeta};\sigma)$ , or an equivalent set of CoM coordinates, the Jacobian $\mathcal{J}_{\rm CoM}$ must be accompanied by instructions concerning how to deal with the index $\sigma$ that specifies the sign of the parallel velocity. First, it must be clarified whether $\sigma$ is the sign relative to the magnetic field or the plasma current if their directions differ. Second, it must be clarified how particles trapped in a magnetic mirror are to be counted, because $\sigma$ is redundant for such trapped particle orbits. Depending on the choice made, the integral in Eq. (7) can be written in at least three ways:


$\displaystyle\mathcal{N}=$	$\displaystyle\sum\limits_{\sigma=\pm 1}\int{\rm d}E{\rm d}\mu{\rm d}P_{\zeta}\underbrace{\mathcal{J}_{\rm CoM}^{\rm t/p}({\rm trap./pass.})f_{\rm orb}}\limits_{h_{\rm orb}^{\rm t/p}(E,\mu,P_{\zeta};\sigma)},$	(8a)
$\displaystyle=$	$\displaystyle\sum\limits_{\sigma_{\rm pass}=0,\pm 1}\int{\rm d}E{\rm d}\mu{\rm d}P_{\zeta}\underbrace{\mathcal{J}_{\rm CoM}f_{\rm orb}}\limits_{h_{\rm orb}(E,\mu,P_{\zeta};\sigma_{\rm pass})},$	(8b)
$\displaystyle=$	$\displaystyle\sum\limits_{\sigma_{\rm HFS}=\pm 1}\sum\limits_{\sigma_{\rm LFS}=\pm 1}\int{\rm d}E{\rm d}\mu{\rm d}P_{\zeta}\underbrace{\frac{\mathcal{J}_{\rm CoM}}{2}f_{\rm orb}}\limits_{h_{\rm orb}(E,\mu,P_{\zeta};\sigma)}.$	(8c)

In the first case (8a), the summation over $\sigma$ counts each trapped orbit twice, which means that the histogram $h^{\rm t/p}_{\rm orb}$ has two sets of identical entries each representing half of the particles in the domain of trapped orbits. In the second case (8b), the summation index $\sigma_{\rm pass}=\{0,\pm 1\}$ is three-valued, identifying trapped ( $0$ ), co- and counter-passing orbits ( $\pm 1$ ), so that $h_{\rm orb}$ represents the full number of particles in all regions of GC phase space without duplicate entries. In mathematical form:

	$\displaystyle h_{\rm orb}(\sigma_{\rm pass}=0)=2\times h_{\rm orb}^{\rm t/p}({\rm trap.}),$		(9)
	$\displaystyle\mathcal{J}_{\rm CoM}=\mathcal{J}_{\rm CoM}^{\rm t/p}({\rm pass.})=2\times\mathcal{J}_{\rm CoM}^{\rm t/p}({\rm trap.}).$		(10)

The values of the density function $f_{\rm orb}$ itself are, of course, independent of the convention used. The point we want to make is that it is crucial to ensure consistency between the summation over $\sigma$ and the Jacobian when counting marker particles by integrating $f_{\rm orb}$ in $(E,\mu,P_{\zeta};\sigma)$ space as in Eq. (8).

In this paper, we adopt the third option (8c), with the single-valued constant Jacobian $\mathcal{J}_{\rm CoM}/2$ . Our CoM phase space mesh defined in Section 4 will double-count all orbits — passing and trapped — when they cross the plasma midplane: once on the high-field side (HFS) and once on the low-field side (LFS). The factor $2$ dividing the Jacobian in Eq. (8c) compensates this.

3 Geometry, coordinates and normalization

System geometry

The magnetic field vector is written

	$\displaystyle{\bm{B}}={\bm{\nabla}}\times{\bm{A}}=$	$\displaystyle{\bm{\nabla}}\zeta\times{\bm{\nabla}}\Psi_{\rm P}+B_{\zeta}{\bm{\nabla}}\zeta$		(11)
	$\displaystyle=$	$\displaystyle{\bm{\nabla}}\zeta\times{\bm{\nabla}}\Psi_{\rm P}+{\bm{\nabla}}\Psi\times{\bm{\nabla}}\vartheta,$		(12)

where $2\pi\Psi=\frac{1}{2\pi}\int{\rm d}\Psi{\rm d}\vartheta{\rm d}\zeta\;\mathcal{J}_{\rm f}^{\bm{x}}{\bm{B}}\cdot{\bm{\nabla}}\zeta$ is the toroidal flux and $2\pi\Psi_{\rm P}=-2\pi A_{\zeta}=\frac{1}{2\pi}\int{\rm d}\Psi{\rm d}\vartheta{\rm d}\zeta\;\mathcal{J}_{\rm f}^{\bm{x}}{\bm{B}}\cdot{\bm{\nabla}}\vartheta$ the poloidal flux. These fluxes are related via the tokamak safety factor $q={\rm d}\Psi/{\rm d}\Psi_{\rm P}$ measuring the mean magnetic field line pitch. $\mathcal{J}_{\rm f}^{\bm{x}}=[{\bm{\nabla}}\Psi\cdot({\bm{\nabla}}\vartheta\times{\bm{\nabla}}\zeta)]^{-1}=1/B^{\zeta}$ is the Jacobian for the transformation from Cartesian position coordinates ${\bm{x}}$ to the toroidal flux coordinates $(\Psi,\vartheta,\zeta)$ , where $\vartheta$ is a poloidal angle and $\zeta=-\varphi$ is the geometric toroidal angle, whose orientations are indicated in Fig. 3. Also shown in Fig. 3 are the axes of the right-handed cylinder coordinate system $(R,z,\zeta)$ , with major radius $R=\sqrt{x^{2}+y^{2}}$ and vertical coordinate $z$ .

The magnetic flux surface geometry and profiles for the MHD equilibrium that we use as a working example in this paper is shown in Fig. 4. This plasma is partly based on recent experiments with central radio-frequency (RF) heating using a 3-ion scheme Nocente20 ; Kazakov21 performed at the Joint European Torus (JET). The magnetic X-point has been removed and the profile of the safety factor $q$ is chosen to increase monotonically from values near unity in the core ( $q_{0}=0.98$ on axis) to $q_{a}=5.44$ at the boundary. This simulation scenario is also used in ongoing studies of MHD instabilities and fast ion physics Bierwage22b . The $q$ profile in the above-mentioned experiments varied dynamically and is likely to differ from ours at most times.

Besides $R$ and $z$ , we will also use the auxiliary coordinates $X=R-R_{0}$ and $Y=z-z_{0}$ measuring horizontal and vertical positions relative to the magnetic axis, which is located here at $(R_{0},z_{0})=(2.98\,{\rm m},0.27\,{\rm m})$ .

Magnetic flux labels

The normalized toroidal and poloidal fluxes ( $\psi$ , $\psi_{\rm P}$ ) or their square roots ( $\rho_{\rm T}$ , $\rho_{\rm P}$ ) can serve as convenient minor radial coordinates:

\rho_{\rm T}^{2}=\psi=\frac{\int_{0}^{V(\Psi)}{\rm d}V/\mathcal{J}_{\rm f}^{\bm{x}}}{\int_{0}^{V(\Psi_{\rm a})}{\rm d}V/\mathcal{J}_{\rm f}^{\bm{x}}},\quad\rho_{\rm P}^{2}=\psi_{\rm P}=\frac{\Psi_{\rm P}-\Psi_{\rm P,0}}{\Psi_{{\rm P},a}-\Psi_{\rm P,0}},

(13)

where $V(\Psi)$ is the volume within the flux surface labeled $\Psi$ . $\Psi_{\rm P,0}$ and $\Psi_{{\rm P},a}$ are the values of the poloidal flux at the magnetic axis ( $r=0$ ) and boundary ( $r=a$ ). For the toroidal flux we have $\Psi_{a}-\Psi_{0}=\int_{\Psi_{{\rm P},0}}^{\Psi_{{\rm P},a}}{\rm d}\Psi_{\rm P}\,q$ . In some occasions, we use a volume-averaged minor radius $0\leq r(\Psi)\leq a$ , which is defined here as

r^{2}=\frac{R_{0}}{\pi}\int_{0}^{V(r)}\frac{{\rm d}V}{\mathcal{J}_{\rm f}^{\bm{x}}B_{\zeta}}=\frac{R_{0}}{\pi}\int_{0}^{\Psi(r)}\frac{{\rm d}\Psi}{B_{\zeta}}\approx a^{2}\psi.

(14)

Our $r$ reduces the geometric minor radius in the limit of a cylindrical plasma with circular cross-section. For $B_{\zeta}\approx R_{0}B_{0}$ , we have $r/a\approx\rho_{\rm T}$ , and $a\approx[2(\Psi_{a}-\Psi_{0})/B_{0}]^{1/2}$ is the minor radius of the plasma boundary.

Particle and orbit coordinates

In the guiding center (GC) description, the velocity space consists of the ignorable gyrophase $\xi$ , the perpendicular particle velocity $v_{\perp}$ , and the parallel GC velocity $v_{\parallel}={\bm{v}}_{\rm gc}\cdot\hat{\bm{b}}$ relative to the magnetic field vector ${\bm{B}}$ with unit vector $\hat{\bm{b}}\equiv{\bm{B}}/B$ and magnitude $B=|{\bm{B}}|$ . The canonical momentum can be written in the form ${\bm{P}}={\bm{A}}/B_{0}+\Omega_{0}^{-1}{\bm{v}}$ for a particle with gyrofrequency $\Omega_{0}=QeB_{0}/M$ , charge $Qe$ and mass $M$ . In GC coordinates, its covariant toroidal component,

P_{\zeta}={\bm{P}}_{\rm gc}\cdot\partial_{\zeta}{\bm{x}}=-\frac{\Psi_{\rm P}}{B_{0}}+\frac{v_{\parallel}}{\Omega_{0}}\frac{B_{\zeta}}{B},

(15)

is called the canonical toroidal angular momentum and is chosen here to have units of area per radian, $(2\pi)^{-1}[{\rm m}^{2}]$ . Together with the lowest-order magnetic moment $\mu$ and the signed kinetic energy $\mathcal{E}=\sigma_{\rm ref}E$ ,

\mu=\frac{Mv_{\perp}^{2}}{2B},\quad\mathcal{E}=\sigma_{\rm ref}\frac{M}{2}v^{2},\quad\sigma=\frac{v_{\parallel}J_{\parallel}}{|v_{\parallel}J_{\parallel}|},

(16)

this yields a set of coordinates $(\mathcal{E},\mu,P_{\zeta})$ that belongs to the class of constants of motion (CoM), which are conserved by our equations of motion (presented in Section 5) when they are solved in a time-independent field ${\bm{B}}$ that is symmetric along $\zeta$ (or “axisymmetric” around $\hat{\bm{e}}_{z}$ ), satisfying $\partial_{\zeta}B_{\zeta}=\partial_{\zeta}\Psi_{\rm P}=0$ .

The sign $\sigma$ of the parallel velocity is measured relative to the parallel current density $J_{\parallel}={\bm{J}}\cdot{\bm{B}}/B$ . While redundant or unnecessary for orbits trapped in a magnetic mirror (both signs are present), $\sigma$ identifies whether a passing orbit with coordinates $(E,\mu,P_{\zeta})$ is co- or counter-going. In mathematical treatments, $\sigma$ usually appears as a separate index, $(E,\mu,P_{\zeta};\sigma)$ , but this is inconvenient for numerical representations. In our orbit database, we define $\sigma_{\rm ref}$ at a certain reference point; namely, the starting point used for the orbit calculation. This or any equivalent convention allows us to attach $\sigma_{\rm ref}$ to $E$ as in Eq. (16) and obtain a compact representation where all CoM arrays in the computer code have the same dimensionality: three.

It is also useful to define several different coordinates measuring the velocity pitch:

\Lambda\equiv\frac{\mu B_{0}}{E}=\frac{\cos^{2}\alpha}{\hat{B}},\quad\sin\alpha=\lambda\equiv\frac{v_{\parallel}}{v}=\sigma_{B}\sqrt{1-\Lambda\hat{B}},

(17)

with $\hat{B}=B/B_{0}\approx R_{0}/R$ and $\sigma_{B}\equiv v_{\parallel}/|v_{\parallel}|$ . Each pitch coordinate has its advantages in practical applications:

•

The coordinate $\lambda=v_{\parallel}/v$ is widely used together with $E$ or velocity $v=\sqrt{2E/M}$ . The local velocity space Jacobian of this set is a CoM ( $J_{E\lambda}\propto\sqrt{E}$ , Eq. (62)), and the Fokker-Planck equation can be integrated analytically in these coordinates Weiland18 .
•

The coordinate $\alpha$ is a true angle, so its value measures the velocity pitch in a very intuitive form. Moreover, the domains of co- and counter-passing particles are enlarged along the $\alpha$ -axis, making this coordinate a good choice for sampling the velocity distribution of tangentially injected beam ions Bierwage18 and distributions produced by central RF heating using a 3-ion scheme Nocente20 ; Kazakov21 ; Stancar21 .
•

The coordinate $\Lambda=\hat{B}^{-1}\cos^{2}\alpha$ has the advantage of being a CoM and that (in contrast to $\mu$ ) its upper limit is independent of energy: $0\leq\Lambda\leq\hat{B}_{\rm min}^{-1}$ , where $B_{\rm min}={\rm min}\{\hat{B}\}$ is the minimal field strength in the considered domain. This makes $\Lambda$ a convenient choice for defining the lines along which to divide the phase space as will be done in Section 4.

Unlike the conserved quantities $E$ , $\mu$ and $\Lambda$ , the pitch coordinates $\alpha$ and $\lambda$ vary along a GC orbit. The latter can be turned into orbit coordinates (i.e., CoM), $\alpha_{\rm ref}$ and $\lambda_{\rm ref}$ , by evaluating them at a well-defined reference point, as we did earlier for the sign $\sigma_{\rm ref}$ .

Reference points for orbit coordinates

In cases where we import computed particle distributions (blue workflow in Fig. 2), the initial reference point is taken to be the position of an imported marker particle, which may be located anywhere inside the plasma or in the surrounding vacuum.

In cases where we need to construct a new CoM mesh (red workflow in Fig. 2) the reference points lie in the midplane, which (at least in a usual tokamak plasma) contains all O-type stagnation points and is defined as the curve $z_{\rm m}(X)$ along which the magnetic field is perpendicular to its gradient,

z_{\rm m}(R)\;:\;{\bm{B}}\cdot{\bm{\nabla}}B=0,

(18)

The midplane is shown as a dash-dotted green line in Fig. 4. It contains the magnetic axis, $(R,z)=(R_{0},z_{0})$ , but deviates from the horizontal plane $z=z_{0}$ in up-down asymmetric plasmas like ours. Every GC orbit crosses the midplane twice, once in a region with higher field strength and once in a region with lower field strength. Figure 5 shows an example, where the orbit’s high- and low-field side crossings of the midplane are labeled “HFS” and “LFS”, respectively.

Normalization

Spatial positions and lengths are given in meters unless stated otherwise. Velocities are normalized by a reference value $v_{0}$ . Energy is normalized by $Mv_{0}^{2}$ (for each particle species) and the magnetic field by its on-axis value $B_{0}$ (here $3.7\,{\rm Tesla}$ ):

\hat{v}=\frac{v}{v_{0}},\quad\hat{E}=\frac{E}{Mv_{0}^{2}}=\frac{\hat{v}^{2}}{2},\quad\hat{\mu}=\frac{\mu B_{0}}{Mv_{0}^{2}}=\frac{\hat{v}_{\perp}^{2}}{2\hat{B}}.

(19)

4 CoM mesh and drift orbit types

In this section, we describe how we sample the GC orbit space. The example in Fig. 6 is used for illustration. A low resolution is chosen in order to make all grid points clearly visible. All coordinates shown in Fig. 6 are evaluated at the height $z_{\rm m}(X)$ of the midplane as defined in Eq. (18). This means that all coordinates appearing in this section are constants of motion (CoM). The set $(\hat{E},\alpha,X)$ used in panels (a) and (b) yields particularly compact and, in our opinion, intuitive plots. Another view of pitch-position space in $(\Lambda,X)$ coordinates is shown in panel (c). The plasma boundary is also taken to be the loss boundary for the GC orbits, so we do not include the vacuum here and our midplane mesh covers only the width of the plasma itself: $-1\,{\rm m}\lesssim X\lesssim 0.84\,{\rm m}$ .

The reference velocity used for normalization is chosen to be $v_{0}=1.3\times 10^{7}\,{\rm m/s}$ , which corresponds to $E_{0}=Mv_{0}^{2}/2\approx 3.5\,{\rm MeV}$ for fusion-born alpha particles ( ${}^{4}_{2}{\rm He}^{+2}$ ) or $1.75\,{\rm MeV}$ for deuterons ( ${}^{2}_{1}{\rm H}^{+}$ ). Both species have the same charge-to-mass ratio $Qe/M$ and their characteristic gyroradius in this example is $\rho_{0}=v_{0}/\Omega_{0}\approx 0.07\,{\rm m}\approx 0.02\times R_{0}$ . We consider the full energy range $0\leq\hat{E}\leq\hat{E}_{0}$ with $\hat{E}_{0}=0.5$ and define a mesh consisting of $N^{\rm(g)}_{E}$ grid points that are the boundaries of $N_{E}=N^{\rm(g)}_{E}-1$ cells. The value of energy at the center of cell $i$ is denoted $\hat{E}_{i}$ . The example in Fig. 6(a) has $N_{E}=15$ cells with equal sizes $\Delta\hat{E}_{i}=\hat{E}_{0}/N_{E}={\rm const}$ .

The remaining phase space is divided along the coordinate lines of the conserved pitch $\Lambda$ . These lines appear parabolic in the $(\alpha,X)$ -plane of panel (b) and horizontal in the $(\Lambda,X)$ -plane of panel (c). Two different procedures are used in the domains above and below the red line in Fig. 6, which represents $\Lambda\hat{B}_{\rm max}=1$ :

•

$0<\Lambda<\hat{B}_{\rm max}^{-1}$ : At $X=X_{\rm min}$ , namely the left vertical axis of panel (b), we define a grid in the pitch angle coordinate $-\pi/2\leq\alpha\leq\pi/2$ . In the present example, it consists of $N_{\alpha}=16$ equally sized cells with $\Delta\alpha_{j}=\pi/N_{\alpha}={\rm const}$ . The circles at $X=X_{\rm min}$ in panel (b) indicate the cell centers $\alpha_{j}$ with $j=1...N_{\alpha}$ . These points $(\alpha_{j},X_{\rm min})$ are the starting points for lines $\Lambda=\cos^{2}\alpha/\hat{B}={\rm const}$ . plotted gray in Fig. 6(b). Along each of these lines, we define a grid in the radial coordinate $X$ . In the present example, it consists of $N_{X}=16$ cells with $\Delta X_{k}=(X_{\rm max}-X_{\rm min})/N_{X}={\rm const}$ . The upper part of panel (c) shows the resulting samples $(\Lambda_{j},X_{k})$ at the cell centers in the $\alpha>0$ portion of panel (b).
•

$B_{\rm max}^{-1}<\Lambda<\hat{B}_{\rm min}^{-1}$ : At $\alpha=0$ , namely the horizontal axis at the center of panel (b), we define a grid in the radial coordinate $X$ . Our example has $N_{X}=16$ equally sized cells. The circles at $\alpha=0$ indicate the cell centers $X_{j}$ with $j=1...N_{X}$ , which form a diagonal ( $\Lambda\hat{B}=1$ ) in the lower part of panel (c). They are the origins of grid lines $\Lambda^{\rm(g)}_{j}={\rm const}$ ., which are plotted gray in panel (b). Here, we choose to have an approximately uniform sample density. For this purpose, we define an auxiliary coordinate $0\leq d_{j}\leq L_{j}$ that measures the distance along a grid line $\Lambda^{\rm(g)}_{j}$ in the $(\hat{\alpha},\hat{X})$ -plane, with $\hat{\alpha}=\alpha/(2\pi)$ and $\hat{X}=X/(X_{\rm max}-X_{\rm min})$ . Along each grid line, we create $N_{j}$ cells with roughly equal sizes $\Delta d_{j,k}=L_{j}/N_{j}\approx{\rm const}$ . The number of cells satisfies $2\leq N_{j}\leq 2N_{X}$ and the cell index $1\leq k\leq N_{j}$ covers positive and negative pitches, so there is at least 1 cell and at most $N_{X}$ cells in each domain, $\alpha\gtrless 0$ , respectively. The lower part of panel (c) shows the resulting samples $(\Lambda_{j},X_{k})$ at the cell centers in the $\alpha>0$ portion of panel (b).

The above procedure samples the phase space of confined GC orbits twice, once on the LFS and once on the HFS. This ensures full coverage of the phase space and eliminates the need to distinguish between orbit types inside the mesh-cutting algorithm. The double-counting will be corrected by a factor $1/2$ in the volume elements derived in Section 6 below.

The color of each sample in Fig. 6(b) identifies the type of an orbit for a given energy. Here, we have chosen $\hat{E}_{i=7}=0.217$ . In the case of alpha particles with $E_{0}=Mv_{0}^{2}/2=3.5\,{\rm MeV}$ , this corresponds to an energy of $1.5\,{\rm MeV}$ . Orbits are divided into two groups: orbits that are trapped in a magnetic mirror, and passing orbits (co- and counter-current). Trapped orbits contain a point where $v_{\parallel}$ changes sign, while passing orbits do not. Each group is further divided into two sub-categories as defined in Table 1. This orbit classification scheme is relatively simple; see Ref. Rome79 for more elaborate discussions. Note that for every tuple $(E,\Lambda)$ , the midplane contains a pair of O-type stagnation points, where co- or counter-passing stagnation orbits are point-like in the poloidal plane. Since these O-type stagnation points alone are insignificant, we have chosen to expand the class of “stagnation orbits” to include all passing orbits that do not encircle the magnetic axis.

	Encircle magnetic axis at $(X,Y)=(0,0)$
	Yes	No
Trapped:	Potato orbits	Banana orbits
Passing:	Circulating orbits	Stagnation orbits

Table 1: Orbit classification.

The classes of potato and stagnation orbits arise from magnetic drifts. The magnitude of the magnetic drift velocity ${\bm{v}}_{\rm dB}$ (see Eq. (23) below) is proportional to the inverse aspect ratio $a/R_{0}$ and the kinetic energy. Furthermore, its $(R,z)$ -component is inversely proportional to the poloidal field $B_{\rm P}\propto I_{\rm p}$ , which tends to divert it into the toroidal direction. This means that magnetic drifts tend to be large for energetic particles in tokamak plasmas with relatively low plasma current and in compact tori. In such cases, the overall fraction of potato and stagnation orbits can be significant.

The mesh in Fig. 6(b) is rather sparse, and our automatic classification algorithm identified only 3 samples as potato orbits (see Fig. 8(a) for a denser mesh and higher energy). The $\alpha>0$ legs of potato orbits are found around the upper rim of the trapped-passing boundary on the $X>0$ side of Fig. 6(b). On the $X<0$ side, they are located slightly below the $\alpha=0$ line.

It is important to note that the trapped-passing (t-p) boundary on the $X<0$ side also lies below the $\alpha=0$ line. This is due to the magnetic drifts and has the implication that the turning points — where the magnetic mirror force causes a sign reversal in $v_{\parallel}$ — are near but not identical to the V-type stagnation points of trapped particle orbits.³³3The discussion below Eq. (19) in Ref. Porcelli94 is related to this. The example in Fig. 7(a) illustrates this for the case of an orbit very close to a t-p boundary.

The t-p boundary refers to a separatrix in the orbit topology; i.e., it marks the topological transition of two V-type stagnation points via an X-type stagnation point as illustrated schematically in the box at the top left of Fig. 7. The marginally trapped orbit in panel (a) has two V-type stagnation points that nearly form an X. Across the t-p boundary, this orbit decomposes into a counter-passing orbit and a trapped orbit, each with a single V-type stagnation point that turns into a smooth part of the orbit contour as one departs from the t-p boundary.

We use the word “stagnation” to indicate that GCs spend a relatively large amount of time in a relatively small region of the poloidal plane. Around O-type stagnation points, this is the case for the entire orbit. In the case of X- and V-type stagnation points as in Fig. 7, GCs spend a large amount of time in a small portion of the orbit contour. Spatial stagnation also means that $v_{\parallel}/v=\sigma_{B}\sqrt{1-\Lambda\hat{B}}$ stagnates as shown in panel (b).

The pitch angle scans of the poloidal and toroidal transit times in panels (c)–(f) of Fig. 7 show true singularities in $\tau_{\rm pol}$ and $\tau_{\rm tor}$ , and a topological discontinuity (see also Fig. 3 of Ref. Zhao97 ). The role of the transit time singularities will be discussed later in Section 6. The discontinuity is a matter of perspective, as it merely reflects how a larger volume element of GC phase space is divided into two components represented by a pair of orbits that were unified on the other side of the t-p boundary. In some sense, the transition is smooth if one views the disconnected orbits as a pair.

Finally, we note that $v_{\parallel}=0$ constitutes a degeneracy in CoM space, since $\Lambda$ and $P_{\zeta}(X)$ cannot be varied independently for fixed $E$ at that location. This is why, geometrically, the parabolic $\Lambda={\rm const}$ . curves have their pole at that location in Fig. 6. Since we are defining our CoM mesh in the midplane, we are able to avoid these degenerate $v_{\parallel}=0$ points entirely by locating them along the $X_{k}$ cell boundary on $(E,\Lambda)={\rm const}$ . lines. All cell centers are, thus, ensured to have $v_{\parallel}\neq 0$ .

5 Equations of motion

As indicated by the red and blue arrows between Steps 2 & 3 in our workflow Fig. 2, the GC phase space samples that are provided by

(a)

the mesh constructed in Section 4, or
(b)

particle codes like OFMC Tani81 ; Tani08 and MEGA Todo98 ; Todo05 (when operated in full- $f$ mode Bierwage13a ; Todo14b )

serve as initial conditions for the GC equations of motion, which we solve once around each orbit.

VisualStart solves the same equations with the same numerical scheme as its current target code MEGA. Following the formalism introduced by Littlejohn Littlejohn83 and reviewed by Cary & Brizard Cary09 , we use the effective electromagnetic fields ${\bm{B}}^{*}\equiv{\bm{\nabla}}\times{\bm{A}}^{*}$ and ${\bm{E}}^{*}\equiv-{\bm{\nabla}}\Phi^{*}-\partial_{t}{\bm{A}}^{*}$ , with the potentials ${\bm{A}}^{*}\equiv{\bm{A}}+\frac{v_{\parallel}}{\Omega}{\bm{B}}$ and $Qe\Phi^{*}\equiv Qe\Phi+\mu B$ . Here, we consider the motion of GCs in an unperturbed magnetic field and in the absence of any electric field ( $\partial_{t}{\bm{B}}=\Phi=0$ ), so the effective fields reduce to

{\bm{B}}^{*}={\bm{B}}+\frac{v_{\parallel}}{\Omega}B{\bm{\nabla}}\times\hat{\bm{b}},\qquad{\bm{E}}^{*}=-\frac{\mu{\bm{\nabla}}B}{Qe},

(20)

and the equations for the GC velocity ${\bm{v}}_{\rm gc}=\dot{\bm{x}}_{\rm gc}$ and its parallel acceleration become

	$\displaystyle\dot{\bm{x}}_{\rm gc}=\;$	$\displaystyle v_{\parallel}\frac{{\bm{B}}^{}}{B_{\parallel}^{}}+\frac{{\bm{E}}^{}\times\hat{\bm{b}}}{B_{\parallel}^{}}={\bm{v}}_{\parallel}^{}+{\bm{v}}_{\rm dB}^{},$		(21)
	$\displaystyle\dot{v}_{\parallel}=\;$	$\displaystyle\frac{\Omega}{B}\frac{{\bm{B}}^{}\cdot{\bm{E}}^{}}{B_{\parallel}^{*}}.$		(22)

$\hat{\bm{b}}={\bm{B}}/B$ is the unit vector along ${\bm{B}}$ , and $B_{\parallel}^{*}\equiv{\bm{B}}^{*}\cdot\hat{\bm{b}}$ is Littlejohn’s GC Jacobian for the noncanonical set of GC coordinates $(\mu,u,\xi,{\bm{x}}_{\rm gc})$ with $u\equiv\hat{\bm{b}}({\bm{x}}_{\rm gc})\cdot\dot{\bm{x}}_{\rm gc}\approx v_{\parallel}$ (see also the discussion on p. 704 of Ref. Cary09 ). The GC velocity is composed of a parallel component (modified by curvature effects) and magnetic drifts due to ${\bm{\nabla}}B$ and curvature ${\bm{\nabla}}\times\hat{\bm{b}}$ :

{\bm{v}}_{\parallel}^{*}=v_{\parallel}\frac{\bm{B}}{B_{\parallel}^{*}},\quad{\bm{v}}_{\rm dB}^{*}=\frac{\mu}{QeB_{\parallel}^{*}}\hat{\bm{b}}\times{\bm{\nabla}}B+\frac{v_{\parallel}^{2}}{\Omega}\frac{B}{B_{\parallel}^{*}}{\bm{\nabla}}\times\hat{\bm{b}}.

(23)

In the normalization of VisualStart (Section 4):

\hat{\bm{v}}_{\parallel}^{*}=\hat{v}_{\parallel}\frac{\hat{\bm{B}}}{\hat{B}_{\parallel}^{*}},\quad\hat{\bm{v}}_{\rm dB}^{*}=\rho_{0}\frac{\hat{\mu}}{\hat{B}_{\parallel}^{*}}\hat{\bm{b}}\times{\bm{\nabla}}\hat{B}+\rho_{0}\frac{\hat{v}_{\parallel}^{2}}{\hat{B}_{\parallel}^{*}}{\bm{\nabla}}\times\hat{\bm{b}},

(24)

with $\rho_{0}=v_{0}/\Omega_{0}$ . The computed GC orbits are recorded in a database, following accuracy checks (B).

The example in Fig. 6, gives a total of about 8000 orbit samples, 52% of which lie inside the region bounded by the red curve ( $\Lambda\hat{B}_{\rm max}\leq 1$ ) for the present setup. About 1100 orbits (14%) are lost as they hit the (artificial) boundary of the simulation domain due to magnetic drifts. Using 3 processes on an 8th Generation Intel CORE i7-8565U processor ( $1.80\,{\rm GHz}\times 8$ ), it took about 35 minutes to compute these 8000 orbits with our simple Matlab code in which there is much room left for optimization. A proper orbit database for the entire plasma volume may require at least 100 times more samples, which would take 150 core-hours with the present implementation. An optimized code in a language such as Fortran, C or Julia is expected to need only a small fraction of this time.⁴⁴4An experimental Julia version of our orbit solver runs 50 times faster. Making use of GPUs, the time required to compute such an orbit database should be reasonably short for routine work.

6 Volume elements and weighting

6.1 Derivation

Having computed the GC orbits using the unperturbed equations of motion in Section 5 starting from the midplane mesh defined in Section 4, the next step is to represent the distribution function $f$ in discretized form using an ensemble of marker particles labeled $n=1...N_{\rm mrk}$ , which sample the orbits in our database. In other words, we seek the Klimontovich representation⁵⁵5A describes in detail how we evaluate Eq. (25).

\displaystyle f({\bm{Z}})\approx\sum_{n=1}^{N_{\rm mrk}}w_{n}\times\delta({\bm{Z}}_{{\rm gc},n}-{\bm{Z}}).

(25)

At this point, ${\bm{Z}}_{{\rm gc},n}$ is the marker position in an arbitrary set of GC coordinates including ignorable angles, for instance $(v_{\perp}^{2},v_{\parallel},\xi,{\bm{x}}_{\rm gc})$ . The Dirac $\delta$ distribution has units of inverse phase space volume $[{\rm m}^{-6}{\rm s}^{3}]$ and may be viewed as a proxy for any sort of particle shape factor needed to map the weights $w_{n}$ onto a discrete mesh. Each marker follows the trajectory of a GC as indicated by the subscript in ${\bm{Z}}_{{\rm gc},n}$ . Its weight factor $w_{n}$ represents a certain number $[\Delta\mathcal{N}]_{n}$ of GCs of physical particles:

	$\displaystyle w_{n}=[\Delta\mathcal{N}]_{n}=$	$\displaystyle\;\Delta\mathcal{V}_{n}\times\frac{\mathcal{N}(\Delta\mathcal{V}_{n})}{\Delta\mathcal{V}_{n}},$
	$\displaystyle=$	$\displaystyle\;\Delta\mathcal{V}_{n}\times f.$		(26)

where $\Delta\mathcal{V}_{n}$ is a 6-D volume element of GC phase space in Cartesian coordinates that is attached to marker $n$ .

Since we are considering only exact equilibrium distributions, the orbit time $\tau$ is an ignorable angle coordinate with period $\tau_{\rm pol}$ and all factors in Eq. (26) are independent of the three angle coordinates $(\tau,\xi,\zeta)$ . These distributions depend only on three CoM coordinates, whose values we will represent by cell indices $(i,j,k)$ . When we load $N_{\tau}$ markers uniformly in time along an unperturbed GC orbit, neighboring markers are indistinguishable. Each marker carries precisely an $N_{\tau}$ ’s fraction of their orbit’s volume $\Delta\mathcal{V}_{\rm orb}$ and weight $W$ :

	$\displaystyle\Delta\mathcal{V}_{n}\rightarrow\Delta\mathcal{V}_{ijkl}=$	$\displaystyle\left[\frac{\Delta\mathcal{V}_{\rm orb}}{N_{\tau}}\right]_{ijk}=\left[\frac{\Delta\mathcal{V}_{\rm orb}\Delta\tau_{l}}{\tau_{\rm pol}}\right]_{ijk},$		(27)
	$\displaystyle w_{n}\rightarrow w_{ijkl}=$	$\displaystyle\left[\frac{W}{N_{\tau}}\right]_{ijk}\quad\text{for}\;l=1,...,N_{\tau,ijk}.$		(28)

As indicated in Eq. (27), the markers also represent equal fractions $\Delta\tau$ of an orbit’s poloidal period $\tau_{\rm pol}$ :

\Delta\tau_{ijkl}=\left[\frac{\tau_{\rm pol}}{N_{\tau}}\right]_{ijk}.

(29)

Recalling from Eqs. (26) that the volume element is eventually applied to a distribution function $f$ , it is clear that Eq. (29) actually implies an orbit time average:

\sum_{l=1}^{N_{\tau,ijk}}\Delta\tau_{ijkl}f({\bm{Z}}_{{\rm gc},ijkl})\approx\oint_{\tau_{\rm pol}}{\rm d}\tau\,f\equiv\tau_{\rm pol}\underbrace{\left<f\right>_{\rm orb}}\limits_{f_{\rm orb}}.

(30)

The average is trivial for true equilibrium distributions, which are independent of $\tau$ , so that $f=f_{\rm orb}$ . However, $f$ in Eq. (30) may also be replaced by an arbitrary quasi-distribution $G$ , so this equation also constitutes a recipe for transforming arbitrary models into true equilibria. In summary, Eq. (26) becomes

w_{ijkl}=\left[\frac{W}{N_{\tau}}\right]_{ijk}=\left[\frac{\Delta\mathcal{V}_{\rm orb}}{N_{\tau}}\times f_{\rm orb}\right]_{ijk}.

(31)

In applications where we merely import particle data computed by another code (blue workflow in Fig. 2), the weights $w_{l}$ are provided with the data. Each sample is interpreted as the initial position of a new GC orbit, so the input weight is interpreted as the weight of an orbit: $[w_{\rm import}]_{l}\rightarrow W_{ijk}$ . This step enforces the equilibrium constraint and the first equality of Eq. (31) yields the marker weights. The volume element $\Delta\mathcal{V}$ and distribution function $f$ do not appear explicitly in such a case.

In applications where we construct a new distribution function from a model or where we reexpress a mesh-based distribution function in the orbit-based representation (red workflow in Fig. 2), it is necessary to determine the GC orbit volume element $[\Delta\mathcal{V}_{\rm orb}]_{ijk}$ for a given mesh in CoM space. In our case, the mesh has the form shown in Fig. 6 and the associated volume elements $[\Delta\mathcal{V}_{\rm orb}]_{ijk}$ in units of $[{\rm m}^{6}{\rm s}^{-3}]$ can be readily obtained when written in terms of the canonical action coordinates $(E,\mu,P_{\zeta})$ associated with the ignorable angles $(\tau,\xi,\zeta)$ , taking advantage of the fact that the Jacobian for a canonical transformation is constant in both space and time (see, for instance, Ref. Porcelli94 ).

With our definition $P_{\zeta}=-\Psi_{\rm P}/B_{0}+v_{\parallel}B_{\zeta}/(\Omega_{0}B)$ , which has units of $(2\pi)^{-1}[{\rm m}^{2}]$ , the transformation between Cartesian coordinates $({\bm{v}},{\bm{x}})$ for the phase space of physical particles and canonical coordinates $(E,\tau,\mu,\xi,P_{\zeta},\zeta)$ for the phase space of GCs has the form

	$\displaystyle{\rm d}^{3}{\bm{v}}{\rm d}^{3}{\bm{x}}=$	$\displaystyle\frac{1}{2}\times\left(\left[\frac{{\rm d}E}{M}\frac{{\rm d}\mu B_{0}}{M}{\rm d}P_{\zeta}{\rm d}\tau{\rm d}\xi{\rm d}\zeta\right]_{\sigma_{\rm HFS}}\right.$
		$\displaystyle\quad\,\left.+\left[\frac{{\rm d}E}{M}\frac{{\rm d}\mu B_{0}}{M}{\rm d}P_{\zeta}{\rm d}\tau{\rm d}\xi{\rm d}\zeta\right]_{\sigma_{\rm LFS}}\right).$		(32)

The index $\sigma=v_{\parallel}J_{\parallel}/|v_{\parallel}J_{\parallel}|$ determines in which domain (co- or counter-going) a volume element lies, irrespective of the orbit type (trapped or passing). Using the fact that canonical volume elements are identical anywhere on a GC orbit $(E,\mu,P_{\zeta})$ , Eq. (32) adds the volume elements at the high- and low-field side (HFS, LFS) midplane crossings and divides the result by 2, in accordance with our double-counting convention in Eq. (8c). Integrating over the angles $\xi$ and $\zeta$ , using the pitch coordinate $\Lambda=\mu B_{0}/E$ instead of the magnetic moment, and applying our normalizations $\hat{v}=v/v_{0}$ and $\hat{E}=E/(Mv_{0}^{2})$ , we obtain

	$\displaystyle(2\pi)^{2}\hat{E}{\rm d}\hat{E}{\rm d}\Lambda\times\left(\left[{\rm d}P_{\zeta}{\rm d}\tau\right]_{\sigma_{\rm HFS}}+\left[{\rm d}P_{\zeta}{\rm d}\tau\right]_{\sigma_{\rm LFS}}\right)\times\frac{v_{0}}{2}$
	$\displaystyle\approx(2\pi)^{2}\hat{E}_{i}\Delta\hat{E}_{i}\Delta\Lambda_{j}\Delta P_{\zeta,k}v_{0}\Delta\tau_{ijkl}\times\frac{1}{2},$		(33)

for our CoM mesh with indices $(i,j,k)$ as defined in Section 4, where the cell index $k$ covers positive and negative signs of $\sigma$ for both the HFS and LFS midplane crossings. It remains to specify the increment $\Delta P_{\zeta}$ .

Since our CoM mesh in Fig. 6 is defined in the midplane, the canonical toroidal momentum $P_{\zeta}$ for fixed $E$ and $\Lambda$ is a function of major radius only, which is here expressed in terms of $X=R-R_{0}$ :

[P_{\zeta}]_{E\Lambda}(X)=-\frac{\Psi_{\rm P}(X)}{B_{0}}+\rho_{0}\frac{\hat{v}_{\parallel}(X)}{B(X)}B_{\zeta}(X).

(34)

Its increment can then be evaluated as

\Delta P_{\zeta,ijk}\approx P_{\zeta}\left(\left.X^{\rm(g)}_{k+1}\right|\hat{E}_{i},\Lambda_{j}\right)-P_{\zeta}\left(\left.X^{(g)}_{k}\right|\hat{E}_{i},\Lambda_{j}\right),

(35)

where $X^{(g)}_{k}$ and $X^{(g)}_{k+1}$ are the grid points adjacent to cell $X_{k}$ . Alternatively, using $v_{\parallel}=\sigma_{B}\sqrt{2E(1-\Lambda\hat{B})}$ , the increment $\Delta P_{\zeta}$ can be expressed in terms of $\Delta X$ as

\Delta P_{\zeta}=\underbrace{\left|-\frac{\Psi_{\rm P}^{\prime}}{B_{0}}-\frac{\rho_{0}}{\hat{v}_{\parallel}}\left(\hat{E}+\frac{\hat{v}_{\parallel}^{2}}{2}\right)\frac{B_{\zeta}B^{\prime}}{B^{2}}+\rho_{0}\frac{\hat{v}_{\parallel}}{B}B_{\zeta}^{\prime}\right|}\limits_{|P_{\zeta}^{\prime}|}\Delta X,

(36)

where radial derivatives like $P_{\zeta}^{\prime}\equiv[\partial_{X}P_{\zeta}]_{E\Lambda}$ are taken along the midplane with fixed $E$ and $\Lambda$ , and they are evaluated at cell centers $X_{k}$ .

In summary, substituting Eq. (29) into (33) and multiplying by the number $N_{\tau}$ of markers used to represent an orbit, we obtain the GC phase space volume element represented by the orbit sample in cell $(E_{i},\Lambda_{j},X_{k})$ in a form that can be readily evaluated numerically:⁶⁶6Eq. (37) differs from the formula for the volume element given in Eq. (40) of Ref. Bierwage12a . We consider the present derivation to be more rigorous. This affects mainly trapped particle orbits, which were effectively absent in previous applications of the code VisualStart.

[\Delta\mathcal{V}_{\rm orb}]_{ijk}=4\pi^{2}v_{0}\hat{E}_{i}\Delta\hat{E}_{i}\Delta\Lambda_{j}\underbrace{[\tau_{\rm pol}\Delta P_{\zeta}]_{ijk}}\limits_{\parbox{51.21504pt}{\centering\scriptsize contains zeros \& singularities\@add@centering}}\times\frac{1}{2}.

(37)

The final factor $1/2$ originates from Eq. (8c) via Eq. (32) and compensates our double-counting of all orbits in the domain $X_{\rm min}\leq X\leq X_{\rm max}$ in Fig. 6(b) via index $k$ .

6.2 Discussion

The formula for the volume element in Eq. (37) may look harmless, but the factor $\tau_{\rm pol}\Delta P_{\zeta}$ should be evaluated with care so as to minimize numerical inaccuracies. We have already seen in Fig. 7(c,d) that $\tau_{\rm pol}$ possesses singularities at t-p boundaries, although not all of them are problematic as we will see shortly. The existence of zeros and singularities in the increment $\Delta P_{\zeta}$ can be readily anticipated from the derivative $|P_{\zeta}^{\prime}|_{E\Lambda}$ in Eq. (36), which can be approximated at leading order as

P_{\zeta}^{\prime}(X)\approx-\frac{\Psi_{\rm P}^{\prime}(X)}{B_{0}}+\frac{\rho_{0}}{\hat{v}_{\parallel}(X)}\left(\hat{E}+\frac{1}{2}\hat{v}_{\parallel}^{2}(X)\right).

(38)

Clearly, $v_{\parallel}=0$ constitutes a singularity and there are at least two points where the terms on the right-hand side cancel, giving $|P_{\zeta}^{\prime}|_{E\Lambda}=0$ . Figure 8 shows an example for $3.5\,{\rm MeV}$ alpha particles. We consider orbits along the black parabolic line in Fig. 8(a), where $\Lambda=0.865$ .

Figure 8(b) shows the form of the increment $\Delta P_{\zeta}(X)$ evaluated using Eq. (35) (solid curves) and Eq. (36) (dotted). The values agree nicely nearly everywhere, except in the vicinity of the $v_{\parallel}=0$ singularity labeled (E).⁷⁷7We expect that the red and blue curves in Fig. 8(b,c) connect smoothly near $X\approx-0.4\,{\rm m}$ if one adds grid points closer to $v_{\parallel}=0$ . A verification exercise similar to that performed later in Section 8.4 showed systematic errors when using Eq. (36). These errors seemed to be effectively absent when $\Delta P_{\zeta}$ is evaluated using Eq. (35), which is therefore our default choice.

The zeros of $\Delta P_{\zeta}(X)$ correspond to stagnation points and imply that these and nearby orbits have a small weight factor $W_{ijk}$ . In the example in Fig. 8 there are three such points: (C) and (F) are the O-type stagnation points in the domains of counter- and co-passing particles; whereas (D) is an X-type stagnation point situated on the separatrix of the t-p boundary.

The zero of $\Delta P_{\zeta}(X)$ at point (D) plays an important role: it cancels one of the three singularities of the poloidal transit time $\tau_{\rm pol}$ in Fig. 8(c). In fact, a glance at the nearby barely trapped orbit in Fig. 8(e) shows that points (B,D,G) actually refer to the same transit time singularity; i.e., the same t-p separatrix. Strictly speaking, our algorithm counts such separatrices thrice, but one of these instances — namely the X-point at (D) — vanishes, owing to $\Delta P_{\zeta}(X)$ being zero at that location. This confirms that all orbits are effectively double-counted, so that the overall factor $1/2$ appearing in Eq. (37) is justified even near the t-p boundary. Of course, on a discrete mesh the cancellation of the third sample at point (D) is unlikely to be perfect, so some amount of numerical inaccuracy must be expected.

The remaining transit time singularities — namely points (B) and (G) in our example — have the consequence that orbits representing those and nearby portions of GC phase space may accumulate a relatively large weight. One could compensate these singularities by demanding that the value of the distribution function $f_{\rm orb}$ vanishes at the t-p boundary; e.g., by including a factor $\tau_{\rm pol}^{-1}$ in $f_{\rm orb}$ . However, this may not be physically meaningful as the following arguments show.

The particles populating the barely trapped orbit in Fig. 8(e) as well as nearby potato and counter-passing orbits spend a relatively large amount of time in a small portion of pitch-position space $(v_{\parallel}/v,R,z)$ near V-type stagnation points. This fact is highlighted by the bunching of the small black circle symbols in Fig. 8(e,f), which represent 16 marker particles distributed uniformly in time. However, unlike the tiny O-type stagnation orbits around points (C) and (F), the large orbits near t-p boundaries like (B) and (G) can actually become densely populated in real plasmas. This becomes evident if one thinks in terms of particle sources: particles deposited on orbits located near a t-p boundary, like our example in Fig. 8(e), will be conveyed to the V-type stagnation points and spend a long time in that region, producing a local spike in the particle density.

Of course, the finite particle supply rate along with instabilities and collisions will prevent the formation of a true singularity in a real plasma. Nevertheless, the above arguments show that, within certain limits, the spikes in the transit times in Fig. 8(c) and resulting spikes in the orbit volume elements (37) and marker weights (31), are physical meaningful.

Such spikes can however cause numerical problems. In some cases, the resulting inaccuracies may be tolerable. In the worst case, the presence of singularities may prevent the attainment of numerical convergence, since a finer mesh will have grid points closer to the singularity. If necessary, one may adjust the weights as mentioned a few paragraphs earlier; e.g., by imposing an upper limit on their values (or their gradients) as in the flux limiter scheme of fluid dynamics. Although this trick may ensure convergence, it is likely to converge to a somewhat inaccurate value. On the other hand, one may argue that the correct value is unknown as it depends on physical mechanisms that have been ignored, such as scattering by collisions and field fluctuations. Perhaps the most elegant solution in the collisionless limit is to use an analytical estimate of a volume element near the singularity, provided that a well-behaved solution exists. Anyhow, we believe that the optimal approach to dealing with singularities depends on the application at hand, so we leave it for future research. It seems meaningful to deal with this problem in the context of smoothing techniques that will be needed to evaluate phase space gradients for instability analyses and may be based on physical scattering mechanisms.

7 Modeling CoM distributions with large drifts

A computed or modeled equilibrium distribution $f_{\rm orb}$ is needed to determine the weight factor in Eq. (31). A model that is a function of CoM only, such as $G_{\rm mdl}(\mathcal{E},\mu,P_{\zeta})$ , would directly yield an exact equilibrium distribution $f_{\rm orb}$ . While simple models may suffice for physics studies at a fundamental level, predictive simulations and simulations used to interpret observations require more realism. It would be convenient to have $G_{\rm mdl}$ given in an explicit, purely analytical form. However, for fast ions with significant magnetic drifts, a large number of control parameters is usually necessary to mimic realistic distribution functions, which tends to make the modeling job cumbersome, at least for non-artificial intelligence.⁸⁸8With a human-manageable amount of control parameters, it can be difficult to design a CoM distribution that has both the desired radial profile and the desired pitch distribution, because both vary simultaneously when one varies parameters controlling the dependence of $f_{\rm mdl}$ on the canonical toroidal angular momentum $P_{\zeta}$ .

For instance, consider the problem of constructing a model for fusion-born alpha particles. Even if one assumes that their pitch angle distribution is uniform at birth time,⁹⁹9Precisely speaking, when an anisotropy exists in the distribution of fusion fuel particles, this anisotropy may also be imprinted on the fusion products due to angular correlations in nuclear reactions (e.g., see Ref. Deutsch51 and Fig. 10 in Ref. Hanson49 ). Here this effect is neglected. strong pitch angle nonuniformities can promptly arise within one poloidal transit time. The loss cone, which contains orbits whose trajectories intersect plasma-facing components, is an obvious example. The magnetic drifts that cause those losses in toroidal geometry also affect the pitch angle distribution of confined orbits. For instance, for fusion-born alpha particles in ITER, the anisotropy associated with drift orbit topology was estimated to be around $10\%$ (see Section 7 of Ref. Salewski18b ). Pitch anisotropies and bump-on-tail structures along the energy axis are also a well-known cause of velocity space instabilities that can be observed as ion cyclotron emissions (ICE) (e.g., see Refs. Sumida19 ; Sumida21 ).

Let us consider a more concrete example. Figure 9 shows the GC orbit of a barely trapped $3.5\,{\rm MeV}$ alpha particle, which traverses both the hot plasma core and the cooler plasma periphery. If one uses a slowly varying function $G_{\rm mdl}(\mathcal{E},\mu,P_{\zeta})$ to model the overall distribution of alpha particles in this range of energies and pitch angles, the density will be similar on the barely trapped orbit shown in Fig. 9 and on the neighboring large potato and small counter-passing orbits indicated by dashed lines, because their coordinates have similar values.¹⁰¹⁰10In our example: $\mu B_{0}/E=0.875...0.892$ and $\hat{P}_{\zeta}=0.215...0.219$ with normalized $\hat{P}_{\zeta}\equiv\psi_{\rm P}-\hat{v}_{\parallel}\frac{B_{\zeta}}{B}\frac{\rho_{0}B_{0}}{\Psi_{{\rm P},a}-\Psi_{{\rm P},0}}$ and $0\leq\psi_{\rm P}\leq 1$ . In other words, the particle density will be similar on both sides of the t-p boundary.

However, this may not be an accurate representation of a real distribution of fusion-born alphas, in particular when the fusion reactivity varies substantially in different portions of the orbit. To illustrate this, the $3.5\,{\rm MeV}$ alpha particle orbit in Fig. 9 has been underlaid with the contours of the sharply peaked MHD pressure field from of our working example (cf., Fig. 4), which is characteristic for the bulk component of a JET plasma where 3-ion RF heating was applied near the magnetic axis Nocente20 ; Kazakov21 . In such a case, fusion-born alphas from the plasma center will populate the entire orbit shown in Fig. 9, albeit in a nonuniform manner as indicated by the small black circles. Just across the t-p boundary, however, most of the newly born alphas will be concentrated on the small counter-passing orbit, while the large potato orbit would be populated only sparsely, since it is topologically disconnected from the particle source. In principle, this can produce a steep gradient in the distribution function $f_{\rm orb}(\mathcal{E},\mu,P_{\zeta})$ across the t-p boundary. In reality, collisions or instabilities can be expected to have a smoothing effect, but depending on the relevant time scales, a strong nonuniformity may be sustained. Such information may be necessary for interpreting experimentally observed signals or for reproducing them in simulations, and it is advantageous to have the ability to construct suitable models.

A viable recipe for doing so can be found in Eq. (30). Such an orbit time average can be readily realized using our orbit database and constitutes a physically meaningful way to convert an arbitrary model function $G_{\rm mdl}({\bm{Z}})$ into a true equilibrium distribution $f_{\rm orb}$ :

	$\displaystyle f_{\rm orb}=$	$\displaystyle\;\left<G_{\rm mdl}(\rho_{\rm P})\right>_{\rm orb}\equiv\frac{1}{\tau_{\rm orb}}\oint{\rm d}\tau\,G_{\rm mdl}({\bm{Z}}(\tau))$
	$\displaystyle\approx$	$\displaystyle\;f_{\rm orb}(i,j,k)=\frac{1}{N_{\tau,ijk}}\sum_{l=1}^{N_{\tau,ijk}}G_{\rm mdl}({\bm{Z}}_{{\rm gc},ijkl});$		(39)

where ${\bm{Z}}_{{\rm gc},ijkl}$ is the position of marker $l$ on orbit $(i,j,k)$ . For instance, in the example discussed in the previous paragraphs, one may let $G_{\rm mdl}\propto P(\Psi_{\rm P})$ , so that the MHD pressure field in Fig. 9 acts as a particle source.

This approach can be viewed as taking the path of “self-organization”, where we leave it to the equations of motion to determine a distribution function that is consistent with the magnetic geometry and system boundaries at hand, and the constraints encoded in $G_{\rm mdl}$ determine the particular spatial profile and velocity distribution for a certain case. This idea also underlies Monte-Carlo codes that follow GC orbits many times around the torus in the presence of sources and collisions. We mimic a part of this process with reduced computational effort: the GC orbits in our database are followed for only one poloidal turn in the absence of any perturbations (field fluctuations and collisions).

Form of	Sampling parameters				$E_{0}=35\,{\rm keV}$		$E_{0}=350\,{\rm keV}$		$E_{0}=3.5\,{\rm MeV}$
$E$ -distrib.	$N_{E}$	$N_{\alpha}$	$N_{X}$	$N_{\tau}$	$N_{\rm orb}$	$N_{\rm mrk}$	$N_{\rm orb}$	$N_{\rm mrk}$	$N_{\rm orb}$	$N_{\rm mrk}$
Mono- $E_{0}$	$1$	$2\times 64$	$2\times 64$	$128\frac{L_{\rm orb}}{L_{\rm bnd}}$ (min.: 16)	33,335	1.88M	31,841	1.70M	27,684	1.29M
Flat- $E\leq E_{0}$	$24$	$2\times 24$	$2\times 24$	$48\frac{L_{\rm orb}}{L_{\rm bnd}}$ (min.: 16)	112,483	2.74M	—	—	99,308	2.18M

Table 2: Phase space sampling parameters in our five test cases.

N_{E}

N_{\alpha}

N_{X}

determine the number of cells in

(E,\alpha,X)

-space defined in Section 4. The number

N_{\tau}

of markers sampling an orbit is (with a lower bound of 16) chosen to be proportional to the orbit length

L_{\rm orb}

, the reference value being the number of markers per plasma boundary length, here

L_{\rm bnd}\approx 7.50\,{\rm m}

. There are

N_{\rm mrk}

markers sampling

N_{\rm orb}

confined orbits.

$E\,[{\rm keV]}$	Circ.	Stagn.	Banana	Potato	Lost $\approx$
$35$	25,780	359	7,158	38	780
$350$	24,430	1,066	6,191	154	2,280
$3500$	20,531	3,165	3,505	483	6,430
$0...35$	87,368	612	24,456	47	1,710
$0...3500$	75,178	7,208	15,861	1,061	14,880

Table 3: Partition of orbit types (circulating, stagnation, banana, potato, lost) as defined in Section 4. These numbers are not weighted by particle densities, so they only characterize the partition of orbit samples on our mesh. The actual prompt losses of real particles are negligible in our large plasma with centrally peaked density profile.

8 Application and verification

In this section, we demonstrate how the methods that we presented in this paper can be applied to a practical example. Section 8.1 describes the setup and how we construct — from a simple model $G_{\rm mdl}$ — an equilibrium distribution $f_{\rm orb}$ of alpha particles that are assumed to be born near the center of a large tokamak plasma with an isotropic distribution of pitch angles at birth time. The results are discussed in Sections 8.2 and 8.3. In the final Section 8.4, we demonstrate how we use our orbit database to convert a binned 4-D distribution function $f^{(0)}(E,\lambda,R,z)$ (or any other coordinates) back to an orbit-based representation $f_{\rm orb}^{(1)}$ in CoM space. The result is binned again to give $f^{(1)}$ for comparison with the original $f^{(0)}$ . The procedure is then iterated one more time, yielding $f^{(2)}$ , in order to reveal systematic errors. The contents of this section may be summarized in diagrammatic form like this:

\begin{array}[]{r@{\hspace{0.2cm}}c@{\hspace{0.cm}}c@{\hspace{0.2cm}}l}\parbox{28.45274pt}{\centering\text{\scriptsize Modeling}\@add@centering}\hskip 5.69046pt&\parbox{56.9055pt}{\centering\text{\scriptsize CoM space}\@add@centering}&\hfil\hskip 5.69046pt&\parbox{56.9055pt}{\centering\text{\scriptsize Arbitrary coordinates}\@add@centering}\\ \framebox{$G_{\rm mdl}$}\rightarrow\hskip 5.69046pt&\framebox{$f_{\rm orb}^{(0)}(i,j,k)$}&\rightarrow\hfil\hskip 5.69046pt&\framebox{$f^{(0)}(E,\lambda,R,z)$}\\ \hskip 5.69046pt&&\swarrow\hfil\hskip 5.69046pt&\;\;\rotatebox{90.0}{$\approx$}\;?\;\;\text{\scriptsize 1st verification}\\ \hskip 5.69046pt&\framebox{$f_{\rm orb}^{(1)}(i,j,k)$}&\rightarrow\hfil\hskip 5.69046pt&\framebox{$f^{(1)}(E,\lambda,R,z)$}\\ \hskip 5.69046pt&&\swarrow\hfil\hskip 5.69046pt&\;\;\rotatebox{90.0}{$\approx$}\;?\;\;\text{\scriptsize 2nd verification}\\ \hskip 5.69046pt&\framebox{$f_{\rm orb}^{(2)}(i,j,k)$}&\rightarrow\hfil\hskip 5.69046pt&\framebox{$f^{(2)}(E,\lambda,R,z)$}\end{array}

(40)

Strictly speaking, we transform only between 2-D and 3-D spaces, $(A_{j},P_{\zeta}(X_{k}))\leftrightarrow(\lambda,R,z)$ , since both sets share the energy coordinate $E$ , which is a conserved quantity, so the mapping in that direction is identical. We will work here with an arbitrary normalization (henceforth indicated by hats), so the values of $\hat{f}_{\rm orb}$ of $\hat{f}$ must match but bear no particular meaning.

Our distribution functions will be inspected here in the same way as an experimental diagnostician would: measure spatial distributions by integrating over certain portions of velocity space, and measure velocity distributions by integrating over certain portions of position space. The required binning operations and Jacobians are described in A.

8.1 Setup and parameters

In order to highlight how toroidal geometry affects the form of the constructed equilibrium distribution function via the mirror force and magnetic drifts, we start from a quasi-distribution $G_{\rm mdl}$ that is independent of pitch and has a relatively sharp peak in minor radius:

G_{\rm mdl}(\rho_{\rm P})=C\times\hat{n}_{\rm mdl}(\rho_{\rm P})\times H_{\rm mdl}(E),

(41)

where $C$ is a normalization constant that is arbitrary here.

The radial density profile $\hat{n}_{\rm mdl}(\rho_{\rm P})$ is shown as a solid black line in Fig. 10(a). Although this is primarily meant to be a toy model for testing our algorithms with a localized particle source and large magnetic drifts, we note that this profile is physically feasible if one considers scenarios with a sharply peaked ion temperature profile as produced in 3-ion-heated JET plasmas Nocente20 ; Kazakov21 .

The energy dependence of particle distributions in a real plasma is governed by processes such as collisional drag and instabilities (e.g., see Fig. 4 in Ref. Bierwage18 ). These processes are not considered here. We use energy only to exemplify the effects of magnetic drifts, which increase with $E$ ; or, rather, with its square root $\sqrt{E}$ , because the effect of magnetic drifts is roughly proportional to the ratio $v_{\rm dB}/v$ of the drift velocity to the particle velocity. For this purpose, we use two types of distributions: one mono-energetic with $E=E_{0}$ , and one that is flat in the range $0\leq E\leq E_{0}$ and zero beyond:

\displaystyle H_{\rm mdl}(E)

\displaystyle=\left\{\begin{array}[]{l}\delta(E-E_{0})\quad\;\;\;=\text{ mono-energetic},\\ \frac{1-\int_{0}^{E}{\rm d}E^{\prime}\,\delta(E^{\prime}-E_{0})}{E_{0}}=\frac{1}{E_{0}}\,{\rm for}\,0\leq E\leq E_{0},\end{array}\right.

(44)

where $\delta$ is the Dirac delta distribution.

Table 2 shows the parameters for five test cases that we will consider here: three mono-energetic alpha distributions with $E=E_{0}=(35,350,3500)\,{\rm keV}$ and two distributions that are uniform in the energy range $0\leq E\leq E_{0}$ with $E_{0}=(35,3500)\,{\rm keV}$ . Table 3 summarizes the numbers of circulating, stagnation, banana, potato and lost orbit samples (CoM mesh points) in each case. Panels (b) and (c) of Fig. 10 show how the contours of a banana orbit (A) and a circulating orbit (B) change due to magnetic drifts when the kinetic energy is increased from $35\,{\rm keV}$ to $3.5\,{\rm MeV}$ .

While the model in Eq. (41) has a simple appearance, the actual form of the alpha particle distributions we construct here is determined by the combination of Eq. (41) and the equations of GC motion (21) and (22), whose solutions are folded into the simple model via the integral (39). In other words, the geometric effects are captured by the GC orbits in our database, whereas the model $G_{\rm mdl}$ effectively describes the structure of the particle sources that populate those orbits, and the integral (39) combines all that information.

The geometric effects that we expect to see in this setup may be summarized as follows. These properties of a toroidally confined and radially bounded plasma are well-known and are listed here only for completeness and to simplify the discussion of the results:

•

The peripheral plasma will be populated by particles on orbits that also pass through the core as in Fig. 10(c). Such orbits are found in a certain range of pitch angles around the t-p boundary.
•

Flux surface-averaged profiles tend to be smoothed on the scale length of the magnetic drifts. The gyroradii have a similar effect and may be included via satellite particles Bierwage16c (not done here).
•

Magnetic drifts shift co-passing orbits outward in $R$ and counter-passing orbits inward.
•

The mirror force confines some particles to the low-field side of the plasma and causes them to accumulate near V-type stagnation points.
•

The energy- and pitch-dependence of the above-mentioned spatial nonuniformities gives rise to net flows both poloidally (due to radial nonuniformity) and toroidally (due to a local imbalance between co- and counter-going particles).

8.2 Properties of mono-energetic distributions

Figure 10(a) shows the computed density profiles $\hat{n}(\rho_{\rm P})$ on a mesh that is uniformly spaced in normalized poloidal flux $\psi_{\rm P}=\rho_{\rm P}^{2}$ . The profile for $35\,{\rm keV}$ alphas is similar to the model $\hat{n}_{\rm mdl}(\rho_{\rm P})$ . The small discrepancy at $\rho_{\rm P}=0$ can be eliminated with a finer mesh as shown in the inset. The magnetic drifts increase with increasing energy, causing the profile $\hat{n}(\rho_{\rm P})$ to broaden and the central value $\hat{n}(0)$ to drop. The total number of particles in the plasma is nearly the same in all three cases, since prompt losses are relatively small: less than $1\%$ for $350\,{\rm keV}$ and about $2.5\%$ for $3.5\,{\rm MeV}$ .

Figures 11 and 12 show the form of various moments of the alpha distributions with kinetic energies $35\,{\rm keV}$ and $3.5\,{\rm MeV}$ in the poloidal plane $(X,Y)$ :

•

Particle density (a): Even for $3.5\,{\rm MeV}$ , our model yields an alpha density that is peaked at the axis, but the peak is wider than in the $35\,{\rm keV}$ case, especially in the horizontal ( $X=R-R_{0}$ ) direction.
•

Markers/cell (b): Since we have set a lower bound of 16 markers per orbit, the smallest stagnation orbits stand out in the plots showing the number of markers per cell, while the surroundings are fairly uniform. At $35\,{\rm keV}$ , the co- and counter-passing stagnation points effectively coincide with the magnetic axis. At $3.5\,{\rm MeV}$ , their mean separation along $X$ is about $0.14\,{\rm m}$ (cf. Fig. 8(b)).
•

Temperature (c): With the density nonuniformity divided out, the $3.5\,{\rm MeV}$ temperature field clearly shows that our alpha particle distribution effectively consists of two partially overlapping co- and counter-going mono-energetic beams, which are shifted radially out- and inward, respectively.
•

Flows (d–g): The relative shift of the co- and counter-going beams causes net local currents. The poloidal current field (consisting of $j_{R}$ and $j_{z}$ ) forms two vortices, one on the high-field side (HFS) of the magnetic axis rotating clockwise, and one on the low-field side (LFS) rotating counter-clockwise. Toroidal and parallel currents $j_{\zeta}$ and $j_{\parallel}$ are very similar here, since the magnetic pitch is relatively uniform ( $q\approx 1$ ) in the plasma center. The structure and spatial extent of the currents is similar at lower and higher energies, but the magnitude should be proportional to $\sqrt{E}$ and, thus, differ by about a factor of 10. This seems to be the case here if one accounts for the fact that the $35\,{\rm keV}$ case has a smaller signal-to-noise ratio.

Figure 13 shows the pitch angle distribution of our mono-energetic alpha particles. The spatially integrated distributions $\hat{\nu}(\alpha|E_{0})$ in panel (a) show that a pitch nonuniformity exists around $\alpha=0$ , up to intermediate pitch angles, $|\alpha|\lesssim 0.2\pi$ . The deviation from the mean increases from about $10\%$ for $35\,{\rm keV}$ to about $50\%$ for $3.5\,{\rm MeV}$ . At larger pitch angles, $|\alpha|\gtrsim 0.2\pi$ , the distributions are flat (except for aliasing noise, to be discussed later). More detailed views of the nonuniformities in the region $|\alpha|\lesssim 0.2\pi$ of the cases with $E_{0}=35\,{\rm keV}$ and $3.5\,{\rm MeV}$ are provided in panels (b) and (c), where we plot the contours of $\hat{\nu}(X,\alpha|E_{0})$ to show the radial dependence of the pitch angle distribution in the range $-0.48\,{\rm m}<X<0.48\,{\rm m}$ . Here, the summation over markers was performed over the $Y$ coordinate for radial bins of size $\Delta X=0.04\,{\rm m}$ .

The structures seen in Fig. 13 can be readily explained on the basis of our model function $G_{\rm mdl}$ , the effects of toroidal geometry (mirror force and magnetic drifts), and the loss boundary. Rather than saying anything new, the purpose of the following discussion is to show that our procedures work reliably and produce the expected results. We focus on the $3.5\,{\rm MeV}$ case.

Figure 13(c) clearly shows the opposite radial shift of co- and counter-going particles due to their opposite magnetic drifts: the peak in the region $\alpha>0$ is shifted towards the right ( $X>0$ ) and the peak in the region $\alpha<0$ is shifted to the left ( $X<0$ ). The wedge-shaped structure of the pitch anisotropy around $\alpha=0$ represents the widening of the domain of trapped particles towards the low-field side as expected from Fig. 8(a).

In order to explain the structure of the pitch anisotropy in the region $|\alpha|\lesssim 0.2\pi$ , it is helpful to consider the trajectories of a few GC orbits in that region. For this purpose, Fig. 13(d) shows the GC orbits of four simulation particles with energy $E=3.5\,{\rm MeV}$ . Two of them (solid lines) are launched near the midplane ( $Y_{\rm start}\approx 0$ ): the small orange orbit starts from $(X,\alpha)_{\rm start}=(0.05\,{\rm m},0.09\pi)$ , and the blue orbit starts from $(X,\alpha)_{\rm start}=(-0.05\,{\rm m},-0.09\pi)$ , as indicated by the crosses in panels (c) and (d). The blue orbit is a large potato orbit that passes both near the center of the plasma and through the periphery. The small orange orbit is a stagnation orbit. Now recall Eq. (39), where we defined the orbit distribution $f_{\rm orb}$ to be the orbit time-average of the model $G_{\rm mdl}(\rho_{\rm P})$ that represents our alpha particle birth profile. Since the latter is sharply peaked near the plasma center, the values of $G_{\rm mdl}(\rho_{\rm P})$ can differ substantially on the above-mentioned pair of orbits.

Generally speaking, orbits that spend most of their time inside the peak of the source profile $G_{\rm mdl}(\rho_{\rm P})$ are densely populated with physical particles, so our GC markers representing them have larger weights $w_{l}$ . Orbits that spend most of the time outside are sparsely populated. The mirror force and magnetic drifts have the consequence that particles on counter-passing orbits spend less time inside the central peak of $G_{\rm mdl}(\rho_{\rm P})$ than co-passing orbits. In addition, the losses of co- and counter-going particles are asymmetric, especially in the vicinity of the trapped-passing boundary. These effects are manifested in the minima and maxima of $\hat{\nu}(\alpha|E_{0})$ and $\hat{\nu}(X,\alpha|E_{0})$ in Fig. 13, which are most pronounced in the high-energy case due to its larger drifts.

Note that the distribution $\hat{\nu}(X,\alpha|E_{0})$ in Fig. 13(c) was computed by integrating along the vertical coordinate $Y$ , so it also includes particles on orbits like those shown by dashed lines in panel (d). The dashed orbits were launched from $(X,\alpha)_{\rm start}=(0.05\,{\rm m},\pm 0.09\pi)$ , but more than half a meter above the midplane, from $Y_{\rm start}=0.56\,{\rm m}$ , as indicated by pluses in panels (c) and (d). The orbit starting with a positive pitch (orange) passes through the central plasma, where it acquires a relatively large weight due to the centrally peaked source profile $G_{\rm mdl}(\rho_{\rm P})$ . In contrast, the orbit starting with a negative pitch travels outward and nearly hits the boundary (light green curve).

In fact, the entire peak of $\hat{\nu}(\alpha|E_{0})$ around $\alpha/\pi\approx+0.1$ can be explained in that way. We have verified this with a test case where the plasma boundary had been placed at $\psi_{\rm P}=0.25$ ( $\rho_{\rm P}=0.5$ ) instead of $\psi_{\rm P}=1$ , so that orbits are constrained to the inner 25% of the present flux space Bierwage22b . In that case, the plasma boundary is located around $X\approx 0.38\,{\rm m}$ at the height of the outer midplane, so that the three large orbits in Fig. 13(d) are all lost. The resulting pitch angle profile is shown as a black dotted line in Fig. 13(a). It has a deeper and wider loss cone, and there is no peak at $\alpha/\pi\approx+0.1$ .

8.3 Properties of flat-energy distributions

Realistic alpha particle distributions are nonuniform in energy. The reason for our use of a model that is flat in energy is that the $E$ -dependence of geometrically induced nonuniformities can be easily seen, and it is also easier to verify that our method works as expected.

The densities and current fields are very similar to those in the mono-energetic cases shown in Figs. 11 and 12 above. A few different features can be observed in the high energy case, $E=0...3.5\,{\rm MeV}$ , whose moments are shown in Fig. 14. The presence of particles with a wide range of energies can be inferred primarily from the blurring of the stagnation points in panel (b), and from the nonuniformity of the temperature field in (c). The latter resembles somewhat the letters “IO” as fast particles dominate outside the central density peak; especially the co-passing ones (forming the “O” shifted to the right) for the same reasons as discussed in the last paragraphs of Section 8.2 above.

Figure 15 shows the velocity distributions integrated over position space. One can see that, apart from noise, our result has the desired uniformity in both energy and pitch in the domains populated by passing particles. As discussed in Section 8.2 above, the pitch anisotropy around $\alpha=0$ is induced by toroidal geometry (mirror force, magnetic drifts and resulting losses) in combination with the centrally peaked source profile. The anisotropy clearly increases with increasing energy.

Our orbit database is relatively sparse ( $N_{E}\times N_{\alpha}\times N_{X}=24\times 48\times 48$ , see Table 2), and the bins used in Figs. 14 and 15 have effectively the same resolution: $N_{E}\times N_{\alpha}\times N_{Rz}=24\times 48\times 48^{2}$ . The binned results are smooth along energy, because the orbit samples coincide with the bins. Most of the noise in Fig. 15 is in the pitches $\alpha$ and $\Lambda$ , where the distribution of bins is similar but not identical to the distribution of orbit samples (CoM space). In other words, the noise in the range $0.3\lesssim|\alpha|/\pi\lesssim 0.5$ and $\Lambda\lesssim 0.7$ is effectively an aliasing effect. The noise is largely suppressed in the poloidal plane $(X,Y)$ in Fig. 14, and in $(\hat{v}_{\parallel},\hat{\mu})$ -space shown in the central column of Fig. 15, since the bins in these coordinates have a shape very different from the mesh we used to sample the CoM space (cf. Figs. 6 and 8).

Figure 15 also highlights potential difficulties in attempts to convert the binned distributions from one set of coordinates to another. In particular, the small region around the origin $(\mu,v_{\parallel})=(0,0)$ in the central column of Fig. 15 is expanded to the full length of an axis in $(E,\alpha)$ and $(\Lambda,\mathcal{E})$ spaces on the left and right sides of Fig. 15, so some features of the distribution functions are likely to get lost in direct conversions unless a polar mesh is used in $(\mu,v_{\parallel})$ coordinates, which is equivalent to using another set of coordinates, such as $(E,\alpha)$ . We emphasize that we have not performed direct conversions between binned distributions in Fig. 15. They were all obtained by binning, on the respective mesh, the marker particles of the orbit-based representation $f_{\rm orb}$ .

The binning in Figs. 10–15 was performed using a 1st-order (linear) interpolation scheme as commonly used in particle-in-cell (PIC) codes, where each marker is assumed to have the shape of a top-hat function whose width is equal to the bin size. Naturally, the noise is larger with a 0th-order (binary) binning scheme, where markers resemble $\delta$ functions. The difference between 1st- and 0th-order binning can be seen by comparing panels (a) and (b) of Fig. 16, where we show the spatially integrated velocity distribution $\hat{\nu}(E,\lambda)$ with $\lambda=v_{\parallel}/v=\sin^{-1}(\alpha)$ . Note also that the result in Fig. 16(a) appears smoother than that in 15(d) simply by using $\lambda$ instead of $\alpha$ , because this change in coordinates reduces the above-mentioned aliasing effect. Another difference is that Fig. 16 was obtained with 10 times more markers per orbit, but this does not have a strong influence on the aliasing in the domain of deeply passing particles, since those markers rarely cross the bin boundaries.

8.4 Verification of reversibility and accuracy of coordinate conversions via an orbit-based representation

The results discussed in the previous sections were obtained by binning the marker weights of our orbit-based representation in CoM space to a mesh in arbitrary coordinates. The reverse transformation is also easily performed as we will demonstrate in this section. On platforms like ITER IMAS, we expect this to be an essential part of the transformation toolbox, especially if an orbit-based representation is used as a reference. Here, this operation serves us primarily as a verification exercise to test with what degree of accuracy we can perform and reverse coordinate conversions without applying any smoothing algorithms.

We consider the flat-energy case with $E_{0}=3.5\,{\rm MeV}$ and follow the procedure that was outlined in Eq. (40). Let us now define the relevant steps in detail:

1.

Using Eqs. (39) and (41), we construct a distribution function $f_{\rm orb}$ in the orbit-based representation:

$G_{\rm mdl}\rightarrow f_{\rm orb}(i,j,k).$ (45)

The CoM mesh that underlies our orbit database has the same number of cells as specified in Table 2: $N_{E}\times(N_{\alpha}\times N_{X})=24\times 48^{2}$ . The number of markers per boundary length is increased by a factor 10 to $N_{\tau}=480\times L_{\rm orb}/L_{\rm bnd}$ (min. 16), giving a total of $N_{\rm mrk}=18\,{\rm M}$ markers.

The orbit-based CoM distribution $f_{\rm orb}$ is binned on a uniform 4-D mesh in the widely used coordinates $(E,\lambda,R,z)$ . In compact form, this operation can be written as

f_{\rm orb}^{(\eta)}(i,j,k)\rightarrow f^{(\eta)}(E,\lambda,R,z)

(46)

where $\eta=0,1,2$ is the iteration number. The initial result $f^{(0)}$ is referred to as the “input” distribution for the subsequent back-and-forth transformations. More precisely, the binning operation we perform is given by Eq. (53) in the form

	$\displaystyle f(E_{i},\lambda_{j},R_{k},z_{m})=\$
	$\displaystyle\sum_{n=1}^{\rm N_{\rm mrk}}\frac{w_{n}\times\Pi({\bm{Z}}_{{\rm gc},n}-{\bm{Z}}_{ijkm})}{(2\pi)^{2}\left[\frac{B^{*}_{\parallel}}{B}\hat{v}R\right]_{n}\Delta E_{i}\Delta\lambda_{j}\Delta R_{k}\Delta z_{m}}.$		(47)

where the shape factor $0\leq\Pi\leq 1$ is a top-hat function, which means that the markers resemble $\delta$ functions. That is, we use the 0th-order binning scheme as in Fig. 16(b) and the number of bins (= cells) is $N_{E}\times N_{\lambda}\times(N_{R}\times N_{z})=24\times 24\times 48^{2}$ . The pitch resolution was reduced from $N_{\lambda}=48$ to $24$ in order to reduce noise and aliasing (cf. Section 8.3).

For each orbit in the database, we compute new weight factors $W_{ijk}^{(\eta+1)}$ by integrating $f^{(\eta)}(E,\lambda,R,z)$ in time $\tau$ along the orbit contour as

W_{ijk}^{(\eta+1)}=\left[\Delta\mathcal{V}_{\rm orb}\sum_{l=1}^{N_{\tau}}f^{(\eta)}(E,\lambda(\tau_{l}),R(\tau_{l}),z(\tau_{l}))\right]_{ijk},

(48)

which gives a new orbit-based distribution:

f^{(\eta)}(E,\lambda,R,z)\rightarrow f_{\rm orb}^{(\eta+1)}(i,j,k).

(49)

4.

Repeat step 2 and compare the result with the previous iteration:

$f^{(\eta+1)}\stackrel{{\scriptstyle?}}{{\approx}}f^{(\eta)}.$ (50)

We have then iterated steps 3 and 4 one more time, looking for increasing deviations that may hint at systematic errors. The results are summarized in Figs. 17 and 18.

Figure 17(a) shows that the overall shape of the radial density profile is reproduced well. A closer inspection of the relative errors in Fig. 17(b) reveals that the first iteration changes the particle density by about $+2\%$ in the densely populated region $r/a\lesssim 0.3$ . Around $r/a\approx 0.7$ the increase is about $+10\%$ , and near the boundary we have a reduction by about $-20\%$ . The relative error roughly doubles in the second iteration, indicating that it is not random but systematic. As indicated in the legend of Fig. 17(a), the total number of GCs $\mathcal{N}$ increases by $1.9\%$ in the first and $2.4\%$ in the second iteration, where the error accumulates to $4.3\%$ .

Figure 17(b) shows that our procedure has the tendency to increase the particle density in the plasma interior and reduce it near the plasma boundary. We have not investigated why the systematic errors have this particular form and it may, in fact, depend on the particular magnetic configuration at hand. However, the larger errors in the outer region of the plasma are likely related to the fact that, in our test case, this region is populated almost exclusively by particles on large orbits near the trapped-passing boundary (cf. Fig. 10(c)), where the poloidal transit time $\tau_{\rm pol}$ exhibits singular behavior (see Fig. 8) and inaccuracies can be expected. In the present case, the magnitude of these errors can be reduced by increasing the number of phase space samples as shown in Appendix B.4.

Figure 18 shows the results for the velocity space distributions $\hat{\nu}^{(\eta)}(E,\alpha)$ and $\hat{\nu}^{(\eta)}(E,\lambda)$ for iterations $\eta=0,1,2$ . No significant differences can be seen in the energy direction, which is to be expected since it is a conserved quantity and our energy bins are identical in both representations. But even along the pitch coordinates $\alpha$ and $\lambda$ , excellent agreement is obtained after dividing out the systematic increase of $\mathcal{N}$ . The noise at large pitch values retains its form as expected for an aliasing effect, and its magnitude increases systematically from one iteration to the next. The minimum associated with the loss cone tends to deepen with each iteration, which is also to be expected due to the irreversibility of losses. As was noted in connection with Fig. 16, we see again in Fig. 18 that binning on a uniform mesh in $\lambda$ instead of $\alpha$ reduces aliasing and, thus, yields smoother results in the deeply passing domain ( $|\alpha|\gtrsim 0.3\pi$ or $|\lambda|\gtrsim 0.8$ ).

Before proceeding to the conclusion of this paper, we shall add a few encouraging words for readers who may find the results in this section somewhat unsatisfactory. First, it should be noted that we have chosen a fairly challenging setup for our application example: Our density peak is localized deeply in the plasma core, where derivatives of the poloidal magnetic flux $\Psi_{\rm P}$ are small so that $P_{\zeta}$ is dominated by drift terms. This means the so-called “non-standard” orbits are actually prevalent here. Second, we have sampled the CoM space only relatively sparsely (Table 2) and used no smoothing procedure whatsoever. Even the binning of particles was done using a 0th-order scheme, which is why the noise in Fig. 18 is larger than in Fig. 13.

In other words, one may say the we have been looking here at a worst-case scenario with the least degree of sophistication. This should be taken in a positive way, because it means that there is much room for improvement even with well-established techniques. Higher resolution in CoM space, higher-order binning schemes and other smoothing techniques are obvious solutions to reduce noise to a level where the distribution becomes sufficiently smooth to measure local phase space gradients of a particle distribution. The application of such techniques in combination with the orbit-based representation method described here is left for future work.

9 Summary and conclusion

This paper was motivated by the need to model and process distributions of charged particles in tokamak plasmas in a unified framework like the ITER Integrated Modelling & Analysis Suite (IMAS), which has to address the needs of a diverse community of users dealing with experimental and numerical data. One of the tasks IMAS will need to perform is the conversion of distribution functions between various coordinates. Our proposal is to construct a database of unperturbed GC drift orbits and use it as a basis for representing distribution functions in constants-of-motion (CoM) space. A working prototype of a concrete workflow as implemented in the code VisualStart BierwageVSTART12 ; Bierwage12a is described here.

It should be noted that dedicated IMAS data structures for storing GC orbit properties and quantities along GC trajectories of marker distributions have recently been added into the ‘distributions’ IDS (Interface Data Structures) ITER_IDS , making IMAS readily compatible with the representation and methods proposed here.

An orbit-based representation maps any input distribution to an exact equilibrium that is consistent with the specified magnetic field configuration. While the result may generally differ from the original input, the equilibrium constraint ensures that all subsequent mappings from 3-D CoM orbit space back to arbitrary 4-D sets of coordinates become unique and straightforward. In fact, further operations like visualization or exporting are reduced to binning marker particles on user-defined meshes in arbitrary coordinates and arbitrary subspaces.

A new orbit database needs to be set up for each MHD equilibrium, and one may even have to set up one for different energy ranges (say, $1...\text{tens}\,{\rm keV}$ , and $0.1...\text{few}\,{\rm MeV}$ ). The required computational effort should be manageable if the code is optimized and capable of exploiting state-of-the art computing technology. Once the database has been constructed, it can be used to represent and model various populations of charged particles that may be present in the same plasma, such as beam ions, RF-heated ions and fusion products. One orbit database may be shared between particle species with identical charge-to-mass ratio $Qe/M$ by scaling their energies.

The time spent on preparing the orbit database is also justified by the fact that the acquired information can be used for other important tasks, such as resonance analyses and modeling, which must also be part of the IMAS toolbox. The degrees of freedom for modeling are vast, since each orbit can be assigned an independent weight. Rather than prescribing parametric functions of CoM — whose limitations were briefly discussed in Section 7 — it is worth exploring unconventional ways that make better use of the degrees of freedom available and take advantage of the fact that our orbit-based representation automatically enforces the equilibrium condition and physical constraints associated with toroidal geometry; namely, the mirror force, magnetic drifts and the resulting velocity space nonuniformities and boundary losses. As an example, we have demonstrated in Section 8 how the orbit database can be used to model a distribution of alpha particles with a core-localized birth profile that is not a function of the fast alpha’s CoM. This method to design distribution functions was inspired by orbit-following Monte-Carlo simulations, but is much cheaper computationally. We expect that similar strategies can be used to model distribution functions based on various types of experimentally measured data or from results of bounce-averaged Fokker-Planck codes that ignore magnetic drifts Bierwage17c .

Finally, we demonstrated in Section 8.4 how a distribution of energetic alpha particles with large magnetic drifts can be converted back and forth between a 4-D mesh in arbitrary coordinates and the orbit-based representation in 3-D CoM space. Besides being an essential part of the coordinate transformation toolbox needed by platforms like IMAS, this procedure served us as a method to verify the numerical accuracy and logical integrity of the workflow. In fact, this procedure enabled us to identify and correct subtle problems that arose in early versions of our implementation of the algorithm. The lessons we learned are documented in this paper.

The methods presented here can serve as a starting point that already provides much of the required functionality required by IMAS. One of the next steps should be the adoption or development of smoothing techniques with the goal to mitigate the effect of discretization noise, while minimizing systematic errors associated with such irreversible data manipulations. This is a prerequisite for performing stability analyses of resonant modes and for preparing initial conditions for delta- $f$ simulations, which require explicit information about gradients in the phase space density.

We suggest the reader to also follow related and complementary developments in orbit tomography Stagner22 ; StagnerPhdThesis ; Jaerleblad21 ; BenjaminBachelorThesis . For instance, those studies demonstrate an alternative way to sample the orbit space and employ methods for adaptive meshing and reducing binning errors.

Acknowledgments

One of the authors (A.B.) acknowledges insightful discussions with Yasushi Todo concerning the handling of the Jacobian factor $B^{*}_{\parallel}$ , and with Stuart Benjamin and David Pfefferle concerning methods for phase space discretization and mapping techniques. On a broader note, A.B. would like to thank Liu Chen, Zhihong Lin, Fulvio Zonca as well as Sergio Briguglio, Claudio DiTroia and Gregorio Vlad for their support during the early development stages of VisualStart in 2008–2010 for the hybrid code HMGC. Subsequent extensions enabling VisualStart to initialize long-time simulations performed for JT-60U with the code MEGA benefited from support by Kouji Shinohara, Yasushi Todo and Masatoshi Yagi, to whom A.B. would also like to express his gratitude. Finally, A.B. is grateful to Jeronimo Garcia and Shunsuke Ide for encouraging and supporting a fruitful collaboration between JET and QST. The working example used in the present paper was inspired by plasmas produced in recent JET experiments with central 3-ion RF heating, which posed a suitably challenging and practically relevant scenario for testing our numerical algorithms and workflow. We would like to thank all JET contributors for making these interesting experiments possible and allowing us to participate.

The workstation used for the numerical calculations reported here was funded by QST President’s Strategic Grant (Creative Research). This work has been partially carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A Binning operations and Jacobians

In this section, we describe how we evaluated the Klimontovich representation (25) to visualize our results in Section 8. For clarity, we break the procedure up into (i) the computation of a histogram $h$ and (ii) its conversion to a density function $f$ . As a concrete example, we use the set of noncanonical coordinates $(E,\lambda,R,z)$ . The cell center positions are ${\bm{Z}}_{ijkm}=(E_{i},\lambda_{j},R_{k},z_{m})$ and the volume elements are $\Delta V_{ijkm}=\Delta E_{i}\Delta\lambda_{j}\Delta R_{k}\Delta z_{m}$ .¹¹¹¹11In contrast to the canonical 6-D phase space volume elements $\Delta\mathcal{V}_{\rm orb}$ and $\Delta\mathcal{V}$ of GC orbits and their markers as defined in Section 6, the symbol $\Delta V$ appearing here denotes an element in an arbitrary mesh in an arbitrary number of dimensions. The corresponding histogram of GCs is given by

h({\bm{Z}}_{ijkm})=\sum_{l=1}^{N_{\rm mrk}}w_{l}\frac{S({\bm{Z}}_{{\rm gc},l}-{\bm{Z}}_{ijkm})}{\Delta V_{ijkm}}.

(51)

Marker $l$ is located at GC position ${\bm{Z}}_{{\rm gc},l}$ , and its weight $w_{l}$ determines how many physical particles this marker is meant to represent (cf. Eq. (26)). The shape function $0\leq S\leq 1$ prescribes how we map these weights onto our mesh. The phase space density function $f$ is then

f=\frac{h}{\mathcal{J}^{({\bm{x}},{\bm{v}})}_{\rm gc}}\approx\frac{h({\bm{Z}}_{ijkm})}{\left[\frac{B^{*}_{\parallel}}{B}\mathcal{J}^{\bm{x}}_{Rz}\mathcal{J}^{\bm{v}}_{E\lambda}\right]_{ijkm}};

(52)

where $\mathcal{J}^{({\bm{x}},{\bm{v}})}_{\rm gc}$ is the Jacobian for the conversion from Cartesian to noncanonical GC coordinates. Substituting Eq. (51) into (52) and assuming, within the accuracy limits of our discrete mesh, that $[...]_{ijkm}\approx[...]_{l}$ , the formula for mapping the phase space density function $f$ to our 4-D mesh becomes¹²¹²12Conventional PIC codes sample the phase space by loading marker particles in a uniformly randomized manner at the beginning of a simulation ( $t=0$ ). In that case, the initial marker weights (via their volume elements) carry a factor $B^{*}_{\parallel}(t=0)/B(t=0)$ . When evaluating moments of the distribution function during a simulation, that factor is effectively multiplied by the inverse factor $B(t)/B^{*}_{\parallel}(t)$ that appears in Eq. (53). These two factors cancel approximately and may, thus, be omitted, if one ignores corrections of order $\mathcal{O}(\rho_{0}^{2}/a^{2})$ .

f(E_{i},\lambda_{j},R_{k},z_{m})=\sum_{l=1}^{N_{\rm mrk}}\frac{w_{l}\times S({\bm{Z}}_{{\rm gc},l}-{\bm{Z}}_{ijkm})}{\left[\frac{B^{*}_{\parallel}}{B}\mathcal{J}^{\bm{x}}_{Rz}\mathcal{J}^{\bm{v}}_{E\lambda}\right]_{l}\Delta V_{ijkm}}.

(53)

Based on Eq. (53), we compute density fields in 1-D and 2-D by binning marker weights $w_{l}$ as follows:

	$\displaystyle n(\psi_{{\rm P},i})=$	$\displaystyle\sum_{l=1}^{N_{\rm mrk}}\frac{w_{l}}{\left[\frac{B^{*}_{\parallel}}{B}\mathcal{J}^{\bm{x}}_{\psi_{\rm P}}(R,z)\right]_{l}}\frac{S({\bm{Z}}_{{\rm gc},l}-{\bm{Z}}_{i})}{\Delta\psi_{{\rm P},i}},$		(54)
	$\displaystyle n(R_{i},z_{j})=$	$\displaystyle\sum_{l=1}^{N_{\rm mrk}}\frac{w_{l}}{\left[\frac{B^{*}_{\parallel}}{B}\mathcal{J}^{\bm{x}}_{Rz}(R)\right]_{l}}\frac{S({\bm{Z}}_{{\rm gc},l}-{\bm{Z}}_{ij})}{\Delta R_{i}\Delta z_{j}}.$		(55)

The spatially integrated velocity distributions are

\displaystyle\nu(a_{i},b_{j})=

\displaystyle\sum_{l=1}^{N_{\rm mrk}}\frac{w_{l}}{\left[\frac{B^{*}_{\parallel}}{B}\hat{\mathcal{J}}_{ab}^{\bm{v}}\right]_{l}}\frac{S({\bm{Z}}_{{\rm gc},l}-{\bm{Z}}_{ij})}{\Delta a_{i}\Delta b_{j}}.

(56)

In Sections 8.2 and 8.3, we used 1st-order shape functions $S$ for binning (linear interpolation), and in Section 8.4 we used 0th-order (binary) binning.

Finally, we present expressions for the Jacobian factors multiplying $B^{*}_{\parallel}/B$ . Converting the right-hand side of Eq. (6) to polar GC velocity coordinates by substituting ${\rm d}\mu\rightarrow B^{-1}M{\rm d}v_{\perp}^{2}/2$ , we obtain

{\rm d}\hat{v}_{x}{\rm d}\hat{v}_{y}{\rm d}\hat{v}_{z}\times{\rm d}^{3}{\bm{x}}=\frac{B^{*}_{\parallel}}{B}\times 2\pi\hat{v}_{\perp}{\rm d}\hat{v}_{\perp}{\rm d}\hat{v}_{\parallel}\times{\rm d}^{3}{\bm{x}}.

(57)

The Jacobian $\mathcal{J}^{\bm{x}}_{\psi_{\rm P}}$ for transforming between ${\rm d}^{3}{\bm{x}}$ and the normalized poloidal flux $\psi_{\rm P}$ is evaluated numerically. The Jacobian for cylinder coordinates is $\mathcal{J}_{Rz}^{\bm{x}}=2\pi R$ . Analytical expressions can also be derived for the Jacobians of all velocity coordinates that we use here. Namely, from the relations

	$\displaystyle\hat{v}_{\perp}^{2}=$	$\displaystyle 2\hat{E}\Lambda\hat{B}=2\hat{E}\cos^{2}\alpha=2\hat{E}(1-\lambda^{2}),$		(58)
	$\displaystyle\hat{v}_{\parallel}=$	$\displaystyle\sqrt{2\hat{E}}\sin\alpha=\sqrt{2\hat{E}}\lambda=\sigma_{B}\sqrt{2\hat{E}(1-\Lambda\hat{B})},$		(59)

with $\sigma_{B}\equiv v_{\parallel}/|v_{\parallel}|$ , one readily obtains the Jacobians

\mathcal{J}_{ab}^{\bm{v}}=2\pi\left|\begin{array}[]{cc}\frac{1}{2}\partial_{a}v_{\perp}^{2}&\partial_{a}v_{\parallel}\\ \frac{1}{2}\partial_{b}v_{\perp}^{2}&\partial_{b}v_{\parallel}\end{array}\right|

(60)

for locally (in ${\bm{x}}$ ) transforming velocity space elements $2\pi v_{\perp}{\rm d}v_{\perp}{\rm d}v_{\parallel}$ to other sets $(a,b)$ , such as

	$\displaystyle\hat{\mathcal{J}}_{\mathcal{E}\mu}^{\bm{v}}=\pi\hat{B}/\|\hat{v}_{\parallel}\|,\quad\hat{\mathcal{J}}_{\mathcal{E}\Lambda}^{\bm{v}}=2\pi{\hat{E}\hat{B}}/\|\hat{v}_{\parallel}\|,$		(61)
	$\displaystyle\hat{\mathcal{J}}_{E\alpha}^{\bm{v}}=\hat{\mathcal{J}}_{v_{\perp},v_{\parallel}}=2\pi\hat{v}_{\perp},\quad\hat{\mathcal{J}}_{E\lambda}^{\bm{v}}=2\pi\hat{v}.$		(62)

Appendix B Numerical accuracy and data storage

Researchers tend to invest significant effort in the preparation of their simulation scenarios, balancing accuracy against speed. However, the data on platforms like IMAS are likely to be of variable quality, since some of it originates from automated or semi-automated workflows. The tools used to process those data should hence have a certain degree of tolerance with respect to inaccuracies and be able to detect and handle exceptions. The workflow in VisualStart also contains such sensibility checks and ways to mitigate the effects of inaccuracies. Some of the inaccuracies arising in VisualStart could be avoided with more sophisticated techniques, while others are practically inevitable due to numerical resolution being necessarily finite. This Appendix contains notes concerning some issues that are relevant for the subject of this paper or, at least, for the current implementation of VisualStart.

B.1 Particle pushing

At present, the particle pushing routine used to compute GC orbits in VisualStart uses the same finite-difference algorithm as the MHD-PIC hybrid code MEGA Todo98 ; Todo05 : a 4th-order Runge-Kutta scheme advances the equations of motion expressed in cylinder coordinates. Spatial derivatives are also taken through simple finite differences. This choice has been made deliberately in order to have the same level of accuracy as the hybrid code for which VisualStart prepares initial conditions. This allows us to detect possible problems, if any, at an early stage, before launching expensive hybrid simulations. Of course, accuracy can be improved with various techniques, including a Hamiltonian formulation like that used in codes like ORBIT White84 ; ChenY99 ; White19 and HAGIS Pinches98 ; Pinches04 , albeit with the complication of having to work with a specialized set of coordinates. The benefits of such techniques become noticeable when following orbits for long periods of time (as in Poincaré analyses). For evaluating a single poloidal transit as we do here, simple techniques often suffice, but due attention must always be paid to regions near singularities and stagnation points like those shown in Figs. 7 and 8.

B.2 Orbit database

Although one may load GC phase space markers as soon as an orbit contour has been computed, it is advantageous to store orbits in a database and load the markers later as needed. One reason is that the orbits can be recycled when preparing distribution functions for different sets of particles: alpha particles and deuterons are indistinguishable for identical velocities, and one may have multiple beams as well as RF-heated populations of a certain particle species in the same plasma. These sets of particles may all share the same orbit database.

In order to keep the size of the database manageable, the representation of the orbits must be efficient. VisualStart saves all orbits in discretized form, typically using 32–64 samples per orbit. The user can choose between three options for the distribution of samples: (i) uniform in space, (ii) uniform in time, or (iii) a 50:50 combination of (i) and (ii). We usually choose option (i) for beam-like distributions that consist only of passing particles. Option (iii) is preferred when trapped particles are present in order to ensure good resolution both near and far away from V-type stagnation points as the examples in Fig. 19 show.

The recorded time arrays must include at least $R(t)$ , $Z(t)$ , $v_{\parallel}(t)$ and the time $t$ itself. Our databases contain many other arrays that have been stored for various purposes, such as $\zeta(t)$ , $v_{R}(t)$ , $v_{z}(t)$ , $v_{\zeta}(t)$ and $\psi_{\rm P}(t)$ . The database file from our low-resolution example in Fig. 6 has a size of 33 Mbyte in NetCDF format. The databases with the parameters in Table 2 occupy about 150 Mbyte each in the mono-energetic cases and 500 Mbyte in the flat-energy cases. These files also contain useful auxiliary information, such as the numerical parameters, poloidal and toroidal transit times, the orbit length, bounce angle, the orbit type, and some measures of the computational accuracy: energy conservation, and spatial mismatch between start- and end-point.

Orbits that have not been completed within a certain upper bound of time steps are labeled as “incomplete”. Before proceeding to Step 4 of our workflow in Fig. 2, we deal with such incomplete orbits and those that have insufficient accuracy. Bad orbits can be corrected if necessary, or simply discarded if sufficient accuracy is difficult to achieve and unnecessary (as may happen for tiny orbits very close to an O-type stagnation point).

Data related to marker weighting are also present; in particular, the volume element sizes $\Delta\mathcal{V}_{\rm orb}$ . There is an option to undo the factor $1/2$ in the volume element given by Eq. (37) for stagnation orbits that are smaller than the cell size and, hence, cannot be double-counted. Their volume elements may then be rescaled in proportion to the orbit’s radial diameter, but the effect is usually negligible since such orbits carry little weight (and one should use a finer mesh if they are of interest).

Interpolation is necessary when converting the raw data from the particle pushing algorithm to the samples used to represent an orbit in the database. Further interpolation is necessary when loading markers on an orbit that has been retrieved from the database. As shown in Fig. 20, we found that a tiny discontinuity in an orbit contour (mismatch between start and end points) can cause relatively large “ringing” in spline interpolations. Subsequent interpolations in that region — namely, when we load marker particles — can produce corrupt samples that may even lie outside the plasma. When applying the 1-D interpolation function ${\tt interp1}$ of Matlab for processing our orbits (and other potentially “noisy” data), we avoid such problems by using the pchip option instead of spline.

B.3 Marker loading

In order to avoid crowding at the midplane, the markers on each orbit are loaded with a random time offset $0\leq\tau_{0}<\tau_{\rm pol}$ relative to the starting point at the midplane. When using the markers in a simulation, they are also spread (randomly) along the toroidal angle $\zeta$ .

B.4 Convergence test

Figure 17(b) showed indications of relatively large errors that increased systematically with each iteration of our back-and-fourth transformation between a 4-D mesh in the coordinates $(E,\lambda,R,z)$ and the orbit-based representation in 3-D CoM space. In our discussion in Section 8.4, we attributed this problem to inaccuracies in the estimation of the volume elements for orbits sampling the region near the trapped-passing boundary, where some factors exhibit singular behavior.

Further tests show that these errors can be reduced by simply increasing the number of phase space samples via the parameters $N_{\alpha}$ , $N_{X}$ and $N_{\tau}$ in Table 2.¹³¹³13Since kinetic energy is conserved and our 3-D and 4-D distributions both share the same energy mesh, the parameter $N_{E}$ has no impact on the numerical accuracy in our verification exercise. In Fig. 21 we demonstrate this for the case where we increased the number of radial grid points $N_{X}$ . While the flat-energy case with $E=0...3.5\,{\rm MeV}$ was used in Fig. 17, we have chosen only mono-energetic alpha particles with $E=3.5\,{\rm MeV}$ for the convergence test in Fig. 21. For an identical mesh, this results in larger values of the relative errors in Fig. 17, presumably because the signal-to-noise ratio is reduced.

The comparison between panels (a) and (b) of Fig. 21 shows a significant reduction of the numerical errors in both the local values of the alpha particle density profile $\hat{n}^{(\eta)}(r)$ and the total number of GCs $\mathcal{N}$ . As before, the errors accumulate with each iteration, which confirms their systematic character. Large errors persist only at the artificial loss boundary ( $r/a\approx 1$ ), which is to be expected and requires dedicated solutions that are not covered here.

Besides increasing the number of samples, there may exist other ways to correct or otherwise suppress the inaccuracies associated with singular volume elements. This is left for future work. A preliminary discussion was given in Section 6.2 of the main text.

References

[1] S.D. Pinches et al. Progress in the ITER integrated modelling programme and the ITER scenario databases. 27th IAEA Fusion Energy Conference (FEC 2018), 2018. https://nucleus.iaea.org/sites/fusionportal/ Shared%20Documents/FEC%202018/fec2018-preprints/ preprint0741.pdf.
[2] M. Salewski, M. Nocente, G. Gorini, A.S. Jacobsen, V.G. Kiptily, S.B. Korsholm, F. Leipold, J. Madsen, D. Moseev, S.K. Nielsen, J. Rasmussen, M. Stejner, M. Tardocchi, and JET Contributors. Velocity-space observation regions of high-resolution two-step reaction gamma-ray spectroscopy. Nucl. Fusion, 55(9):093029, 2015.
[3] M. Salewski, M. Nocente, G. Gorini, A.S. Jacobsen, V.G. Kiptily, S.B. Korsholm, F. Leipold, J. Madsen, D. Moseev, S.K. Nielsen, J. Rasmussen, M. Stejner, M. Tardocchi, and JET Contributors. Fast-ion energy resolution by one-step reaction gamma-ray spectrometry. Nucl. Fusion, 56(4):040009, 2016.
[4] M. Salewski, M. Nocente, A.S. Jacobsen, F. Binda, C. Cazzaniga, J. Eriksson, B. Geiger, G. Gorini, C. Hellesen, V.G. Kiptily, T. Koskela, S.B. Korsholm, T. Kurki-Suonio, F. Leipold, D. Moseev, S.K. Nielsen, J. Rasmussen, P.A. Schneider, S.E. Sharapov, M. Stejner, M. Tardocchi, and JET Contributors, ASDEX Upgrade Team & EUROfusion MST1 Team. Bayesian integrated data analysis of fast-ion measurements by velocity-space tomography. Fusion Sci. Technol., 74(1-2):23, 2018.
[5] M. Salewski, M. Nocente, B. Madsen, I. Abramovic, G. Gorini, A.S. Jacobsen, V.G. Kiptily, S.B. Korsholm, D. Moseev, S.K. Nielsen, A.F.L. Poulsen, J. Rasmussen, M. Tardocchi, B. Geiger, J. Eriksson, the JET Contributors, the ASDEX Upgrade Team, and the EUROfusion MST1 Team. Diagnostic of fast-ion energy spectra and densities in magnetized plasmas. JINST, 14:C05019, 2019.
[6] L. Stagner, W.W. Heidbrink, M. Salewski, A.S. Jacobsen, B. Geiger, and The DIII-D and ASDEX Upgrade Teams. Orbit tomography of energetic particle distribution functions. Nucl. Fusion, 62(2):026033, 2022.
[7] L. Stagner. Inference of the fast-ion distribution function. PhD thesis, University of California, Irvine, 2018.
[8] H. Järleblad, L. Stagner, M. Salewski, J. Eriksson, S. Benjamin, B. Madsen, M. Nocente, J. Rasmussen, and B.S. Schmidt. Fast-ion orbit sensitivity of neutron emission spectroscopy diagnostics. Rev. Sci. Instr., 92(4):043526, 2021.
[9] J. Breslau and D. Liu. Reconstructing fast ion distribution functions from NUBEAM. TRANSP Users Group Meeting, May 4, 2017.
[10] Ph. Lauber et al. Energetic particle dynamics induced by off-axis neutral beam injection on ASDEX Upgrade, JT-60SA and ITER. 28th IAEA Fusion Energy Conference (FEC 2020), May 10–15, 2021.
[11] E. Tholerus and M. Fitzgerald. Calculating the fast ion distribution function in constants-of-motion space for MHD instability analysis. Joint T17-07 fast ion physics group & M15-24 meeting, May 31, 2017.
[12] A. Bierwage. VisualStart: GUI-Aided Unified Initialization Tool for Hybrid (MHD + Particle) Simulations. 17th NEXT Meeting, Kashiwa, Japan, March 15-–16, 2012.
[13] A. Bierwage, C. Di Troia, S. Briguglio, and G. Vlad. Orbit-based representation of equilibrium distribution functions for low-noise initialization of kinetic simulations of toroidal plasmas. Comp. Phys. Comm., 183:1107, 2012.
[14] A. Bierwage, K. Shinohara, N. Aiba, and Y. Todo. Role of convective amplification of $n=1$ energetic particle modes for N-NB ion dynamics in JT-60U. Nucl. Fusion, 53(7):073007, 2013.
[15] A. Bierwage, Y. Todo, N. Aiba, and K. Shinohara. Dynamics of low- $n$ shear Alfvén modes driven by energetic N-NB ions in JT-60U. Nucl. Fusion, 54(10):104001, 2014.
[16] A. Bierwage, K. Shinohara, Y. Todo, N. Aiba, M. Ishikawa, G. Matsunaga, M. Takechi, and M. Yagi. Self-consistent long-time simulation of chirping and beating energetic particle modes in JT-60U plasmas. Nucl. Fusion, 57(1):016036, 2017.
[17] A. Bierwage, K. Shinohara, Y. Todo, N. Aiba, M. Ishikawa, G. Matsunaga, M. Takechi, and M. Yagi. Simulations tackle abrupt massive migrations of energetic beam ions in a tokamak plasma. Nature Comms., 9:3282, 2018.
[18] A. Bierwage and K. Shinohara. Orbit-based analysis of resonant excitations of Alfvén waves in tokamaks. Phys. Plasmas, 21(11):112116, 2014.
[19] A. Bierwage and K. Shinohara. Orbit-based analysis of nonlinear energetic ion dynamics in tokamaks. I. Effective mode number profile and resonant frequency tracking. Phys. Plasmas, 23(4):042511, 2016.
[20] A. Bierwage and K. Shinohara. Orbit-based analysis of nonlinear energetic ion dynamics in tokamaks. II. Mechanisms for rapid chirping and convective amplification. Phys. Plasmas, 23(4):042512, 2016.
[21] F. Zonca, L. Chen, S. Briguglio, G. Fogaccia, G. Vlad, and X. Wang. Nonlinear dynamics of phase space zonal structures and energetic particle physics in fusion plasmas. New J. Phys., 17(1):013052, 2015.
[22] L. Chen and F. Zonca. Physics of Alfvén waves and energetic particles in burning plasmas. Rev. Mod. Phys., 88:15008, 2016.
[23] M.V. Falessi and F. Zonca. Transport theory of phase space zonal structures. Phys. Plasmas, 26(2):022305, 2019.
[24] M. V. Falessi F. Zonca, L. Chen and Z. Qiu. Nonlinear radial envelope evolution equations and energetic particle transport in tokamak plasmas. J. Phys.: Conf. Ser., 1785:012005, 2021. Varenna-Lausanne 2020, 12-16 October, Lausanne, Switzerland.
[25] S. Benjamin. Distribution transforms for orbit tomography. Bachelor Thesis, The Australian National University, 2021. Paper in preparation.
[26] J.A. Rome and Y-K.M. Peng. The topology of tokamak orbits. Nucl. Fusion, 19(9):1193, 1979.
[27] Ph. Lauber and Z. Lu. Analytical finite-lamor-radius and finite-orbit-width model for the LIGKA code and its application to KGAM and shear Alfvén physics. J. Phys.: Conf. Ser., 1125:012015, 2018. Varenna-Lausanne 2018, 27-31 August, Varenna, Italy.
[28] A. Bierwage, M. Toma, and K. Shinohara. MHD and resonant instabilities in JT-60SA during current ramp-up with off-axis N-NB injection. Plasma Phys. Control. Fusion, 59(12):125008, 2017.
[29] R.G. Littlejohn. Variational principles of guiding centre motion. J. Plasma Physics, 29:111, 1983.
[30] D. Moseev and M. Salewski. Bi-Maxwellian, slowing-down, and ring velocity distributions of fast ions in magnetized plasmas. Phys. Plasmas, 26(2):020901, 2019.
[31] M. Nocente et al. Generation and observation of fast deuterium ions and fusion-born alpha particles in JET D–³He plasmas with the 3-ion radio-frequency heating scenario. Nucl. Fusion, 60(12):124006, 2020.
[32] Y.O. Kazakov et al. Physics and applications of three-ion ICRF scenarios for fusion research. Phys. Plasmas, 28(2):020501, 2021.
[33] A. Bierwage, K. Shinohara, Ye.O. Kazakov, V. Kiptily, Ph. Lauber, M. Nocente, Ž. Štancar, S. Sumida, M. Yagi, J. Garcia, S. Ide, and JET Contributors. Energy-selective confinement of fusion-born alpha particles during internal relaxations in a tokamak plasma. Submitted, 2022. Preprint: http://arxiv.org/abs/2109.03427.
[34] M. Weiland, R. Bilato, R. Dux, B. Geiger, A. Lebschy, F. Felici, R. Fischer, D. Rittich, M. van Zeeland, the ASDEX Upgrade Team, and the Eurofusion MST1 Team. RABBIT: Real-time simulation of the NBI fast-ion distribution. Nucl. Fusion, 58(8):082032, 2018.
[35] Ž. Štancar, Z. Ghani, J. Eriksson, A. Žohar, S. Conroy, Ye. O. Kazakov, T. Craciunescu, M. Nocente, L. Garzotti, V. Radulović, P. Sirén, V. Kiptily, K. Kirov, Y. Baranov, G. Szepesi, M. Dreval, M. Gorelenkova, H. Weisen, E. Militello-Asp, L. Snoj, and JET Contributors. Experimental validation of an integrated modelling approach to neutron emission studies at JET. Nucl. Fusion, 61:126030, 2021.
[36] F. Porcelli, R. Stankiewicz, W. Kerner, and H.L. Berk. Solution of the drift-kinetic equation for global plasma modes and finite particle orbit widths. Phys. Plasmas, 1(3):470, 1994.
[37] Y. Zhao and R.B. White. Simulation of $\alpha$ -particle redistribution due to sawteeth on the Tokamak Fusion Test Reactor. Phys. Plasmas, 4(4):1103, 1997.
[38] K. Tani, M. Azumi, H. Kishimoto, and S. Tamura. Effect of toroidal field ripple on fast ion behavior in a tokamak. J. Phys. Soc. Japan, 50(5):1726, 1981.
[39] K. Tani and M. Azumi. Simulation studies on alpha-particle-driven current in tokamak reactors. Nucl. Fusion, 48(8):085001, 2008.
[40] Y. Todo and T. Sato. Linear and nonlinear particle-magnetohydrodynamic simulations of the toroidal Alfvén eigenmode. Phys. Plasmas, 5(5):1321, 1998.
[41] Y. Todo, K. Shinohara, M. Takechi, and M. Ishikawa. Nonlocal energetic particle mode in a JT-60U plasma. Phys. Plasmas, 12(1):012503, 2005.
[42] Y. Todo, M.A. Van Zeeland, A. Bierwage, and W.W. Heidbrink. Multi-phase simulation of fast ion profile flattening due to Alfvén eigenmodes in a DIII-D experiment. Nucl. Fusion, 54:104012, 2014.
[43] J. R. Cary and A. J. Brizard. Hamiltonian theory of guiding-center motion. Rev. Mod. Phys., 81:693, 2009.
[44] M. Deutsch. Angular correlations in nuclear reactions. Rep. Progr. Phys., 14(1):196, 1951.
[45] A.O. Hanson, R.F. Taschek, and J. H. Williams. Monoergic neutrons from charged particle reactions. Rev. Mod. Phys., 21(4):635, 1949.
[46] M. Salewski, M. Nocente, B. Madsen, I. Abramovic, M. Fitzgerald, G. Gorini, P.C. Hansen, W.W. Heidbrink, A.S. Jacobsen, T. Jensen, V.G. Kiptily, E.B. Klinkby, S.B. Korsholm, T. Kurki-Suonio, A.W. Larsen, F. Leipold, D. Moseev, S.K. Nielsen, S.D. Pinches, J. Rasmussen, M. Rebai, M. Schneider, A. Shevelev, S. Sipilä, M. Stejner, and M. Tardocchi. Alpha-particle velocity-space diagnostic in ITER. Nucl. Fusion, 58(9):096019, 2018.
[47] S. Sumida, K. Shinohara, R. Ikezoe, M. Ichimura, M. Sakamoto, M. Hirata, and S. Ide. Characteristics of fast ³He ion velocity distribution exciting ion cyclotron emission on JT-60U. Plasma Phys. Control. Fusion, 61(2):025014, 2019.
[48] S. Sumida, K. Shinohara, M. Ichimura, T. Bando, A. Bierwage, and S. Ide. Identification of slow-wave ion cyclotron emission on JT-60U. Nucl. Fusion, 61(11):116036, 2021.
[49] A. Bierwage, Y. Todo, N. Aiba, and K. Shinohara. Sensitivity study for N-NB-driven modes in JT-60U: Boundary, diffusion, gyroaverage, compressibility. Nucl. Fusion, 56(10):106009, 2016.
[50] ITER physics data model documentation (restricted access). https://sharepoint.iter.org/departments/ POP/CM/IMDesign/Data%20Model/CI/imas-3.33.0/ html_documentation.html. Status: Aug. 2021.
[51] R.B. White and M.S. Chance. Hamiltonian guiding center drift orbit calculation for plasmas of arbitrary cross section. Phys. Fluids, 27(10):2455, 1984.
[52] Y. Chen, R.B. White, G.-Y. Fu, and R. Nazikian. Numerical study of the nonlinear evolution of toroidicity-induced Alfvén eigenmodes. Phys. Plasmas, 6:226, 1999.
[53] R.B. White, V.N. Duarte, N.N. Gorelenkov, E.D. Fredrickson, M. Podesta, and H.L. Berk. Modeling of chirping toroidal Alfvén eigenmodes in NSTX. Phys. Plasmas, 26:092103, 2019.
[54] S.D. Pinches, L.C. Appel, J. Candy, S.E. Sharapov, H.L. Berk, D. Borba, B.N. Breizman, T.C. Hender, K.I. Hopcraft, G.T.A. Huysmans, and W. Kerner. The HAGIS self-consistent nonlinear wave-particle interaction model. Comp. Phys. Comm., 111:133, 1998.
[55] S.D. Pinches, H.L. Berk, M.P. Gryaznevich, S.E. Sharapov, and JET-EFDA Contributors. Spectroscopic determination of the internal amplitude of frequency sweeping TAE. Plasma Phys. Control. Fusion, 46(7):S47, 2004.