Good news: The amount of spurious heating caused by N-point gyroaveraging in our guiding center simulations is significantly smaller than what was shown in Figs. 14(i)–14(x) of the original published paper.1 

The additional spurious heating was due to an error that we discovered in our simulation code; specifically, in the equation of motion for the toroidal angle φgc of the guiding center,
(1)
where êφ is the unit vector along the toroidal angle φ in right-handed cylinder coordinates (R,φ,z),vgc is the instantaneous guiding center velocity, and R is the major radius. The Jacobian factor 1/R on the right-hand side of Eq. (1) should be evaluated at the location of the guiding center, Rgc(t). However, when gyroaveraging was enabled, our code had accidentally assigned a corrupted value of the radial position Rsat(t) of the last satellite particle on the gyroaveraging circle. The impact this error had on spurious heating seems to have varied from case to case, and this mistake was also responsible for the mysterious inward or outward shifts of the resonance that were visible in Figs. 14(i)–14(x) of Ref. 1.

We have corrected the error and the new results are shown in Figs. 1(i)–1(x) of this Erratum. One can see from the good conservation of the rotating frame energy E(t) that spurious heating caused by gyroaveraging is now ignorable in our working example on the timescale of a few milliseconds.

FIG. 1.

(Replaces Fig. 14 of the original published paper.1) Comparison of results from simulations using the full orbit model (a)–(d), the GC model (e)–(h), and the GC model with N-point gyroaveraging (i)–(x). We followed co-passing 400keV deuterons, and the nonnormal mode of case (c) was used as a perturbation. Results are presented as Poincaré plots (columns 2 and 4) and time traces of the rotating frame energy Ê(t) (columns 3 and 5) for scenarios (i) and (ii) with electrostatic and electromagnetic perturbation, respectively. The first column shows schematically the differences between the three methods for representing gyration in the (R, z) plane. In the full orbit case, particles travel along a helix, whose poloidal projection is drawn here as a red circle of radius ϱg0.1m around the guiding center (black dot). In the GC model, fields are evaluated at the GC only. Gyroaveraging is done by placing Navg satellite particles (small red circles) on the gyrocircle and taking the average of the fluctuating fields E and δB at those locations to push the GC. Here we compare results obtained with Navg=2,4,8 and different satellite positions. As a visual aid, the Poincaré plots are overlaid with a pair of vertical dashed lines, which indicate the limits of the resonance's separatrix in the full orbit simulations.

FIG. 1.

(Replaces Fig. 14 of the original published paper.1) Comparison of results from simulations using the full orbit model (a)–(d), the GC model (e)–(h), and the GC model with N-point gyroaveraging (i)–(x). We followed co-passing 400keV deuterons, and the nonnormal mode of case (c) was used as a perturbation. Results are presented as Poincaré plots (columns 2 and 4) and time traces of the rotating frame energy Ê(t) (columns 3 and 5) for scenarios (i) and (ii) with electrostatic and electromagnetic perturbation, respectively. The first column shows schematically the differences between the three methods for representing gyration in the (R, z) plane. In the full orbit case, particles travel along a helix, whose poloidal projection is drawn here as a red circle of radius ϱg0.1m around the guiding center (black dot). In the GC model, fields are evaluated at the GC only. Gyroaveraging is done by placing Navg satellite particles (small red circles) on the gyrocircle and taking the average of the fluctuating fields E and δB at those locations to push the GC. Here we compare results obtained with Navg=2,4,8 and different satellite positions. As a visual aid, the Poincaré plots are overlaid with a pair of vertical dashed lines, which indicate the limits of the resonance's separatrix in the full orbit simulations.

