7th Asia-Pacific Conference on Plasma Physics, 12-17 Nov, 2023 at Port Messe Nagoya

On gyroaveraging in drift-kinetic simulations of fast particle dynamics

Andreas Bierwage1, Kouji Shinohara2,1, Chang Liu3
1
National Institutes for Quantum Science and Technology (QST),
2
Dept. Complexity Science & Engineering, U. Tokyo, 3 Princeton Plasma Physics Laboratory (PPPL)
e-mail (speaker): bierwage.andreas@qst.go.jp
Gyroaveraging in numerical simulations is often necessaryfor obtaining a realistic spectrum of resonant instabilities
and realistic saturation amplitudes because the effectivestrength of wave-particle interactions in a magnetizedplasma depends on the ratio ρ / λ⊥ of gyroradius andperpendicular wavelength. A rigorous recipe is provided
by gyrokinetic theory, which splits the background field
Bref from the fluctuations δB. Only the latter are
gyroaveraged and the motion of charged particles is thenrepresented by the motion of gyrocenters. Here, we areinterested in formulations that make use of physicalelectric and magnetic fields, E and B (rather than theirpotentials Φ and A), which requires that Faraday's law∂B/∂t = – ∇ × E and the solenoidal condition ∇·B = 0 are
both enforced explicitly. In this case, rigorous
gyroaveraging in a gyrokinetic model requires that E andδB are subject to two-dimensional gyrodisk integrals(Stokes’ theorem), which is computationally expensive.
In practice, it can be meaningful to sacrifice some rigor
for convenience and computational speed, if the errors
made are tolerable. For instance, the gyrodisk integral ina gyrokinetic model formulated in terms of E and B may
be approximated by a reduced one-dimensional gyrocircle
integral, although this is justified only when thegyroradius is guaranteed to be small compared to theperpendicular wavelength for all simulation particles [1].
Here, we consider another choice that is often used in
kinetic orbit-following codes and MHD-PIC hybrid simulations of fast particles: instead of gyrokinetics, one solves
drift-kinetic equations of motion while performing gyrocircle-averaging in the poloidal (R, z) plane around guiding
centers. Hamiltonian guiding center theory does notincorporate gyroaveraging, and the impact of ad hoc gyroaveraging on conservation laws was recently analyzed [2].
The example studied was a JT-60U plasma containing
beam ions (deuterons) with a kinetic energy of Mv02/2 =400 keV and fairly large gyroradius: ρ0 = Mv0 / eB0 ~ 10 cm.
Conventional gyroaveraging using ρ(B) = Mv⊥(B) / eBref,with B = |B| yields fields <E> and <B> that violate Faraday's law: ∂<B>/∂t ≠ <∂B/∂t> = –<∇ × E> ≠ –∇ × <E>.
However, we found that the errors in the form of spurious
heating were tolerable or hardly noticeable at all as in theexample shown in Fig.1(c,d). For comparison, Fig.1 also
shows results for (a,b) Hamiltonian guiding center motion,and (e,f) conservative gyroaveraging [2], where <δB> isreplaced by a modified field satisfying ∂tδBmod = –∇ × <E>(a procedure that one may refer to as “Faraday cleaning”).
Gyroaveraging in a guiding center simulation can beviewed as a physically motivated smoothing operation on
the fields. Different choices will have different effects ona simulation’s conservation properties. Without in-depth
analysis, the outcome may even seem counter-intuitive atfirst glance as in Fig.1(g,h), which shows the results of asimulation where we used a constant gyroradius ρ0 insteadof ρ(B) with B taken at the current guiding center position.
In this case, one may expect to be making only a smallerror |<∇ × E> – ∇ × <E>| ~ |E(x+Δx) – E(x)| / ρ0, but
Fig.1(g,h) shows that this error tends to accumulate here
near the resonance, yielding noticeable spurious heating.
We also found that, in the presence of resonant fluctuations that resemble “normal” (eigen)modes of the system,
spurious heating tends to be slow and becomes significant
only on a time scale of 10 ms, where collisions andsources must also be considered. In contrast, spurious
heating was found to be much faster in the presence ofstrongly distorted “nonnormal” fluctuation patterns [2].
While this amplification effect remains to be understood,
it can be useful for code verification (as done in Fig.1).
Of course, in reality, such nonnormal patterns tend to be
short-lived and have rapidly varying phases, whichshould reduce their net spurious heating effect.

Figure 1. Conservation properties of a guiding center
simulation without (a,b) and with gyroaveraging (c-h).
The central column shows Poincaré plots of a portion ofa resonance. The right column shows time traces of therotating frame energy ℰ’ = E + ωPφ / n, which should beconserved in the present setup. See Ref. [2] for details.
References[1] Porazik & Lin, Commun. Comput. Phys. 10 (2011) 899.
[2] Bierwage & Shinohara, Phys. Plasmas 29 (2022) 113906.
Fig.1(c,d) above is a corrected version of Fig.14(w,x) of [2],where a coding error caused enhanced spurious heating.