NNonresonant Diffusion in Alpha Channeling
Ian E. Ochs and Nathaniel J. Fisch
Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08540, USA (Dated: December 14, 2020)The gradient of fusion-born alpha particles that arises in a fusion reactor can be exploited toamplify waves, which cool the alpha particles while diffusively extracting them from the reactor. Thecorresponding extraction of the resonant alpha particle charge has been suggested as a mechanismto drive rotation. By deriving a coupled linear-quasilinear theory of alpha channeling, we show that,for a time-growing wave with a purely poloidal wavevector, a current in the nonresonant ions cancelsthe resonant alpha particle current, preventing the rotation drive but fueling the fusion reaction.
Introduction:
A particle gyrating in a magnetic fieldwith a velocity v ⊥ greater than the phase velocity ω r /k ⊥ of an electrostatic wave will become Landau-resonant atsome point in its orbit, allowing for efficient wave-particleenergy exchange. Each time energy is exchanged with thewave, the particle gyrocenter position changes as well,leading to a wave-induced diffusion along a specified pathin energy-gyrocenter space [1, 2]. If the diffusion alongthis path on average pushes particles from higher to lowerenergy, the wave will amplify. This effect is known as al-pha channeling, so named because of its utility in simul-taneously cooling and extracting alpha particles from thehot core of a fusion reactor, and channeling the resultingenergy into useful wave power.For alpha channeling in a slab geometry, the diffusionpath slope in energy-gyrocenter space has the simple form ∂ X /∂K = k × ˆ b/m α ω Ω α , where K is the perpendicularkinetic energy, X is the gyrocenter position, m α and Ω α are the alpha particle mass and gyrofrequency, and ω and k are the wave frequency and wavenumber. Thus,the condition for alpha channeling is: ∂∂K + k × ˆ bm α ω Ω α · ∂∂ X ! F α > , (1)where F α is the zeroth-order distribution function inenergy-gyrocenter space.What remains unknown is whether or not the alphaparticles carry net charge out of the plasma as a resultof the wave-induced diffusion. If charge is in fact carriedout, then alpha channeling can be used drive E × B ro-tation in the plasma, providing an advantageous mecha-nism for shear rotation drive and centrifugal confinementin mirror fusion reactors [3]. Understanding whether suchschemes are possible at all requires evaluation of the effectof the wave on the nonresonant particles, which has neverbeen examined for alpha channeling. Such reactions inthe nonresonant particles are extremely important in en-forcing momentum and energy conservation [4, 5], mak-ing theories that ignore them liable to error.The reason the nonresonant response has proved elu-sive is that there is no existing linear theory of alphachanneling. Typically, a coupled linear wave and quasi-linear particle system is necessary to calculate the non-resonant particle response. The elusivity of the lineartheory is related to the fact that Landau damping can- not be derived from the magnetized dispersion relation,a conundrum sometimes termed the Bernstein-Landauparadox [6, 7]. Derivation of the diffusion thus requiresa nonlinear calculation, which allows for stochastic diffu-sion of the particle throughout phase space above a cer-tain wave amplitude at which Landau-resonant particlesdephase from the wave [8–12].In this Letter, we show that a linear-quasilinear systemcan be derived by assuming the wave-particle dephas-ing, which we accomplish by transforming the familiar unmagnetized linear-quasilinear kinetic theory to gyro-center coordinates. To show that the system describesalpha channeling, we show that it recovers both the am-plification condition (Eq. (1)) and the nonlinear diffusioncoefficient [9] for channeling by lower hybrid (LH) waves.Treating the channeling problem in this way positionsus to answer the question of whether alpha channeling ex-tracts charge from the fusion reactor. We find that for theinitial value problem in a slab geometry, where an elec-trostatic wave with purely poloidal wavenumber grows intime, the charge flux from the resonant alpha particlesis canceled by an equal and opposite charge flux in thenonresonant particles, so that no reactor charging occurs.We also determine which particles carry this nonresonantreturn current, in both single- and multi-ion-species plas-mas. For LH waves, the nonresonant return current iscarried exclusively by fuel ions, so that alpha channelinghas the added benefit of fueling the fusion reaction whileextracting alpha particles. Linear theory:
For any electrostatic wave, the dis-persion relation obtained by linearizing and Fouriertransforming Poisson’s equation can be expressed as:0 = 1 + X s D s ; D s ≡ − πq s k ˜ n s ˜ φ , (2)where q s and n s are the charge and density of species s , φ is the potential, and tildes denote Fourier transforms.In the standard limit | ω i /ω r | (cid:28)
1, this becomes:0 = 1 + X s D r,s (3)0 = X s (cid:18) iω i ∂D r,s ∂ω r + iD i,s (cid:19) , (4)where D r,s and D i,s are the real and imaginary parts ofthe dispersion function Eq. (2) evaluated at real ω, k . a r X i v : . [ phy s i c s . p l a s m - ph ] D ec To treat the alpha channeling initial value problem, wetake the wavenumber k k ˆ x , the magnetic field B k ˆ z , andthe gradient of the gyrocenter distribution function to bealong ˆ y . Thus, y corresponds to the “radial” direction,and x to the “poloidal” direction. For simplicity, we spe-cialize to a LH wave, assuming cold fluid populations ofmagnetized electrons e and unmagnetized ions i , and ahot, unmagnetized population of alpha particles α . Ourdispersion components are thus given by: D r,e = ω pe / Ω e , D r,i = − ω pi /ω , D i,e = D i,i = 0, and: D α = − ω pi k x Z dv y dv x ∂f α /∂v x v x − ω/k x − iν/k x , (5)where ω ps is the plasma frequency of species s , and ν → + determines the pole convention. We further take D r,α = 0, which is a good approximation when the alphaparticles are hot and sparse compared to the ions.Plugging these dispersion components into the real dis-persion Eq. (3), we find the familiar LH wave: ω r = ± ω LH ≡ ± (cid:0) ω − pi + | Ω e Ω i | − (cid:1) − / . (6)Taking ω r >
0, the imaginary dispersion yields: ω i = π S k x ω pα ω pi ω LH k x Z dv y ∂f α ∂v x | y = y ,v x = v p , (7)where S k x = sgn( k x ) determines the direction of thephase velocity v p ≡ ω r /k x .To recover alpha channeling, we need to transform thedistribution function and derivatives from phase spacecoordinates x i ≡ ( x, y, v x , v y ) to gyrocenter-energy coor-dinates X i ≡ ( X, Y, K, θ ): X = x + v y Ω α Y = y − v x Ω α (8) K = 12 m (cid:0) v x + v y (cid:1) θ = arctan ( − v y , v x ) . (9)In this coordinate system, the phase space density func-tion transforms as: F α = p | g | f α (10) | g | = | g ij | = (cid:12)(cid:12)(cid:12)(cid:12) ∂x m ∂X i ∂x n ∂X j δ mn (cid:12)(cid:12)(cid:12)(cid:12) = m − α . (11)Thus, we can rewrite our derivatives in terms of X i as: ∂f α ∂v x (cid:12)(cid:12)(cid:12)(cid:12) y , ωrkx = ∂X i ∂v x ∂∂X i ( m α F α ) y , ωrkx (12)= m α ω r k x (cid:18) ∂F α ∂K − k x m α ω r Ω α ∂F α ∂Y (cid:19) Y ∗ ,K ∗ . (13)where Y ∗ ≡ y − v p / Ω a , K ∗ ≡ m (cid:0) v p + v y (cid:1) /
2, and wehave taken F α independent of θ to capture the gyrophasediffusion induced by the stochasticity. Alpha Channeling: Slab B z k x v ⟂ v x ΔY FIG. 1. Schematic of the alpha channeling process. A hotparticle with v ⊥ > v p resonates with the wave at some point inthe orbit, leading to a change in both energy and gyrocenter.Thus the particles diffuse along a specified path in gyrocenter-energy space. The amplification condition Eq. (1) depends onthe derivative of the distibution function along this path. The quantity in parentheses in Eq. (13) can be recog-nized as derivative along the diffusion path in Eq. (1),and thus describes wave amplification from alpha chan-neling. Interestingly, this means that the condition thatthere be (on average) a population inversion along thechanneling diffusion path in gyrocenter-energy space isidentical to the condition that there be a bump-on-tailinstability in the local hot alpha particle distribution.Under this new formalism, in contrast to the nonlinearformalism, it is possible to straightforwardly calculate thewave amplification rate. For instance, for a Maxwellianwith a gradient in Y , F α = e − K/T α e − Y/L / πT α , withthermal velocity v thα = p T α /m α , we have at y = 0: ω i = −| ω LH | r π (cid:12)(cid:12)(cid:12)(cid:12) v px v thα (cid:12)(cid:12)(cid:12)(cid:12) e − (cid:16) vpxvthα (cid:17) × " ω pα ω pi (cid:18) − v thα v px ρ thα L (cid:19) e − ( y − ρ pα ) /L , (14)where ρ thα = v thα / Ω α , and ρ pα = v p / Ω α . Resonant Diffusion:
Identifying alpha channelingwith the unmagnetized bump-on-tail instability allows usto compactly derive the diffusion tensor, by performingthe same coordinate transformation. From unmagnetizedquasilinear theory, we have for a single wave mode [5, 13]: ∂f∂t = ∂∂x i (cid:18) D ij x ∂f α ∂v j (cid:19) (15) D v x v x x = ω pα m α n α W ( y ) k x i (cid:18) v x − ω r /k x − iω i /k x − i(cid:15) − c.c. (cid:19) , (16)where W ( y ) = E ( y ) / π is the wave electrostaticenergy at y . The diffusion equation (15) transforms as: ∂F α ∂t = ∂∂X i p | g | D ij X ∂∂X j F α p | g | !! , (17)where D ij X is determined from D ij x by the same tensortransformation law as for the metric in Eq. (11).Performing the coordinate transformation, taking F α independent of θ , and averaging Eq. (17) over θ , we find: (cid:28) ∂F α ∂t (cid:29) t = ∂∂ ¯ X i (cid:18)D D ij ¯X E θ ∂F α ∂ ¯ X j (cid:19) , (18)where the gyro-averaged coordinates ¯ X ≡ ( K, Y ), and: D D ij ¯X E θ = 12 π ω pα m α n α W ( y ) kv ⊥ i (cid:18) Ω − s I d − m α v ⊥ Ω α I d − m α v ⊥ Ω α I d m α v ⊥ I d (cid:19) (19)with v ⊥ = p K/m s , and I ad = I a − − I a + , where: I a ± ≈ (cid:18) ∓ iω i ∂∂ω r (cid:19) Z π dθ sin a θ sin θ − ω r /kv ⊥ ± i(cid:15) . (20)This integral can be evaluated with the u -substitution u = sin θ . We will focus on the resonant diffusion byignoring for now the nonresonant contribution from ω i ,which gives at Y = y − v p / Ω α : (cid:28) ∂F α ∂t (cid:29) t = ddK (cid:12)(cid:12)(cid:12)(cid:12) path " D KK ddK (cid:12)(cid:12)(cid:12)(cid:12) path F α (21) D KK ≡ m α (cid:18) q α E ( y ) m α (cid:19) v p p k v ⊥ − ω r H ( v ⊥ − v p )(22) ddK (cid:12)(cid:12)(cid:12)(cid:12) path ≡ (cid:18) ∂∂K − km s ω r Ω α ∂∂Y (cid:19) . (23)Eq. (22) is the same as Karney’s [9] diffusion coefficentin v ⊥ in the limit of large k x ρ as used in [2] (see supple-mental material). Furthermore, the diffusion is seen tooccur along the diffusion path in Eq. (1), confirming thatthis approach recovers alpha channeling.The diffusion coefficient corrects the energy-space dif-fusion coefficient in Ref. [1] and Ref. [14]. This discrep-ancy is discussed in the supplemental material. This er-ror did not affect the study of alpha channeling in toroidalgeometry due to ion-Bernstein waves (IBWs) [15, 16],which relied on a different diffusion coefficient from orbit-averaging the cyclotron-resonant response [17, 18]. Nonresonant Reaction:
Having established that thelinear-quasilinear system recovers alpha channeling, weare now in a position to examine the nonresonant re-sponse. In contrast to the resonant particles, which re-main on largely unperturbed gyro-orbits except at theresonance points and thus have a largely gyrotropic dis-tribution function (Fig. 2a), the nonresonant particlesexperience sloshing motion along v x only, and thus havea nongyrotropic distribution function at O ( E ) (Fig. 2b).Thus, instead of transforming the nonresonant diffu-sion coefficient to the coordinates X i , we find the nonres-onant response by first calculating the total force densityon species s from the field-particle correlation: F sx = q s h E x n i x . (24) FIG. 2. Simulated single-particle trajectories in the x - y planeof (a) hot particles and (b) cold particles (relative to the phasevelocity) in a growing electrostatic wave. The color of thetrajectory indicates time, light to dark. The cold particleshave a clearly non-gyrotropic velocity distribution due to theoscillations, and exhibit a clear shift in gyrocenter downward. This approach is equivalent to finding the force from thefull (non-gyro-averaged) quasilinear theory.Linearizing and Fourier transforming the above gives: F sx = q s L Z L/ − L/ dx Z dp x π dk x π ˜ E p x ˜ n s,k x × e i ( p x + k x ) x − i ( ω ( p x )+ ω ( k x )) t , (25)which can be expressed entirely in terms of ˜ φ by using˜ E p x = − ip x ˜ φ p x and ˜ n s,k x from Eq. (2). Then, in thisFourier convention, the wave φ = φ cos( k x x − ωt ) e ω i t corresponds to:˜ φ k x = πφ [ δ ( k x − k x ) + δ ( k x + k x )] . (26)Plugging this in to Eq. (25) and making use of the sym-metry property D ( − k x ) = D ∗ ( k x ) allows us to calculatethe total force on species s as: F sx = φ e ω i t π Im (cid:2) k x k D s (cid:3) (27) ≈ W k x (cid:20) D is + ω i ∂D rs ∂ω r (cid:21) . (28)Here, the first term on the RHS is the force density on theresonant particles, and the second term is the force den-sity on the nonresonant particles. This derivation gener-alizes the result for an unmagnetized plasma in Ref. [13]to any electrostatic wave. Summing over all species, werecover the imaginary component of the dispersion func-tion, Eq. (4). Thus the total force applied to the plasmasums to 0, as demanded by momentum conservation forthe electrostatic wave.The fact that the forces sum to zero in turn meansthat the total cross-field currents from resonant and non-resonant particles sum to zero, as can be seen by calcu-lating the total current from the resulting F × B drifts: X s j sy = X s − q s n s ( F sx /n s ) B z q s B z = − B z X s F sx = 0 . (29)Thus, for a purely poloidal wave mode growing in time,alpha channeling does not charge the plasma.In addition to revealing the conservation of totalcharge, Eq. (28) also tells us which species carries thecancelling nonresonant current. Specializing to an LHwave, with ∂D r,e /∂ω r = 0, we see that the nonresonantreaction is exclusively in the ions, which experience aforce density: F ix = 4 W k x ω i ω pi ω r = n i m i (cid:18) q i E m i (cid:19) k x ω i ω r . (30)Thus, for every alpha particle brought out of the plasmaby the LH alpha channeling instability, two fuel ions arebrought in, fueling the fusion reaction.The total shift in the nonresonant ion gyrocenter dueto the ponderomotive force for the LH wave can be ex-pressed nicely by integrating the F × B drift over thegrowth of the wave, using dE /dt = 2 ω i E . This gives:∆ Y = − q i m i (cid:18) k x ω r (cid:19) ∆ φ B z . (31)In a multi-ion-species plasma, Eq. (31) reveals to whatextent each ion species moves inward as a result of LHalpha channeling. For instance, in a p-B11 fusion plasma,Boron ions would pinch inward Z B /µ B = 5 /
11 as muchas the protons. For other electrostatic waves such as theIBW, the general force density from Eq. (28) can be usedto determine each species’ response. Eq. (31) can also beeasily checked against single-particle (Boris) simulationsin which a wave is ramped up from 0 initial amplitude,which are found to agree well (Fig. 3). Details for thesesimulations are given in the supplemental material.It is important to note that the cancellation of theresonant and nonresonant currents is not locally exact.Because the nonresonant ions are cold, the electric fieldthat enters this equation is evaluated locally at y ≈ Y ,in contrast to the case for the hot resonant particles,where it is evaluated at y = Y + ω r /k x Ω α . Thus, ifthere is variation in the electric field in y on some scalelength L , the slight offset of the resonant and nonresonantcurrents will produce a net current ordered down fromthe resonant current by O ( ρ pα /L ) (cid:28)
1. The resultingcharge accumulation could in principle drive shear flowin the plasma, albeit at a much reduced rate than thatsuggested by the resonant current alone.
Discussion:
The force in Eq. (30) is the same time-dependent force that arises from the unmagnetized pon-deromotive potential in the form:Φ = e E m ( ω − k · v ) , (32)from whence the force is derived via:( mδ ij − Φ v i v j ) dv j dt = Φ v i x j v j + Φ v i t − Φ x i . (33)where the subscripts represent derivatives. The force inEq. (30) appears as the second term on the RHS. FIG. 3. Change in gyrocenter position Y for the particle inFig. 2b due to the slow ramp-up of the electrostatic wave. Thegyro-period-averaged position in the simulations (filled blackdiamonds) agrees well with the predicted shift (gray dashedline) from Eq. (31) due to the nonresonant reaction. Eq. (32) can be derived from unmagnetized plasma sus-ceptibility using the K - χ theorem [19], which relates thelinear susceptibility to the ponderomotive potential. In-terestingly, application of the K - χ theorem to the mag-netized hot plasma dispersion relation [20] does not yieldthe nonresonant force we observe here, as we show in thesupplemental material. Nevertheless, single-particle sim-ulations using the Boris algorithm confirm the effect. Thefailure of the hot plasma dispersion to capture the effectis likely related to the gyro-average in the hot plasma dis-persion, which has been found to obscure the derivationof perpendicular resonant quasilinear forces [21–23].The approach used here, of taking the lowest-order al-pha particle motion to be a straight trajectory, and thenaveraging over a gyroperiod, is similar to how neoclas-sical wave-particle interactions are treated. In those in-teractions, one does not generally use the full constant-of-motion space dispersion relation [17] to calculate thequasilinear diffusion, but rather averages the effect of thediffusion derived from the magnetized dispersion rela-tion over the neoclassical orbit [18, 24]. This destroysresonances associated with the neoclassical orbit period,which are assumed to be destroyed by nonlinearities any-way. In each case, the long-term orbit is ignored in thecalculation of the dispersion, allowing in the neoclassi-cal case cyclotron damping for banana orbits, and in theLH alpha channeling case Landau damping at resonancepoints on the gyro-orbit.Note that, while the charge transport cancellation re-sult is general for any purely poloidal electrostatic wave,the channeling path in Eq. (1) and diffusion coefficientin Eq. (22) apply only to the case of gyro-averaged Lan-dau resonance [1, 2, 14, 25], and not to channeling viacyclotron resonances, as for the IBW [15, 16, 26–29]. Conclusion:
The alpha channeling interaction, whichreleases the free energy of particles through diffusionin coupled energy-space coordinates, can rigorously betransformed to the classic bump-on-tail instability in ve-locity space only. Applying the traditional mathematicalapparatus then shows that, in an initial value problem,where resonant ions are ejected as the electrostatic wavegrows at the expense of the ion energy, those same wavesmust pull in a return current of nonresonant ions so asto draw no current. This unexpected result is relatedto the cancellation of resonant and nonresonant currentsin the bump-on-tail instability [5, 13], except that thesenewly-found currents are perpendicular to the magneticfield, rather than parallel. We also calculated for thefirst time the contribution to the imaginary componentof the dispersion relation due to alpha channeling, a use-ful quantity for ray tracing calculations.We not only prove rigorously the current cancellation,but we also determine the extent to which each speciescontributes to this cancellation. For LH waves, the non-resonant ions are pulled into the plasma core; thus, whileno rotation is driven, the fusion reaction is beneficiallyfueled as ash is expelled. In p-B11 reactor, we showedthat protons are drawn in at twice the rate of Boron.While the nonresonant particles have been ignored inalpha channeling theory up to this point, our analysisshows that they can have important zeroth-order effects on the plasma dynamics. However, the specific prob-lem we considered here is only part of the story; in themost useful scenarios, channeling is driven by a station-ary wave propagating radially inwards from an antennaat the boundary, requiring a fundamentally 2D analy-sis. While the 2D problem is outside the scope of thisLetter, the 1D self-consistent theory of alpha channelinglaid out here provides a sound basis for examining thenonresonant response in more general scenarios.
Acknowledgments:
We would like to thank E.J.Kolmes and M.E. Mlodik for helpful discussions. Thiswork was supported by grants DOE DE-SC0016072 andDOE NNSA DE-NA0003871. One author (IEO) alsoacknowledges the support of the DOE ComputationalScience Graduate Fellowship (DOE grant number DE-FG02-97ER25308). [1] N. J. Fisch and J.-M. Rax, Physical Review Letters ,612 (1992).[2] N. J. Fisch and J. M. Rax, Nuclear Fusion , 549 (1992).[3] A. J. Fetterman and N. J. Fisch, Physical Review Letters , 205003 (2008).[4] A. N. Kaufman, J. Plasma Physics , 1 (1972).[5] N. A. Krall and A. W. Trivelpiece, Principles of PlasmaPhysics (McGraw-Hill, 1973).[6] A. I. Sukhorukov and P. Stubbe, Physics of Plasmas ,2497 (1997).[7] F. Charles, B. Despr´es, A. Rege, and R. Weder, (2020),arXiv:2002.11380.[8] C. F. F. Karney, Physics of Fluids , 1584 (1978).[9] C. F. F. Karney, Physics of Fluids , 2188 (1979).[10] R. Z. Sagdeev and V. D. Shapiro, ZhETF Pisma Redak-tsiiu , 389 (1973).[11] J. M. Dawson, V. K. Decyk, R. W. Huff, I. Jechart,T. Katsouleas, J. N. Leboeuf, B. Lembege, R. M. Mar-tinez, Y. Ohsawa, and S. T. Ratliff, Physical ReviewLetters , 1455 (1983).[12] G. Zaslavskii, M. Zakharov, R. Sagdeev, D. Usikov, andA. Chernikov, Sov. Phys. JETP , 294 (1986).[13] I. E. Ochs and N. J. Fisch, Physics of Plasmas , 062109(2020).[14] I. E. Ochs, N. Bertelli, and N. J. Fisch, Physics of Plas-mas , 112103 (2015).[15] M. C. Herrmann, Cooling Alpha Particles With Waves ,Ph.D. thesis, Princeton (1998). [16] D. S. Clark and N. J. Fisch, Physics of Plasmas , 2923(2000).[17] A. N. Kaufman, Physics of Fluids , 1063 (1972).[18] L. G. Eriksson and P. Helander, Physics of Plasmas ,308 (1994).[19] J. R. Cary and A. N. Kaufman, Physical Review Letters , 402 (1977).[20] T. H. Stix, Waves in Plasmas (AIP Press, 1992).[21] J. Lee, F. I. Parra, R. R. Parker, and P. T.Bonoli, Plasma Physics and Controlled Fusion (2012),10.1088/0741-3335/54/12/125005.[22] X. Guan, H. Qin, J. Liu, and N. J. Fisch, Physics ofPlasmas , 022502 (2013).[23] X. Guan, I. Y. Dodin, H. Qin, J. Liu, and N. J. Fisch,Physics of Plasmas , 102105 (2013).[24] M. C. Herrmann and N. J. Fisch, Physical Review Letters , 1495 (1997).[25] J. W. Cook, S. C. Chapman, and R. O. Dendy, PhysicalReview Letters , 1 (2010).[26] E. J. Valeo and N. J. Fisch, Phys. Rev. Lett. , 3536(1994).[27] F. Cianfrani and F. Romanelli, Nuclear Fusion ,106005 (2019).[28] C. Castaldo, A. Cardinali, and F. Napoli, PlasmaPhysics and Controlled Fusion , 084007 (2019).[29] F. Romanelli and A. Cardinali, Nuclear Fusion ,036025 (2020). upplementary Material: Nonresonant Diffusion in Alpha Channeling Ian E. Ochs and Nathaniel J. Fisch
Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08540, USA (Dated: December 14, 2020)
CONTENTS
I. Comparison of Resonant Perpendicular DiffusionCoefficients 1II. Particle Simulations 2III. K - χ theorem for hot magnetized plasma 2References 4 I. COMPARISON OF RESONANTPERPENDICULAR DIFFUSION COEFFICIENTS
The perpendicular diffusion equation given by Karney[1] in the highly stochastic limit is [see [1] Eqs. (71), (69),(3a-c)]: ∂f∂t = 1 v ⊥ ∂∂v ⊥ v ⊥ D v ⊥ v ⊥ ∂f∂v ⊥ (1) D v ⊥ v ⊥ = π E /B ) ω Ω i k ⊥ v ⊥ (cid:12)(cid:12)(cid:12) H (1) ν ( r ) (cid:12)(cid:12)(cid:12) (2) ν ≡ ω/ Ω i (3) r ≡ k ⊥ v ⊥ / Ω i . (4)It should be noted that this form of the diffusion equa-tion treats arises by treating f as a scalar function onphase space, so that probability integrals must includethe metric of the transformation, i.e. P ( v ⊥ − v ⊥ < ∆ v ⊥ ) = Z π dθ Z v ⊥ +∆ v ⊥ v ⊥ − ∆ v ⊥ dv ⊥ v ⊥ f ( v ⊥ ) . (5)This is in contrast to our approach in the main text,where the metric is included in the distribution trans-form.For large r/ν , i.e. large v ⊥ relative to v p = ω/k ⊥ , wehave [see [1] Eq. (68a)]: H (1) ν ( r ) ≈ (cid:18) π (cid:19) / ( r − ν ) − / (6) ≈ (cid:18) π (cid:19) / Ω / i ( k ⊥ v ⊥ − ω ) − / . (7)Since Karney works in SI, with Ω i = q i B /m i , plug-ging this approximation for ν into the diffusion coefficientyields: D v ⊥ v ⊥ = 12 (cid:18) q i E m i (cid:19) ( ω/k ⊥ v ⊥ ) ( k ⊥ v ⊥ − ω ) / . (8) This simplified form is used for the velocity-space coeffi-cient in the later channeling literature [2].To transform this form to energy space, we note thatthe factors of the metric in the diffusion equation arealready taken care of in the inclusion of the factors of v ⊥ in Eq. (1). Thus, all we need to do is transform thediffusion coefficient according to: D KK = (cid:18) ∂K∂v ⊥ (cid:19) D v ⊥ v ⊥ (9)= ( m i v ⊥ ) D v ⊥ v ⊥ (10)= m i (cid:18) q i E m i (cid:19) ( ω/k ⊥ ) ( k ⊥ v ⊥ − ω ) / . (11)This agrees with the form of the diffusion coefficient inthe main text.Ref. [3] also employs an energy-space diffusion opera-tor, which is used later in Ref. [4]. Focusing on energy-space, their diffusion equation was of the form: ∂F α ∂τ = ∂∂(cid:15) D (cid:15)(cid:15) ∂F α ∂(cid:15) (12) D (cid:15)(cid:15) = ( (cid:15) w ) a ων (cid:18) v osc v α (cid:19) (cid:15) − (cid:15) w ) / (13)with a = −
3, and: τ ≡ νt (14) (cid:15) ≡ v ⊥ v α (15) (cid:15) w ≡ (cid:18) ωk ⊥ v α (cid:19) (16) v osc ≡ q i Em i ω , (17)where ν is the collision frequency and v α is the alphaparticle birth energy.In terms of D (cid:15)(cid:15) , D KK is given by: D KK = (cid:18) ∂τ∂t (cid:19) (cid:18) ∂K∂(cid:15) (cid:19) D (cid:15)(cid:15) (18)= (2 ν ) (cid:18) m i v α (cid:19) D (cid:15)(cid:15) . (19)Thus, plugging in all the definitions, we find: D KK = m i v α (cid:18) ωk ⊥ v α (cid:19) a ω (cid:18) q i Em i ωv α (cid:19) v α ( v ⊥ − ω /k ⊥ ) / (20)= v − aα m i (cid:18) q i Em i (cid:19) ( ω/k ⊥ ) a − ( k ⊥ v ⊥ − ω ) / (21) a r X i v : . [ phy s i c s . p l a s m - ph ] D ec This agrees with the resonant diffusion coefficient in themain text if a = 3 / not if a = −
3, as used in Refs. [3, 4].
II. PARTICLE SIMULATIONS
The particle trajectories presented in the main textresult from Boris [5] simulations of the Lorentz force,given in SI by: m s d v dt = q s ( E + v × B ) . (22)We employ nondimensionalized units where q s = m s = 1.Futhermore, we take B = 1ˆ z , which normalizes t to thegyrofrequency Ω s .The electric field in the simulations is ramped up from0, according to: E = ˆ xk x φ sin ( k x x − ωt ) tanh (2 t/τ ) . (23)For the simulations presented, we had τ = 20 π , i.e. thefield ramped up over the course of 10 gyroperiods. Wealso had k x = 7 . ω = 20 .
12, with a phase velocity ω/k x = 2 .
55, and an amplitude φ = 0 . x = π/k x , y = 0. For the simulation in Fig. 2a, we had v = (5 , v = (0 , . x k +1 − x k ∆ t = v k +1 (24) v k +1 − v k ∆ t = q s m s (cid:18) E k + ( v k +1 + v k ) × B k (cid:19) , (25)where x = x ( t k ), v k = v ( t k − ∆ t/ t k = k ∆ t , E k = E ( x k , t k ), and B k = B ( x k , t k ).This procedure works when ∆ t is constant; however,when ∆ t = ∆ t k is not constant, then the electric fieldwill no longer be evaluated exactly between v k and v k +1 ;instead, v k = v ( t k − ∆ t k − /
2) and v k +1 = v ( t k +∆ t k / x k = x ( t k ).To get the adaptive timestep to work, we split the x step in two. This has the added benefit that output( x k , v k ) from the code represents x and v evaluated at the same timepoint t k . Thus, we define: x k = x ( t k ), v k = v ( t k ). Then, we have: x k +1 / − x k ∆ t k / v k (26) v k +1 − v k ∆ t k = q s m s (cid:18) E k +1 / + ( v k +1 + v k ) × B k +1 / (cid:19) (27) x k +1 − x k +1 / ∆ t k / v k +1 , (28)with E k +1 / = E ( x k +1 / , t k + ∆ t k / B k +1 / = B ( x k +1 / , t k + ∆ t k / x : x k +3 / − x k +1 / (∆ t k + ∆ t k +1 ) / v k +1 , (29)which is the natural generalization of Eq. (24) to an adap-tive timestep, given the redefinition of x k to line up with v k . In effect, we have exchanged asymmetry in the v -stepdue to the adaptive timestep for asymmetry in the x -step,which in simple tests seems to have a less deleteriouseffect on the dynamics, allowing ∼
10x larger timestepsthan simply plugging a variable timestep into the Borisalgorithm.Our adaptive timestep itself was determined by:∆ t k = 120 min (cid:16) Ω − i , ω − , ( k x v k ) − , p /k x a k (cid:17) , (30)with a k the acceleration; this meant that the timestepwas always significantly less than that required for theparticle to traverse the time or length scales of the prob-lem. III. K - χ THEOREM FOR HOT MAGNETIZEDPLASMA
The K - χ theorem [7] relates the ponderomotive poten-tial ( K ) to the dielectric susceptibility ( χ ) of a plasma.The K - χ theorem is written: Z d x d v f s Φ = − π Z d x Z d x E ∗ ω ( x ) χ ω ( x , x ) E ω ( x ) , (31)where E ( x , t ) = E ω ( x ) e − iωt + c.c. (32)Note that this relates to our variables via | E ω | = | E | /
4, yielding an extra factor of 1 /
4, and uses a dif-ferent Fourier convention (in terms of factors of 2 π ) aswell. Furthermore, f s is here defined so that R d v f s = n s ;changing this normalization condition to equal 1 yieldsan extra factor of n s on the LHS. Taking a homogenousplasma so that χ ω ( x , x ) = χ ω ( x − x ), and consideringa Fourier mode, we thus can write the K - χ theorem as: n s Z d v f s Φ = − E π e ∗ ωk χ ωk e ωk , (33)where e ωk is the normalized polarization vector. Fromthis form, the form in the main text quickly follows:Φ = − Wn s e ∗ · δχ s δf s · e . (34)As an example, the Vlasov susceptibility for a hot, un-magnetized plasma is [8]: χ s = − ω ps k Z d v k i ∂f s /∂v i k j v j − ω (35)= − Z d v f s ω ps ( k m v m − ω ) k i k j k . (36)Furthermore, for an electrostatic wave, the polarizationvector is given by e i = k i /k . Thus:Φ = − Wn s k i k − ω ps ( k m v m − ω ) k i k j k ! k j k (37)= q s E m s ( ω − k · v ) . (38)Now, we will apply the K - χ theorem to the hot magne-tized susceptibility [8], with k k = 0. Since we are lookingat electrostatic waves, we will have e = ˆ x , and we willonly care about the χ xx component of the susceptibilitytensor. We will examine the lowest-order thermal cor-rection, where the time-dependent ponderomotive forcemakes its appearance in the unmagnetized case.The hot magnetized susceptibility in this case is givenby (suppressing species labels): χ xx = ω p ω Ω ∞ X n = −∞ Z ∞ πv ⊥ dv ⊥ Ω ω − n Ω S n,xx (39) S n,xx = n z J n ( z ) v ⊥ ∂f ∂v ⊥ (40) z ≡ k ⊥ v ⊥ Ω . (41)We can express this integral entirely in terms of z bytaking noting that dv ⊥ = (Ω /k ⊥ ) dz , and ∂f /∂v ⊥ =( k ⊥ / Ω) ∂f /∂z , so plugging this in above and simplifyingwe find: χ xx = (cid:18) Ω k ⊥ (cid:19) ω p ω Ω 2 π Z ∞ dz ∞ X n = −∞ n a − n J n ( z ) ∂f ∂z , (42)where a = ω/ Ω. We can make use of the Newberger sum rule to derivethe identity [see Ref. [9] pg. 149]: ∞ X n = −∞ n a − n J n ( z ) = πa sin πa J a ( z ) J − a ( z ) − a. (43)Further expanding this in k ⊥ v ⊥ /ω = z/a (cid:28)
1, with a >
1, we find: ∞ X n = −∞ n a − n J n ( z ) ≈ az a −
1) + 38 az a − a + 4 + O (cid:16) za (cid:17) . (44)We plug this expansion in to our integral and integrateby parts to find: χ xx = (cid:18) Ω k ⊥ (cid:19) ω p ω Ω 2 π Z ∞ dz (cid:20) − aza − − az a − a + 4 (cid:21) f . (45)Then, plugging back in for a and z , we find: χ xx = Z ∞ (2 πv ⊥ dv ⊥ ) f × (cid:18) − ω p ω − Ω − ω p ω k ⊥ v ⊥ ω − + 4Ω /ω (cid:19) . (46)In this expression, the first term in brackets is the metricfor the integral over phase space, after the integral over θ . Thus: δχδf = − ω p ω − Ω − ω p ω k ⊥ v ⊥ ω − + 4Ω /ω . (47)Consider first the leading order term, Φ . Plugging into the K - χ theorem Eq. (33), we have:Φ = q E m ( ω − Ω ) . (48)This is the familiar ponderomotive potential for a coldplasma in the electrostatic limit [7]. It yields no x -directed force due to the growth of the field, sinceΦ ,v x t = 0.Next, we have the second-order term:Φ = q E m ( ω − + 4Ω /ω ) k ⊥ v ⊥ ω . (49)Since v ⊥ = v x + v y , the x -directed time-dependent forcedensity is given from this ponderomotive potential by: F x = n Φ v x t (50)= q n m ( ω − + 4Ω /ω ) k ⊥ v x ω dE dt (51)= " nm (cid:18) qE m (cid:19) k ⊥ ω i ω ωk ⊥ v x ( ω − + 4Ω /ω ) , (52)where in the last line we used dE /dt = 2 ω i E . The firstterm in brackets here is the force density from the mainpaper, which was found to be consistent with the single-particle simulations. Thus, the time-dependent pondero-motive force from the K - χ theorem is ordered down (inthe relevant limit ω (cid:29) Ω) by a factor v ⊥ /v p (cid:28) K - χ theorem in this case comes from the fact that the susceptibility tensor isgyro-averaged prior to application of the theorem, thuslosing the angyrotropy of the velocity distribution thatgives rise to the nonresonant reaction. Such gyrophase-dependent structure has been shown to also be importantin evaluating perpendicular forces in resonant diffusionin magnetized plasmas [10–12], so it makes sense that itcould effect the ponderomotive forces as well. [1] C. F. F. Karney, Physics of Fluids , 2188 (1979).[2] N. J. Fisch and J. M. Rax, Nuclear Fusion , 549 (1992).[3] N. J. Fisch and J.-M. Rax, Physical Review Letters ,612 (1992).[4] I. E. Ochs, N. Bertelli, and N. J. Fisch, Physics of Plas-mas , 112103 (2015).[5] J. P. Boris, in Proc. Fourth Conf. Num. Sim. Plasmas,Naval Res. Lab, Wash. DC (1970) pp. 3–67.[6] H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, and W. M.Tang, Physics of Plasmas , 084503 (2013).[7] J. R. Cary and A. N. Kaufman, Physical Review Letters , 402 (1977).[8] T. H. Stix, Waves in Plasmas (AIP Press, 1992).[9] D. G. Swanson,
Plasma waves (Elsevier, 2012).[10] J. Lee, F. I. Parra, R. R. Parker, and P. T.Bonoli, Plasma Physics and Controlled Fusion (2012),10.1088/0741-3335/54/12/125005.[11] X. Guan, H. Qin, J. Liu, and N. J. Fisch, Physics ofPlasmas , 022502 (2013).[12] X. Guan, I. Y. Dodin, H. Qin, J. Liu, and N. J. Fisch,Physics of Plasmas20