Close modal
Notes of caution. This result is, in our opinion, non-trivial and therefore interesting, because the gyroaveraging operation g using a time-dependent gyroradius,
(2)
—with particle mass M, electric charge Ze, field strength B=|B|, and total (reference + fluctuating) magnetic field vector B(t)=Bref+δB(t)—commutes neither with the partial time derivative t nor with our finite-difference-version fd of the operator in Faraday's law,
(3)
Thus, although the results in Fig. 1 are encouraging for guiding center simulations that make use of gyroaveraging, we nevertheless recommend to exercise caution. To illustrate why, let us discuss a counterexample.
Gyroaveraging in guiding center simulations is a physically motivated but still somewhat arbitrary smoothing operation on the fields. Different choices may have a different impact on a simulation's conservation properties and, without in-depth analysis, the outcome may even seem counter-intuitive at first glance. As an example, consider Fig. 2, which shows the results of a simulation where—instead of the conventional ϱg(B) in Eq. (2)—we used a constant gyroradius ϱ0=Mv0/(ZeB0), where B0 is the field strength at the magnetic axis and v0 is the characteristic speed of the deuterons in our simulation (here, Mv02/2=400keV). In this case, the gyroaveraging operation commutes with the time derivative in Eq. (3), so one may expect to be making only a small error of order
(4)
Indeed, panels (c) and (d) of Fig. 2 show that the spurious heating is reduced when the spatial resolution is increased. The results in Figs. 2(a) and 2(b) were obtained with the same resolution as those in Figs. 1(w) and 1(x), and one can see that using ϱ0 instead of ϱg(B) breaks Hamiltonian conservation laws in such a way that errors tend to accumulate near the resonance, yielding noticeable spurious heating within a few milliseconds in the present case. Clarifying which kinds of gyroaveraging operations might cause accumulation of errors and which do not and why is a subject that deserves further study.
FIG. 2.

Counterexample showing that some forms of gyroaveraging can cause an accumulation of spurious heating. Here, we used a constant gyroradius ϱ0=Mv0/(ZeB0) instead of ρg(B(t)) given by Eq. (2). The left column shows Poincaré plots as functions of normalized canonical toroidal angular momentum pφ, poloidal angle ϑ, and normalized rotating frame energy Ê (see the original paper1 for details and definitions). The right column shows time traces of Ê(t). The results in panels (a) and (b) were obtained with the same spatial resolution as in the published paper: NψP×Nϑ=1025×2048 for constructing the mode, and NR=Nz=480×480 in the guiding center simulation. The results in panels (c) and (d) were obtained with twice as many grid points: NψP×Nϑ=2049×4096, and NR=Nz=960×960. Along the toroidal direction, we used the particle-in-Fourier method. The number of satellites was Navg=8 as in Figs. 1(w) and 1(x). The results seem to be converged for Navg4. A comparison between panels (b) and (d) shows that higher resolution reduces the spurious heating, as expected from Eq. (4). The only exception seems to be the fifth sample from the top, which has been colored black for emphasis. This may have various reasons, one of which is that the form and location of the resonant island differs somewhat between panels (a) and (c), so that we are effectively sampling different portions of the resonance. This is most obvious if one looks at the fourth sample from the top (colored magenta), which lies outside the separatrix in (a) and inside the separatrix in (c).

FIG. 2.

Counterexample showing that some forms of gyroaveraging can cause an accumulation of spurious heating. Here, we used a constant gyroradius ϱ0=Mv0/(ZeB0) instead of ρg(B(t)) given by Eq. (2). The left column shows Poincaré plots as functions of normalized canonical toroidal angular momentum pφ, poloidal angle ϑ, and normalized rotating frame energy Ê (see the original paper1 for details and definitions). The right column shows time traces of Ê(t). The results in panels (a) and (b) were obtained with the same spatial resolution as in the published paper: NψP×Nϑ=1025×2048 for constructing the mode, and NR=Nz=480×480 in the guiding center simulation. The results in panels (c) and (d) were obtained with twice as many grid points: NψP×Nϑ=2049×4096, and NR=Nz=960×960. Along the toroidal direction, we used the particle-in-Fourier method. The number of satellites was Navg=8 as in Figs. 1(w) and 1(x). The results seem to be converged for Navg4. A comparison between panels (b) and (d) shows that higher resolution reduces the spurious heating, as expected from Eq. (4). The only exception seems to be the fifth sample from the top, which has been colored black for emphasis. This may have various reasons, one of which is that the form and location of the resonant island differs somewhat between panels (a) and (c), so that we are effectively sampling different portions of the resonance. This is most obvious if one looks at the fourth sample from the top (colored magenta), which lies outside the separatrix in (a) and inside the separatrix in (c).

Close modal

Additional corrections. The minus signs of E=Φ and E=Φ were missing in Fig. 16 of the published paper.1 Note that E=δE.

A.B. thanks Chang Liu (PPPL) for stimulating discussions that led to the discovery of the above-mentioned coding error.

1.
A.
Bierwage
and
K.
Shinohara
, “
Testing the conservative character of particle simulations: II. Spurious heating of guiding centers and full orbits subject to fluctuations expressed in terms of E and B
,”
Phys. Plasmas
29
,
113906
(
2022
).