A Quasi-Linear Diffusion Model for Resonant Wave-Particle Instability in Homogeneous Plasma
Seong-Yeop Jeong, Daniel Verscharen, Robert T. Wicks, Andrew N. Fazakerley
DDraft version October 22, 2020
Typeset using L A TEX twocolumn style in AASTeX63
A Quasi-Linear Diffusion Model for Resonant Wave-Particle Instability in Homogeneous Plasma
Seong-Yeop Jeong, Daniel Verscharen,
1, 2
Robert T. Wicks,
1, 3 and Andrew N. Fazakerley Mullard Space Science Laboratory, University College London, Dorking, RH5 6NT, UK; [email protected], [email protected] Space Science Center, University of New Hampshire, Durham, NH 03824, USA Institute for Risk and Disaster Reduction, University College London, Gower Street, London, WC1E 6BT, UK
AbstractIn this paper, we develop a model to describe the generalized wave-particle instability in a quasi-neutral plasma. We analyze the quasi-linear diffusion equation for particles by expressing an arbitraryunstable and resonant wave mode as a Gaussian wave packet, allowing for an arbitrary direction ofpropagation with respect to the background magnetic field. We show that the localized energy densityof the Gaussian wave packet determines the velocity-space range in which the dominant wave-particleinstability and counter-acting damping contributions are effective. Moreover, we derive a relationdescribing the diffusive trajectories of resonant particles in velocity space under the action of suchan interplay between the wave-particle instability and damping. For the numerical computation ofour theoretical model, we develop a mathematical approach based on the Crank-Nicolson schemeto solve the full quasi-linear diffusion equation. Our numerical analysis solves the time evolutionof the velocity distribution function under the action of a dominant wave-particle instability andcounteracting damping and shows a good agreement with our theoretical description. As an application,we use our model to study the oblique fast-magnetosonic/whistler instability, which is proposed asa scattering mechanism for strahl electrons in the solar wind. In addition, we numerically solve thefull Fokker-Planck equation to compute the time evolution of the electron-strahl distribution functionunder the action of Coulomb collisions with core electrons and protons after the collisionless action ofthe oblique fast-magnetosonic/whistler instability.
Keywords: methods: analytical — instabilities — waves — plasmas — diffusion INTRODUCTIONWave-particle resonances play an important role for theenergy exchange between particles and waves in manyspace and astrophysical plasmas. For example, wave-particle resonances contribute to the acceleration anddeceleration of particles in radiation belts (Ukhorskiy &Sitnov 2014), the deviation of the particle velocity distri-bution function (VDF) from a Maxwellian equilibrium inthe solar wind (Marsch 2006), the thermodynamic stateof the intracluster medium in galaxy clusters (Roberg-Clark et al. 2016), and the scattering and absorptionof the surface radiation in neutron-star magnetospheres(Lyutikov & Gavriil 2006). Therefore, it is of great impor-tance to study the mechanics of wave-particle resonancesin order to advance our understanding of the physics ofastrophysical plasmas throughout the universe.According to kinetic theory, wave-particle resonancescan occur in the form of Landau or cyclotron resonances,which contribute to wave instability or wave dampingdepending on the resonance’s characteristics. The quasi- linear theory of magnetized plasma, first established byYakimenko (1963) and Kennel & Engelmann (1966), pro-vides a mathematical framework to predict the evolutionof the particle VDF under the action of wave-particle res-onances. Quasi-linear theory assumes that the spatiallyaveraged VDF evolves slowly compared to the gyrope-riod of the particles and the wave period. It furthermoreassumes that the fluctuation amplitude is small and thatthe spatial average of the fluctuations vanishes. Basedon this theory, numerous analytical studies have suc-cessfully explained the evolution of VDFs resulting fromwave-particle resonances.Resonant particles diffuse along specific trajectories invelocity space determined by the properties of resonantwave (Kennel & Engelmann 1966; Gendrin 1968, 1981;Gendrin & Roux 1980; Stix 1992; Isenberg & Lee 1996;Summers et al. 1998, 2001). In these models, quasi-lineardiffusion coefficients determine the diffusion rate of theresonant particles (Lyons et al. 1971; Lyons 1974; Albert2004; Glauert & Horne 2005; Summers 2005; Isenberg& Vasquez 2011; Tang et al. 2020). Alternatively, quasi- a r X i v : . [ phy s i c s . s p ace - ph ] O c t linear diffusion models based on a bi-Maxwellian VDF,in which only its moments evolve in time, describe theeffective evolution of particle VDFs under the action ofmicroinstabilities (Seough & Yoon 2012; Yoon & Seough2012; Yoon et al. 2015, 2017; Yoon 2017). Moreover,fully nonlinear simulations based on kinetic theory modelthe evolution of the particle VDF consistently with pre-dictions from quasi-linear theory (Vocks & Mann 2003;Vocks et al. 2005; Gary et al. 2008; Saito et al. 2008, 2010;Saito & Peter Gary 2012). Observations from Helios re-vealed signatures in the proton VDFs consistent with ioncyclotron resonances as predicted by quasi-linear theory(Marsch & Tu 2001; Tu & Marsch 2002; Heuer & Marsch2007; Marsch & Bourouaine 2011).Realistic analytical models must describe the diffusivetrajectory of the resonant particles in velocity space,taking into account the localized (in wavevector space)energy density of the waves that resonate with these parti-cles. These models must also account for non-Maxwellianfeatures in the VDF evolution in order to advance ourunderstanding of plasma observations and kinetic sim-ulation results. A rigorous numerical analysis of thediffusion equation, including both the diagonal and off-diagonal diffusion terms, is necessary to support thetheoretical description through the quantification of thediffusion rates.By analyzing the quasi-linear diffusion equation, wepropose a novel quasi-linear diffusion model for the timeevolution of VDFs under the action of a dominant wave-particle instability and counter-acting damping contri-butions in Section 2. Our model describes the creationand evolution of non-Maxwellian features in the particleVDF. We allow for an arbitrary type of the unstable andresonant wave mode with an arbitrary direction of prop-agation with respect to the background magnetic field.In our analysis, we express the electric field of this waveas a Gaussian wave packet in configuration space. Thelocalization of such a wave packet in configuration spaceis the direct consequence of its generation through alinear instability, which is localized in wavevector space.To investigate the stabilization of the VDF throughquasi-linear diffusion, we apply our analysis of the quasi-linear diffusion equation to Boltzmann’s H -theorem. Inthis scheme, the localized energy density of the Gaussianwave packet in wavevector space defines the velocity-space range in which the dominant wave-particle in-stability and counter-acting damping contributions areeffective. In addition, we derive a relation to describethe diffusive trajectories of resonant particles in velocityspace under the action of such an instability and damping.In this way, our model accounts for the diffusive behaviorof resonant particles in different regions of velocity space. For the numerical evaluation of our theoretical de-scription, we develop a mathematical approach basedon the Crank-Nicolson scheme (for numerical details,see Appendix A) that solves the full quasi-linear dif-fusion equation. Because of its reliable stability, theCrank-Nicolson scheme has been used previously to solvediffusion equations in a variety of fields (Khazanov et al.2002; Albert 2004; Brügmann et al. 2004; Yang et al.2009; Klein & Chandran 2016; Taran et al. 2019). How-ever, most standard Crank-Nicolson schemes ignore theoff-diagonal terms in the diffusion equation. In our case,these off-diagonal terms are important for the descrip-tion of resonant pitch-angle scattering. We note thatour mathematical approach is applicable to all generaltwo-dimensional diffusion equations, including those withoff-diagonal diffusion terms.In Section 3, as an example, we apply our model to thescattering of the electron strahl, which is a field-alignedelectron beam population in the solar wind (Pilipp et al.1987; Štverák et al. 2009). Observations in the solarwind suggest that strahl electrons exchange energy withwhistler waves, which ultimately leads to a scatteringof strahl electrons into the halo population (Pagel et al.2007; Lacombe et al. 2014; Graham et al. 2017). Ourquasi-linear framework confirms that an instability of thefast-magnetosonic/whistler wave in oblique propagationwith respect to the background magnetic field scattersthe electron strahl into the electron halo, as predicted bylinear theory (Vasko et al. 2019; Verscharen et al. 2019).In Section 4, for a more realistic model of the strahlevolution after the collisionless action of the oblique fast-magnetosonic/whistler instability, we numerically solvethe full Fokker-Planck equation for Coulomb collisionswith our mathematical approach (for numerical details,see Appendix B). We model the time evolution of theelectron-strahl VDF through the action of Coulomb col-lisions with core electrons and protons. This combinedmethod allows us to compare the timescales for the strahlscattering and collisional relaxation.In Section 5, we discuss the results of our model for thestrahl scattering and electron-halo formation throughthe instability and Coulomb collisions. In Section 6, wesummarize and conclude our treatment. QUASI-LINEAR DIFFUSION MODELIn this section, we establish our general theoreticalframework for the description of a resonant wave-particleinstability in quasi-linear theory. Because our work fo-cuses on non-relativistic space plasma like the solar wind,we ignore relativistic effects throughout our study.2.1.
Analysis of the Quasi-Linear Diffusion Equation
To investigate the time evolution of the particle VDFthrough wave-particle resonances, we study the quasi-linear diffusion equation, given by Stix (1992): (cid:18) ∂f j ∂t (cid:19) QLD = lim V →∞ ∞ (cid:88) n = −∞ (cid:90) πq j V m j × ˆ G [ k (cid:107) ] (cid:34) v ⊥ (cid:12)(cid:12) v (cid:107) (cid:12)(cid:12) δ (cid:18) k (cid:107) − ω k − n Ω j v (cid:107) (cid:19)(cid:12)(cid:12) ψ nj (cid:12)(cid:12) ˆ G [ k (cid:107) ] f j (cid:35) d k , (1)where ψ nj ≡ √ (cid:20) E Rk e iφ J n +1 ( ρ j ) + E Lk e − iφ J n − ( ρ j ) (cid:21) + v (cid:107) v ⊥ E zk J n ( ρ j ) , (2)and ˆ G [ k (cid:107) ] ≡ (cid:18) − k (cid:107) v (cid:107) ω k (cid:19) v ⊥ ∂∂v ⊥ + k (cid:107) ω k ∂∂v (cid:107) . (3)The integer n determines the order of the resonance,where n = 0 corresponds to the Landau resonance and n (cid:54) = 0 corresponds to cyclotron resonances. In our equa-tions, we label contributions from a given resonanceorder with a superscript n . The subscript j indicatesthe particle species. The particle VDF of species j isdenoted as f j ≡ f j ( v ⊥ , v (cid:107) , t ) which is spatially averagedand gyrotropic, q j and m j are the charge and mass ofa particle of species j , v ⊥ and v (cid:107) are the velocity coor-dinates perpendicular and parallel with respect to thebackground magnetic field. We choose the coordinatesystem in which the proton bulk velocity is zero. Wedenote the n th-order Bessel function as J n ( ρ j ) where ρ j ≡ k ⊥ v ⊥ / Ω j . The cyclotron frequency of species j isdefined as Ω j ≡ q j B /m j c , B is the background mag-netic field, c is the speed of light, k ⊥ and k (cid:107) are theperpendicular and parallel components of the wavevector k , and V is the volume in which the wave amplitude iseffective so that the wave and particles undergo a signifi-cant interaction. We denote Dirac’s δ -function as δ andthe azimuthal angle of wavevector k by φ . The frequency ω is a complex function of k , and we define ω k as itsreal part and γ k as its imaginary part ( ω = ω k + iγ k ).Without loss of generality, we set ω k > . Furthermore,we assume that | γ k | (cid:28) ω k , i.e. the assumption of slowgrowth or damping that is central to quasi-linear theory.The spatially Fourier-transformed electric field has theform of E k = ˆ xE xk + ˆ yE yk + ˆ zE zk and is defined as (Gurnett& Bhattacharjee 2017) E k = 1(2 π ) / (cid:90) E r exp [ − i k · r ] d r , (4) where k · r = k ⊥ x cos φ + k ⊥ y sin φ + k (cid:107) z and r is theposition vector. We take the constant background mag-netic field as B = ˆ zB and define the right- and left-circularly polarized components of the electric field as E Rk ≡ ( E xk − i E yk ) / √ and E Lk ≡ ( E xk + i E yk ) / √ . Thelongitudinal component of the electric field is E zk .Linear instabilities typically create fluctuations acrossa finite range of wavevectors. The Fourier transformationof such a wave packet in wavevector space correspondsto a wave packet in configuration space. For the sake ofsimplicity, we model this finite wave packet by assumingthat the electric field E r of the unstable and resonantwaves has the shape of a gyrotropic Gaussian wave packet E r = E exp (cid:34) − σ ⊥ x + σ ⊥ y + σ (cid:107) z (cid:35) exp [ i k · r ] , (5)where E = ˆ xE x + ˆ yE y + ˆ zE z , k · r = k ⊥ x cos φ + k ⊥ y sin φ + k (cid:107) z , and k is the wavevector of the Gaus-sian wave packet. We allow for an arbitrary angle θ between k and B , which defines the orientation ofthe wavevector at maximum growth of the wave, andassume that k (cid:107) (cid:54) = 0 . The vector E represents the peakamplitude of the electric field. The free parameters σ ⊥ and σ (cid:107) characterize the width of the Gaussian envelope.Quasi-linear theory requires that E r spatially averagesto zero. Therefore, we assume that | k (cid:107) | (cid:29) σ (cid:107) so thatthe spatial dimension of the Gaussian wave packet islarge compared to the parallel wavelength π/ | k (cid:107) | .The spatial Fourier transformation of Eq. (5) accordingto Eq. (4) then leads to E k = E σ (cid:107) σ ⊥ exp (cid:34) − ( k (cid:107) − k (cid:107) ) σ (cid:107) − ( k ⊥ − k ⊥ ) σ ⊥ (cid:35) . (6)We identify V with the volume of the Gaussian envelope, V = 1 / ( σ (cid:107) σ ⊥ ) . Eq. (6) describes the localization ofthe wave energy density in wavevector space. For theinstability analysis through Eq. (6), we define the unsta-ble k -spectrum as the finite wavevector range in which γ k > and argue that resonant waves exist only in thisunstable k -spectrum. We ignore any waves outside this k -spectrum since they are damped.We define k (cid:107) as the value of k (cid:107) at the center of theunstable k spectrum. We then obtain k ⊥ = k (cid:107) tan θ . (7)In the case of a linear plasma instability, we identify k ⊥ and k (cid:107) with the wavevector components at which theinstability has its maximum growth rate as a reasonableapproximation. To approximate the wave frequency ofthe unstable waves at the angle θ of maximum growth,we expand ω k of the unstable and resonant waves around k (cid:107) as ω k ( k (cid:107) ) ≈ ω k + v g (cid:0) k (cid:107) − k (cid:107) (cid:1) , (8)where v g ≡ ∂ω k ∂k (cid:107) (cid:12)(cid:12)(cid:12)(cid:12) k (cid:107) = k (cid:107) . (9)In Eqs. (8) and (9), ω k and v g are the wave frequencyand parallel group velocity of the unstable and resonantwaves, evaluated at k (cid:107) = k (cid:107) . We select the values of σ ⊥ and σ (cid:107) as the half widths of the perpendicularand parallel unstable k -spectrum. In the case of a linearplasma instability, the numerical values for k ⊥ , k (cid:107) , σ ⊥ , σ (cid:107) , ω k and v g can be found from the solutions of thehot-plasma dispersion relation, which thus closes our setof equations.By using Eq. (6) and Eq. (8), we rewrite Eq. (1) as (cid:18) ∂f j ∂t (cid:19) QLD = ∞ (cid:88) n = −∞ (cid:90) ˆ G [ k (cid:107) ] (cid:2) D nj ˆ G [ k (cid:107) ] f j (cid:3) d k , (10)where D nj ≡ πq j v ⊥ σ (cid:107) σ ⊥ m j δ ( k (cid:107) − k n (cid:107) j ) × | ψ nj | | v (cid:107) − v g | exp (cid:34) − ( k (cid:107) − k (cid:107) ) σ (cid:107) − ( k ⊥ − k ⊥ ) σ ⊥ (cid:35) , (11) ψ nj ≡ √ (cid:20) E R e iφ J n +1 ( ρ j ) + E L e − iφ J n − ( ρ j ) (cid:21) + v (cid:107) v ⊥ E z J n ( ρ j ) , (12)and k n (cid:107) j ≡ ω k − k (cid:107) v g − n Ω j v (cid:107) − v g . (13)We set E R = ( E x − i E y ) / √ and E L = ( E x + i E y ) / √ as constant, evaluated at k .Eq. (10) is the quasi-linear diffusion equation describingthe action of the dominant wave-particle instability andco-existing damping contributions from other resonancesin a Gaussian wave packet. We define the n resonance asthe contribution to the summation in Eq. (10) with onlyinteger n . We note that any n resonance can contributeto wave instability or to wave damping depending on theresonance’s characteristics.2.2. Stabilization through a Resonant Wave-ParticleInstability
We define stabilization as the process that createsthe condition in which ( ∂f j /∂t ) QLD → for all v ⊥ and v (cid:107) . For our analysis of the stabilization of a VDF through a resonant wave-particle instability, including co-existing damping effects, we use Boltzmann’s H -theorem,in which the quantity H is defined as H ( t ) ≡ (cid:90) f j ( v , t ) ln f j ( v , t ) d v . (14)By using Eq. (10), the time derivative of H is given by dHdt = ∞ (cid:88) n = −∞ (cid:90) (cid:90) (ln f j +1) ˆ G [ k (cid:107) ] (cid:2) D nj ˆ G [ k (cid:107) ] f j (cid:3) d k d v . (15)The integrand in Eq. (15) is equivalent to ( ln f j + 1) ˆ G [ k (cid:107) ] (cid:2) D nj ˆ G [ k (cid:107) ] f j (cid:3) =ˆ G [ k (cid:107) ] (cid:2) D nj ˆ G [ k (cid:107) ]( f j ln f j ) (cid:3) − D nj (cid:2) ˆ G [ k (cid:107) ] f j (cid:3) (cid:14) f j . (16)Upon substituting Eq. (16) into Eq. (15), the first termon the right-hand side in Eq. (16) disappears after theintegration over v . Then, by resolving the δ -function in D nj through the k (cid:107) -integral, we obtain dHdt = − ∞ (cid:88) n = −∞ (cid:18) dHdt (cid:19) n , (17)where (cid:18) dHdt (cid:19) n ≡ (cid:90) (cid:110) (cid:101) D nj (cid:2) ˆ G [ k n (cid:107) j ] f j (cid:3) (cid:14) f j (cid:111) d v , (18) (cid:101) D nj ≡ W nj πq j v ⊥ σ (cid:107) σ ⊥ m j (cid:90) π (cid:90) ∞ | ψ nj | × exp (cid:34) − ( k ⊥ − k ⊥ ) σ ⊥ (cid:35) k ⊥ dk ⊥ dφ, (19) W nj ≡ | v (cid:107) − v g | exp (cid:34) − k (cid:107) σ (cid:107) (cid:18) v (cid:107) − v n (cid:107) res v (cid:107) − v g (cid:19) (cid:35) , (20) v n (cid:107) res ≡ ω k − n Ω j k (cid:107) , (21) ˆ G [ k n (cid:107) j ] ≡ (cid:20) n Ω j ω k − k (cid:107) v g − n Ω j (cid:21) v (cid:107) − v g v ph v ⊥ ∂∂v ⊥ + 1 v ph ∂∂v (cid:107) , (22)and v ph ≡ ω k ( k n (cid:107) j ) k n (cid:107) j = ( ω k − k (cid:107) v g ) v (cid:107) − n Ω j v g ω k − k (cid:107) v g − n Ω j . (23)The function (cid:101) D nj in Eq. (19) plays the role of a diffusioncoefficient for the n resonance. In (cid:101) D nj , the v (cid:107) -function W nj defined in Eq. (20) serves as a window function thatdetermines the region in v (cid:107) -space in which the quasi-linear diffusion through the n resonance is effective. Thewindow function W nj is maximum at v n (cid:107) res defined inEq. (21), which is the parallel velocity of the particlesthat resonate with the waves at k (cid:107) = k (cid:107) through the n resonance. Our window function W nj is linked to Dirac’s δ -function in the limit lim v n (cid:107) res → v g W nj ≈ √ π σ (cid:107) | k (cid:107) | δ ( v (cid:107) − v g ) , (24)where | k (cid:107) | (cid:29) σ (cid:107) . Through this ordering between | k (cid:107) | and σ (cid:107) , we assume that W nj restricts a finite region in v (cid:107) -space and that the W nj for different resonances do notoverlap with each other in v (cid:107) -space.Only particles distributed within W nj experience the n resonance and contribute to the quasi-linear diffusion,which is ultimately responsible for the stabilization. Be-cause all terms in Eq. (18) are positive semi-definite,all resonances independently stabilize f j through quasi-linear diffusion in the v (cid:107) -range defined by their respective W nj , according to Eq. (17). Therefore, H decreases and dH/dt tends toward zero during the quasi-linear diffu-sion through all resonances while f j is in the process ofstabilization. When f j reaches a state of full stabilizationthrough all n resonances, the instability has saturatedand its growth ends.The v (cid:107) -function k n (cid:107) j defined in Eq. (13) is the resonantparallel wavenumber, fulfilling the condition that k n (cid:107) j = k (cid:107) at v (cid:107) = v n (cid:107) res . It quantifies the k (cid:107) -component ofthe unstable k -spectrum in the v (cid:107) -range defined by W nj .Eq. (23) defines the phase velocity at k n (cid:107) j , which is onlyconstant when v g = ω k /k (cid:107) , in which case v ph = v g for all v (cid:107) . We discuss the diffusion operator ˆ G [ k n (cid:107) j ] inEq. (22) in the next section.2.3. Nature of Quasi-Linear Diffusion in Velocity Space
According to Eq. (18), unless the wave amplitude iszero, the condition for achieved stabilization through the n resonance is ˆ G [ k n (cid:107) j ] F nj ( v ⊥ , v (cid:107) ) = 0 , (25)where F nj ( v ⊥ , v (cid:107) ) represents the stabilized VDF of species j through the n resonance. In Eq. (25), ˆ G [ k n (cid:107) j ] is a di-rectional derivative along the isocontour of F nj evaluatedat a given velocity position. Considering the role of W nj , ˆ G [ k n (cid:107) j ] describes only the diffusion of resonant particleswithin W nj along the isocontour of F nj . Consequently, theparticles experiencing the n resonance diffuse toward the stable state so that ( dH/dt ) n → , while the isocontoursof F nj describe the diffusive velocity-space trajectories forthe n resonance.To find such a trajectory, we express an infinitesimalvariation of F nj along an isocontour as d F nj = ∂ F nj ∂v ⊥ dv ⊥ + ∂ F nj ∂v (cid:107) dv (cid:107) = 0 . (26)Eqs. (22) and (26) allow us to rewrite Eq. (25) as v ⊥ dv ⊥ + (cid:20) n Ω j n Ω j − ω k + k (cid:107) v g (cid:21) ( v (cid:107) − v g ) dv (cid:107) = 0 . (27)By integrating Eq. (27), the diffusive trajectory for the n resonance is then given by v ⊥ + (cid:20) n Ω j n Ω j − ω k + k (cid:107) v g (cid:21) ( v (cid:107) − v g ) = const. (28)Kennel & Engelmann (1966) treat the two limiting casesin which v g = ω k /k (cid:107) and v g = 0 . Using their assump-tions, our Eq. (28) is equivalent to their equation (4.8)if v g = ω k /k (cid:107) , and our Eq. (28) is equivalent to theirequation (4.11) if v g = 0 . Depending on the dispersionproperties of the resonant waves, Eq. (28) is either anelliptic or a hyperbolic equation when n (cid:54) = 0 . In the caseof electron resonances, it is safe to assume that n Ω j n Ω j − ω k + k (cid:107) v g ≥ , (29)in Eq. (28) if v g < ( ω k + n | Ω e | ) /k (cid:107) for all positive n and v g > ( ω k + n | Ω e | ) /k (cid:107) for all negative n . However,in the case of proton resonances, resonant waves are morelikely to violate Eq. (29) because Ω p (cid:28) | Ω e | .Fig. 1 illustrates the diffusive flux of particles experienc-ing two arbitrary resonances: the n and n resonancesfor an unstable wave. The dark shaded areas representisocontours of the VDFs of two particle populations invelocity space. The red and dark-blue solid curves rep-resent the diffusive trajectories according to Eq. (28)with n = n and n = n , assuming that the resonantwave fulfills Eq. (29). The window functions W n j and W n j describe the v (cid:107) -ranges in which the n and n reso-nances are effective. The light-blue dashed semi-circlescorrespond to contours of constant kinetic energy in theproton rest frame, for which v ⊥ + v (cid:107) = const. (30)In general, the diffusive flux is always directed fromhigher to lower phase-space densities during the processof stabilization. At point A, resonant particles in W n j diffuse along the red solid curve toward smaller v (cid:107) . Con-sidering the relative alignment between the diffusive flux Figure 1.
The diffusive flux of resonant particles in velocityspace under the action of two arbitrary ( n and n ) resonances.The dark shaded areas represent isocontours of the VDFsof two particle populations. The red and dark-blue solidcurves show the diffusive trajectories, Eq. (28) with n = n and n = n . W n j and W n j represent the window functionsaccording to Eq. (20), in which the n and n resonancesare effective. The light-blue dashed semi-circles correspondto constant-energy contours. The black solid line indicates v (cid:107) = v g . and the constant-energy contour at point A, the diffusingparticles lose kinetic energy. This energy is transferredto the resonant wave, which consequently grows in am-plitude. Therefore, this situation corresponds to aninstability of the resonant wave. At point B, particles donot diffuse along the red solid curve because this pointlies outside W n j .At point C, resonant particles in W n j diffuse alongthe dark-blue solid curve toward greater v (cid:107) . Consideringthe relative alignment between the diffusive flux andthe constant-energy contour at point C, the diffusingparticles gain kinetic energy. This energy is taken fromthe resonant wave, which consequently shrinks in ampli-tude. Therefore, this situation corresponds to dampingof the resonant wave and counter-acts the driving ofthe instability through the n resonance. Because theresonant wave is unstable, the n resonant instabilitymust overcome the counter-acting n resonant damping.According to Eq. (18), there are three factors that de-termine the diffusion rate for the action of an n resonance.The first factor is the particle density f j within W nj . Thesecond factor is (cid:101) D nj whose magnitude is determined bythe polarization properties of the resonant waves. Thethird factor is the quantity ˆ G [ k n (cid:107) j ] f j /f j which defines therelative alignment between the isocontours of f i and thediffusive flux along the diffusion trajectory within W nj .In Fig. 1, the magnitude of | ˆ G [ k n (cid:107) j ] f j /f j | at point A isgreater than the magnitude of | ˆ G [ k n (cid:107) j ] f j /f j | at point C.Because the diffusive flux is directed from higher tolower values of f j , the quantity ˆ G [ k n (cid:107) j ] f j /f j resolves theambiguity in the directions of the trajectories for reso-nant particles. A careful analysis of ˆ G [ k n (cid:107) j ] using Eq. (29) shows that, if ( k (cid:107) / | k (cid:107) | )( ˆ G [ k n (cid:107) j ] f j /f j ) > at a givenresonant velocity, resonant particles diffuse toward asmaller value of v (cid:107) along the diffusive trajectory, while if ( k (cid:107) / | k (cid:107) | )( ˆ G [ k n (cid:107) j ] f j /f j ) < at a given resonant velocity,resonant particles diffuse toward a greater value of v (cid:107) .2.4. Numerical Analysis of the Quasi-Linear DiffusionEquation
To simulate the VDF evolution and to compare the dif-fusion rates between resonances quantitatively, a rigorousnumerical analysis of Eq. (10) is necessary. For this pur-pose, we develop a mathematical approach based on theCrank-Nicolson scheme and present the mathematicaldetails in Appendix A. Our approach is applicable to alltwo-dimensional diffusion equations with off-diagonal dif-fusion terms. Our numerical solution, given by Eq. (A28),evolves the VDF under the action of multiple resonancesin one time step. We tested our numerical solution byshowing that the diffusive flux obeys the predicted diffu-sion properties discussed in Section 2.3. FAST-MAGNETOSONIC/WHISTLER WAVE ANDELECTRON-STRAHL SCATTERINGAs an example, we apply our model developed in Sec-tion 2 to an electron resonant instability in the solarwind. The fast-magnetosonic/whistler (FM/W) wavepropagating in the anti-sunward direction and with anangle of ∼ ◦ with respect to the background magneticfield scatters the electron strahl (Vasko et al. 2019; Ver-scharen et al. 2019). Because this prediction is based onlinear theory, our quasi-linear framework is appropriatefor demonstrating the action of this instability on theelectron strahl.3.1. Linear Dispersion Relation
To find the characteristics of the unstable obliqueFM/W wave, we numerically solve the hot-plasma dis-persion relation with the NHDS code (Verscharen &Chandran 2018). We use the same plasma parameters asVerscharen et al. (2019), which are, notwithstanding thewide range of natural variation, representative for the av-erage electron parameters in the solar wind (Wilson et al.2019). We assume that the initial plasma consists ofisotropic Maxwellian protons, core electrons, and strahlelectrons. The subscripts p , e , c and s indicate protons,electrons, electron core, and electron strahl, respectively.We choose our coordinate system so that the anti-sunward and obliquely propagating FM/W waves have k (cid:107) > . We set β c = β p = 1 and β s = 0 . , where β j ≡ (8 πn j k B T j ) /B , n j and T j are the density andtemperature of species j , and k B is the Boltzmann con-stant. We set n p = n e , n c = 0 . n p , n s = 0 . n p , Figure 2.
NHDS solutions provide γ k (dashed curves, axis onthe left) and ω k (solid curves, axis on the right) as functionsof the k (cid:107) -component of the wavevector k . We show solutionsfor θ = 51 ◦ , θ = 55 ◦ , and θ = 59 ◦ . T c = T p , and T s = 2 T p . In the proton rest frame, weset n c U c + n s U s = 0 . We initialize the core and strahlbulk velocity with U c /v Ae = − . and U s /v Ae = 2 . where v Ae ≡ B / √ πn e m e is the electron Alfvén speed.NHDS finds that, under these plasma parameters, γ k > at angles between θ = 51 ◦ and θ = 67 ◦ . Our strahlbulk velocity then provides a maximum growth rate of γ k / | Ω e | = 10 − ( ? ).Fig. 2 shows γ k and ω k as functions of the k (cid:107) -component of the wavevector k for three different θ . Theoblique FM/W instability has its maximum growth rateat θ = 55 ◦ , while γ k > for . (cid:46) k (cid:107) v Ae / | Ω e | (cid:46) . which is the parallel unstable k -spectrum. As definedin Section 2.1, we acquire k (cid:107) v Ae / | Ω e | ≈ . . Thisvalue with Eqs. (7)-(9) leads to k ⊥ v Ae / | Ω e | = 0 . , ω k / | Ω e | ≈ . and v g /v Ae ≈ . . We also acquire σ (cid:107) v Ae / | Ω e | ≈ . and σ ⊥ v Ae / | Ω e | ≈ . from theunstable k -spectrum.3.2. Theoretical Description of the Quasi-LinearDiffusion in the FM/W Instability
Using the wave and plasma parameters from the pre-vious section, we describe the electron strahl and corediffusion in velocity space. In our analysis, we only con-sider the n = +1 , − and resonances, ignoring higher- n resonances due to their negligible contributions.Upon substituting our wave parameters into Eq. (20),we quantify the dimensionless window functions W ne v Ae with n = +1 , − and . In Fig. 3, the red, dark-blue andorange lines represent W +1 e v Ae , W − e v Ae and W e v Ae ,which are maximum at v +1 (cid:107) res /v Ae = 4 . , v − (cid:107) res /v Ae = − . and v (cid:107) res /v Ae = 0 . , respectively. We reiteratethat the superscripts indicate the n resonance. The black line indicates v (cid:107) = v g . Each W ne v Ae shows the v (cid:107) -range in which the quasi-linear diffusion through eachresonance is effective. We note that the W ne v Ae forthe three resonances have different widths in v (cid:107) -spaceand maximum values due to the different magnitudes of | v n (cid:107) res − v g | (see Eq. (24)). By substituting our waveparameters into Eq. (28), the diffusive trajectories forthe n = +1 , − and resonances are given by ( v ⊥ /v Ae ) + 1 . (cid:0) v (cid:107) /v Ae − . (cid:1) = const, (31) ( v ⊥ /v Ae ) + 0 . (cid:0) v (cid:107) /v Ae − . (cid:1) = const, (32)and ( v ⊥ /v Ae ) = const. (33)Eqs. (31) and (32) describe ellipses with their axes ori-ented along the v ⊥ - and v (cid:107) -directions. In Eq. (33), theperpendicular velocity of resonant particles is constant.Fig. 4 illustrates the electron diffusion from these threeresonances. We show the v (cid:107) -ranges in which the threeresonances are effective according to W +1 e v Ae , W − e v Ae and W e v Ae from Fig. 3. The red, dark-blue, and orangesolid curves represent the contours given by Eqs. (31)-(33), respectively. The light-blue dashed semi-circlescorrespond to constant-energy contours in the protonrest frame (see Eq. (30)). The black line indicates v (cid:107) = v g . For the initial strahl and core VDF, we apply theplasma parameters in Section 3.1 to the dimensionlessMaxwellian distribution: f Mj = n j v Ae π / n p v th,j exp (cid:34) − v ⊥ + ( v (cid:107) − U j ) v th,j (cid:35) , (34)where v th,j ≡ (cid:112) k B T j /m j . The red and blue areas inFig. 4 represent f Ms and f Mc which are normalized bythe maximum value of f Mc and plotted up to a value of − . In this normalization, Fig. 4 does not reflect therelative density ratio between both electron species.Due to the v (cid:107) -profile of W +1 e , the n = +1 resonance hasa significant effect on f Ms . As discussed in Section 2.3,because ( k (cid:107) / | k (cid:107) | )( ˆ G [ k +1 (cid:107) e ] f Ms /f Ms ) > , this resonanceleads to a diffusion of the resonant strahl electrons in W +1 e along trajectories represented by the red arrows.According to Eq. (30), the phase-space trajectory ofparticles that diffuse without a change in kinetic energyis described by (cid:18) dv ⊥ dv (cid:107) (cid:19) E = − v (cid:107) v ⊥ . (35)According to Eq. (31), the phase-space trajectory of res-onant particles fulfilling the n = +1 resonance, indicated Figure 3.
The red, dark-blue, and orange curves illustrate W +1 e v Ae , W − e v Ae and W e v Ae for the oblique FM/W wave.The black solid line represents v (cid:107) = v g . Each W ne v Ae showsthe v (cid:107) -range in which the corresponding resonance is effective.Each W ne v Ae has a different width in v (cid:107) -space and maxi-mum value due to a different magnitude of | v n (cid:107) res − v g | (seeEq. (24)). Figure 4.
The red, dark-blue, and orange arrows illustratethe diffusive flux for the n = +1 , − and resonances in theoblique FM/W instability. The red and dark-blue filled semi-circles represent isocontours of the strahl and core VDF. Thisfigure does not reflect the relative densities of both electronspecies. The light-blue dashed semi-circles correspond toconstant-energy contours. The black solid line indicates v (cid:107) = v g . by superscript +1 , is described by (cid:18) dv ⊥ dv (cid:107) (cid:19) +1 = − . v Ae v ⊥ (cid:18) v (cid:107) v Ae − . (cid:19) . (36)Comparing Eqs. (35) and (36) in W +1 e shows that | ( dv ⊥ /dv (cid:107) ) +1 | < | ( dv ⊥ /dv (cid:107) ) E | for the resonant electrons.Therefore, resolving the ambiguity in the directions ofthe trajectories, the distance of resonant strahl electronsfrom the origin of the coordinate system decreases. Thisdecrease in v ⊥ + v (cid:107) represents a loss of kinetic energyof the resonant strahl electrons. The n = +1 resonance, therefore, contributes to the driving of the FM/W insta-bility.Due to the v (cid:107) -profile of W − e , the n = − res-onance has a significant effect on f Mc . Because ( k (cid:107) / | k (cid:107) | )( ˆ G [ k − (cid:107) e ] f Mc /f Mc ) < , this resonance leads toa diffusion of the resonant core electrons in W − e alongtrajectories represented by the dark-blue arrows. Accord-ing to Eq. (32), the phase-space trajectory of resonantparticles fulfilling the n = − resonance, indicated bysuperscript − , is described by (cid:18) dv ⊥ dv (cid:107) (cid:19) − = − . v Ae v ⊥ (cid:18) v (cid:107) v Ae − . (cid:19) . (37)Comparing Eqs. (35) and (37) in W − e shows that | ( dv ⊥ /dv (cid:107) ) − | > | ( dv ⊥ /dv (cid:107) ) E | for the resonant electrons.Therefore, resolving the ambiguity in the directions ofthe trajectories, the distance of resonant core electronsfrom the origin of the coordinate system increases. Thisincrease in v ⊥ + v (cid:107) represents a gain of kinetic energyof the resonant core electrons. The n = − resonance,therefore, counter-acts the FM/W instability throughthe n = +1 resonance.Due to the v (cid:107) -profile of W e , the n = 0 resonancehas a significant effect on electrons in the v (cid:107) -rangein which f Mc > f Ms and ∂f Mc /∂v (cid:107) < . Because ( k (cid:107) / | k (cid:107) | )( ˆ G [ k (cid:107) e ] f Mc /f Mc ) < , the resonant electrons in W e diffuse along trajectories represented by the yellowarrows. Because the distance of these electrons fromthe origin of the coordinate system increases, these res-onant electrons diffuse toward greater kinetic energies.This diffusion removes energy from the resonant FM/Wwaves and thus counter-acts the driving of the FM/Winstability through the n = +1 resonance.Fig. 4 only illustrates the nature of the quasi-lineardiffusion through the n = +1 , − and resonances in ve-locity space. It does not give any information regardingthe relative strengths of the diffusion rates between thethree resonances. Because the FM/W wave is unstableaccording to linear theory, the n = +1 resonant instabil-ity must dominate over any counter-acting contributionsfrom the n = − and resonances though.3.3. Numerical Description of the Quasi-LinearDiffusion in the FM/W Instability
We use our numerical procedure from Eq. (A28) tosimulate the quasi-linear diffusion through the n = +1 , − and resonances, predicted in Section 3.2. Accordingto the definitions in Appendix A, we select the discretiza-tion parameters N v = 60 , v ⊥ max /v Ae = v (cid:107) max /v Ae = 7 and | Ω e | ∆ t = 1 . For the computation of Eq. (A28), weuse the same parameters of resonant FM/W waves asthose presented in Section 3.1 and quantify (cid:101) D ± e and (cid:101) D e in Eq. (19).In (cid:101) D ne for each resonance, we only consider the J term in ψ nj , ignoring higher-order Bessel functions dueto their small contributions. Our NHDS solutions showthat | E y | ≈ . | E x | and | E z | ≈ . | E x | in the unsta-ble k -spectrum. Then, we set | E R | ≈ | E L | ≈ . | E x | .Faraday’s law yields E x ≈ [ ω k / ( k (cid:107) c )] B y when ignoringthe small contributions from any E z terms. This allowsus to express E x through B y in ψ nj , where B y representsthe peak amplitude of the wave magnetic-field fluctu-ations. For simplicity, we assume that B y is constantin time during the quasi-linear diffusion. Under theseassumptions, we acquire (cid:101) D ± e ≈ W ± e . π | Ω e | v ⊥ σ (cid:107) σ ⊥ (cid:20) B y B ω k k (cid:107) (cid:21) × (cid:90) ∞ J ( ρ e ) exp (cid:20) − ( k ⊥ − k ⊥ ) σ ⊥ (cid:21) k ⊥ dk ⊥ , (38)and (cid:101) D e ≈ W e . π | Ω e | v (cid:107) σ (cid:107) σ ⊥ (cid:20) B y B ω k k (cid:107) (cid:21) × (cid:90) ∞ J ( ρ e ) exp (cid:20) − ( k ⊥ − k ⊥ ) σ ⊥ (cid:21) k ⊥ dk ⊥ , (39)where the relative amplitude B y /B is a free parameter,and we set B y /B = 0 . . Then, we apply Eqs. (38)and (39) to Eq. (A28).We initialize our numerical computation with the same f Ms and f Mc as defined in Section 3.2. Fig. 5a representsthe normalized f e = f Mc + f Ms , plotted up to a value of − . Fig. 5b shows f e evolved through the n = +1 , − and resonances, resulting from our iterative calculationof Eq. (A28). Considering the maximum value of theinstability’s growth rate, γ k / | Ω e | = 4 . × − in Fig. 2,we terminate the evaluation of our numerical computa-tion at | Ω e | t = 5 × which corresponds to γ k t ∼ andthus a reasonable total growth of the unstable FM/Wwaves.The strahl electrons at around v (cid:107) /v Ae ≈ . diffusethrough the n = +1 resonance, as theoretically predictedin Fig. 4. This diffusion increases the pitch-angle of theresonant strahl electrons and generates a strong pitch-angle gradient at v (cid:107) /v Ae ≈ . . During this process, the v ⊥ of the scattered strahl electrons increases while their v (cid:107) decreases.Because the longitudinal component of the electric-fieldfluctuations is weaker than their transverse components,the diffusion through the n = 0 resonance is only slightlynoticeable over the modeled time interval. The diffusionthrough the n = − resonance is not noticeable even Figure 5.
Panel (a): the initial electron VDF; Panel (b):the electron VDF evolved through the n = +1 , − and resonances. Compared to Fig. 4, the effect of the n = +1 resonance dominates the evolution during the time γ k t ∼ .It causes a significant pitch-angle gradient at v (cid:107) /v Ae ≈ . through the scattering of strahl electrons. An animationof this figure is available. The animation shows the timeevolution of the distribution function from | Ω e | t = 0 to | Ω e | t = 5 × . During this evolution, the strahl scatteringtoward larger v ⊥ is visible. though (cid:101) D − e and (cid:101) D +1 e have similar magnitudes. Thisis because | ˆ G [ k − (cid:107) e ] f e /f e | in W − e is much smaller than | ˆ G [ k +1 (cid:107) e ] f e /f e | in W +1 e , as discussed in Section 2.3, andthe number of core electrons in W − e is very small (seeFig. 4 and 5). THE SECONDARY EFFECT OF COULOMBCOLLISIONSBecause the collisionless action of resonant wave-particle instabilities often forms strong pitch-angle gradi-ents (see, for example, Fig. 5), collisions can be enhancedin the plasma. Therefore, a more realistic evolution ofthe total electron VDF must account for the action ofCoulomb collisions of strahl electrons with core electronsand protons. For this purpose, we adopt the Fokker-Planck equation given by Ljepojevic et al. (1990) withRosenbluth potentials (Rosenbluth et al. 1957) and nor-0malize it in our dimensionless system of units as (cid:18) ∂f j ∂t (cid:19) F okker − P lanck = (cid:88) b Γ jb (cid:26) π m j m b f b f j + ∂h∂v α ∂f j ∂v α + 12 ∂ g∂v α ∂v β ∂ f j ∂v α ∂v β (cid:27) , (40)where g ( v ) ≡ (cid:90) f b ( v (cid:48) ) | v − v (cid:48) | d v (cid:48) , (41) h ( v ) ≡ m b − m j m b (cid:90) f b ( v (cid:48) ) | v − v (cid:48) | − d v (cid:48) , (42)and Γ jb ≡ πn b v Ae | Ω e | (cid:18) Z j Z b q j m j (cid:19) ln Λ jb . (43)The subscript b indicates the species of background parti-cles, with which the particles of species j Coulomb-collide.The quantity ln Λ jb is the Coulomb logarithm and typi-cally ln Λ jb ≈ in space plasmas. The parameters Z j and Z b are the atomic masses of a particle of species j and b . The superscripts α and β indicate the compo-nent of the velocity in cylindrical coordinates and thesummation convention holds.We assume that the timescale of Coulomb collisionsis much longer than the timescale of the quasi-lineardiffusion in the solar wind under our set of parameters.This assumption allows us to model the resonant wave-particle instability first and to use the resulting VDFas the input for the model of the subsequent, secondaryeffects of the collisions.Based on our mathematical approach presented inAppendix A, we present our numerical scheme to solvethe Fokker-Planck equation, Eq. (40), in Appendix B.We tested our numerical solutions, Eq. (B31), by showingthat a set of arbitrary test VDFs diffuses toward f b withtime.For the computation of Eq. (B31), we set isotropicMaxwellian electron-core and proton VDFs as back-ground species, f b = f Mc and f b = f Mp , for which weapply the plasma parameters presented in Section 3.1to Eq. (34). In this numerical computation, we selectthe discretization parameters N v = 60 , v ⊥ max /v Ae = v (cid:107) max /v Ae = 7 and | Ω e | ∆ t = 10 . Moreover, we set B = 5 × − G and n b = 10 cm − in Eq. (43), whichare representative for the conditions in the solar wind ata distance of 0.3 au from the Sun. We initialize f j withthe electron-strahl VDF f s from our quasi-linear analysisof the oblique FM/W instability at time | Ω e | t = 5 × .In this setup, our initial electron VDF for the Coulombcollision analysis is same as the electron VDF shown inFig. 5b. Figure 6.
Panel (a): the electron VDF as initial condi-tion for our collision analysis; Panel (b): the electron VDFevolved through Coulomb collisions of strahl electrons withcore electrons and protons. The strong pitch-angle gradientat v (cid:107) /v Ae ≈ . (shown in Fig. 6a and Fig. 5b) is relaxedthrough Coulomb collisions. However, the required timefor a noticeable collisional effect on that gradient is around times longer than the timescale of the strahl scatter-ing. An animation of this figure is available. The animationshows the time evolution of the distribution function from | Ω e | t = 5 × to | Ω e | t = 7 × . During this evolution, thecollisional smoothing of the pitch-angle gradients is visible. The iterative calculation of Eq. (B31) results in thetime evolution of the electron-strahl VDF under the ac-tion of Coulomb collisions with core electrons and protons.The result of this computation at the time | Ω e | t = 7 × is shown in Fig. 6. A detailed comparison of the dis-tribution function before (Fig. 6a) and after (Fig. 6b)our calculation of the effect of Coulomb collisions revealsthat Coulomb collisions relax the strong pitch-angle gra-dient at v (cid:107) /v Ae ≈ . , which resulted from the action ofthe oblique FM/W instability. However, the Coulombcollisions are only capable of affecting strong pitch-anglegradients in the modified electron VDF under our plasmaparameters. In addition, the required time for a notice-able collisional effect on this pitch-angle gradient is oforder times longer than the characteristic timescaleof the quasi-linear diffusion.1 DISCUSSION OF THE STRAHL SCATTERINGThe numerical computation of Eq. (17) shows that dH/dt is negative and asymptotically tends toward zeroas the electron VDF evolves through the oblique FM/Winstability and the counter-acting damping effects untilthe time | Ω e | t = 5 × , which is presented in Fig. 5.Therefore, our quasi-linear diffusion model reflects thestabilization of the particle VDF through the participat-ing wave-particle resonances.During the action of the oblique FM/W instability,the scattered strahl electrons reduce their collimationalong the B -direction and become more isotropic. Eventhough this instability does not cause significant strahlscattering, we argue that it contributes to the initialformation of the halo population. However, other mech-anisms must be considered to account for the full strahlscattering, in agreement with observations (Gurgioloet al. 2012; Gurgiolo & Goldstein 2016).Alternative models describing Coulomb-collisional ef-fects on the strahl VDF suggest that an anomalous-diffusion process must be considered in order to achievean agreement with observations (Lemons & Feldman1983; Horaites et al. 2018, 2019). We note that our anal-ysis includes the subsequent action of Coulomb collisionsafter the action of collisionless wave-particle resonancesassuming plasma parameters consistent with the solarwind at a distance of about 0.3 au from the Sun. Ourcollisional effects are similar to those proposed by Vockset al. (2005). However, our model predicts that the colli-sional relaxation is so subtle that the strahl scatteringthrough collisions is barely noticeable for the analyzedphase of the VDF evolution.The clear separation of timescales between wave-particle effects and Coulomb-collisional effects com-plicates the description of the VDF evolution on helio-spheric scales, because other processes act on comparabletimescales. These additional processes, which our analy-sis ignores, include turbulence, shocks, plasma mixing,plasma expansion, and magnetic focusing (Feldman et al.1983; Fitzenreiter et al. 2003; Ryu et al. 2007; Yoonet al. 2012; Tang et al. 2020). A complete model for theradial evolution of the VDF must quantify and accountfor these processes as well. In the context of our work,these processes can potentially push a VDF that hasundergone stabilization as shown in Fig. 5b into theunstable regime again. In this case, dH/dt in Eq. (17)returns to a non-zero value, which signifies a new onsetof wave-particle resonances and further scattering ofresonant particles. CONCLUSIONSWave-particle resonances are important plasma-physicsprocesses in many astrophysical plasmas. Often, fullynon-linear simulations with codes solving the equationsof kinetic plasma theory are used to model the evolu-tion of the distribution function under the action ofwave-particle resonances. However, quasi-linear theoryaugments this approach as it allows us to study thecontributions of different processes to these resonances.Therefore, quasi-linear theory is a very helpful tool toimprove our understanding of wave-particle resonancesin astrophysical plasmas.We propose a quasi-linear diffusion model for any gen-eralized wave-particle instability. We analyze the quasi-linear diffusion equation by expressing the electric fieldof an arbitrary unstable and resonant wave mode as aGaussian wave packet. From Boltzmann’s H -theorem inour quasi-linear analysis, we define a window functionthat determines the specific velocity-space range in whicha dominant wave-particle instability and counter-actingdamping contributions are effective. This window func-tion is the consequence of the localized energy densityof our Gaussian wave packet both in configuration spaceand in wavevector space.Moreover, we derive a relation describing the diffu-sive trajectories of the resonant particles for such aninstability in velocity space. These trajectories evolvethe particle VDF into a stable state in which no furtherquasi-linear diffusion occurs. Therefore, our theoreticalmodel illustrates the diffusion and stabilization whichresonant particles, depending on their location in velocityspace, experience in wave-particle resonances.For the computational quantification of our theoreticalmodel, we introduce a mathematical approach basedon the Crank-Nicolson scheme to numerically solve thefull quasi-linear diffusion equation. We highlight thatthis mathematical approach applies to all general two-dimensional diffusion equations, including those withoff-diagonal diffusion terms.As an example, we apply our model to the obliqueFM/W instability that scatters strahl electrons in thesolar wind. Our model shows that the n = +1 resonantinstability of FM/W waves propagating with an angleof ∼ ◦ with respect to the background magnetic fieldscatters strahl electrons toward larger v ⊥ and smaller v (cid:107) . The strahl scattering instability through the n = +1 resonance dominates over the counter-acting dampingcontributions through the n = − and n = 0 resonances.This instability creates a strong pitch-angle gradient inthe electron-strahl VDF.By numerically solving the Fokker-Planck equation,we show that Coulomb collisions of strahl electrons with2core electrons and protons relax this strong pitch-anglegradient on a timescale about times longer thanthe timescale of the collisionless strahl scattering. Thisfinding suggests that collisional effects are negligible inthe strahl-driven oblique FM/W instability, which isa representative example for a resonant wave-particleinstability in the solar wind.Our predicted evolution of the electron VDF is consis-tent with the observed formation of a proto-halo throughstrahl scattering (Gurgiolo et al. 2012). However, furtherobservations are ambiguous regarding the exact sourceof the proto-halo (Gurgiolo & Goldstein 2016). Futurehigh-resolution electron observations with Solar Orbiterand Parker Solar Probe at different distances from theSun may help us resolve these ambiguities.Our general quasi-linear diffusion model applies toall non-relativistic collisionless plasmas, such as plan- etary magnetospheres (e.g. Mourenas et al. 2015). Italso applies to other types of wave-particle instabilitiesin plasmas such as the resonant instabilities driven bytemperature anisotropy or by relative drift. We espe-cially note that our model is also capable of describingion-driven instabilities.ACKNOWLEDGMENTSWe appreciate helpful discussions with Georgios Nico-laou, Konstantinos Horaites, and Jung Joon Seough. D.V.is supported by the STFC Ernest Rutherford FellowshipST/P003826/1. D.V., R.T.W., and and A.N.F. are sup-ported by STFC Consolidated Grant ST/S000240/1.APPENDIX A. NUMERICAL ANALYSIS OF THE QUASI-LINEAR DIFFUSION EQUATIONEq. (10) is a second-order differential equation that includes cross-derivative operators such as ∂ /∂v (cid:107) ∂v ⊥ . In orderto simultaneously evaluate the ∂ /∂v (cid:107) ∂v ⊥ operators with the ∂ /∂v (cid:107) and ∂ /∂v ⊥ operators in Eq. (10), we dividevelocity space into N v × N v steps with equal step sizes of ∆ v/ by defining the outer boundaries of velocity spaceas ± v ⊥ max and ± v (cid:107) max . The v ⊥ -index M and the v (cid:107) -index N both step through 1, 3/2, 2, . . . , N v , N v + 1 / . Wedefine the discrete velocity coordinates as v ⊥ M ≡ − v ⊥ max + ( M − v and v (cid:107) N ≡ − v (cid:107) max + ( N − v . We notethat this definition introduces negative v ⊥ -values that, although they simplify our numerical analysis, we ignore in ourcomputational results. We divide the time t with equal step sizes of ∆ t and the t -index T steps through , , , · · · . Wedefine the discrete time as t T ≡ ( T − t . We then define the discrete VDF as f TM,N ≡ f j ( v ⊥ M , v (cid:107) N , t T ) . For thediscretization of the velocity derivatives, we adopt the two-point central difference operators (Gilat & Subramaniam2011) ∂f j ( v ⊥ M , v (cid:107) N , t T ) ∂v ⊥ ≈ f TM +1 / ,N − f TM − / ,N ∆ v , (A1)and ∂f j ( v ⊥ M , v (cid:107) N , t T ) ∂v (cid:107) ≈ f TM,N +1 / − f TM,N − / ∆ v . (A2)For the discretization of the time derivative, we adopt the forward difference operator ∂f j ( v ⊥ M , v (cid:107) N , t T ) ∂t ≈ f T +1 M,N − f TM,N ∆ t . (A3)By using Eqs. (A1) and (A2), we discretize the right-hand side of Eq. (10) and express it as ( ∂f /∂t ) TM,N (cid:18) ∂f∂t (cid:19)
TM,N ≡ ∞ (cid:88) n = −∞ (cid:90) (cid:34)(cid:32) − k (cid:107) v (cid:107) N ω k + v g ( k (cid:107) − k (cid:107) ) (cid:33) v ⊥ M D nM +1 / ,N [ ˆ Gf ] TM +1 / ,N − D nM − / ,N [ ˆ Gf ] TM − / ,N ∆ v + k (cid:107) ω k + v g ( k (cid:107) − k (cid:107) ) D nM,N +1 / [ ˆ Gf ] TM,N +1 / − D nM,N − / [ ˆ Gf ] TM,N − / ∆ v (cid:35) d k , (A4)where [ ˆ Gf ] TM,N ≡ (cid:32) − k (cid:107) v (cid:107) N ω k + v g ( k (cid:107) − k (cid:107) ) (cid:33) v ⊥ M f TM +1 / ,N − f TM − / ,N ∆ v + k (cid:107) ω k + v g ( k (cid:107) − k (cid:107) ) f TM,N +1 / − f TM,N − / ∆ v , (A5)3and D nM,N ≡ D nj (cid:12)(cid:12) v ⊥ = v ⊥ M v (cid:107) = v (cid:107) N . (A6)According to the Crank-Nicolson scheme (Iserles 2008), the full discretization of Eq. (10) in its time and velocityderivatives is then given by f T +1 M,N − ∆ t (cid:18) ∂f∂t (cid:19) T +1 M,N = f TM,N + ∆ t (cid:18) ∂f∂t (cid:19) TM,N . (A7)By using Eqs. (A4)-(A6) and resolving the δ -functions in D nM ± / ,N and D nM,N ± / through the k (cid:107) -integral, we rewriteEq. (A7) as f T +1 M,N − ∞ (cid:88) n = −∞ (cid:40) µ P nM,N (cid:101) D nM +1 / ,N (cid:104) P nM +1 / ,N (cid:16) f T +1 M +1 ,N − f T +1 M,N (cid:17) + Q nN (cid:16) f T +1 M +1 / ,N +1 / − f T +1 M +1 / ,N − / (cid:17)(cid:105) − µ P nM,N (cid:101) D nM − / ,N (cid:104) P nM − / ,N (cid:16) f T +1 M,N − f T +1 M − ,N (cid:17) + Q nN (cid:16) f T +1 M − / ,N +1 / − f T +1 M − / ,N − / (cid:17)(cid:105) + µ Q nN +1 / (cid:101) D nM,N +1 / (cid:104) P nM,N +1 / (cid:16) f T +1 M +1 / ,N +1 / − f T +1 M − / ,N +1 / (cid:17) + Q nN +1 / (cid:16) f T +1 M,N +1 − f T +1 M,N (cid:17)(cid:105) − µ Q nN − / (cid:101) D nM,N − / (cid:104) P nM,N − / (cid:16) f T +1 M +1 / ,N − / − f T +1 M − / ,N − / (cid:17) + Q nN − / (cid:16) f T +1 M,N − f T +1 M,N − (cid:17)(cid:105) (cid:41) = f TM,N + ∞ (cid:88) n = −∞ (cid:40) µ P nM,N (cid:101) D nM +1 / ,N (cid:104) P nM +1 / ,N (cid:0) f TM +1 ,N − f TM,N (cid:1) + Q nN (cid:16) f TM +1 / ,N +1 / − f TM +1 / ,N − / (cid:17)(cid:105) − µ P nM,N (cid:101) D nM − / ,N (cid:104) P nM − / ,N (cid:0) f TM,N − f TM − ,N (cid:1) + Q nN (cid:16) f TM − / ,N +1 / − f TM − / ,N − / (cid:17)(cid:105) + µ Q nN +1 / (cid:101) D nM,N +1 / (cid:104) P nM,N +1 / (cid:16) f TM +1 / ,N +1 / − f TM − / ,N +1 / (cid:17) + Q nN +1 / (cid:0) f TM,N +1 − f TM,N (cid:1)(cid:105) − µ Q nN − / (cid:101) D nM,N − / (cid:104) P nM,N − / (cid:16) f TM +1 / ,N − / − f TM − / ,N − / (cid:17) + Q nN − / (cid:0) f TM,N − f TM,N − (cid:1)(cid:105) (cid:41) , (A8)where P nM,N ≡ n Ω j [ v (cid:107) N − v g ][( ω k − k (cid:107) v g ) v (cid:107) N − n Ω j v g ] v ⊥ M , (A9) Q nN ≡ ω k − k (cid:107) v g − n Ω j ( ω k − k (cid:107) v g ) v (cid:107) N − n Ω j v g , (A10) (cid:101) D nM,N ≡ (cid:101) D nj (cid:12)(cid:12) v ⊥ = v ⊥ M v (cid:107) = v (cid:107) N , (A11)and µ ≡ ∆ t/ (∆ v ) . Eq. (A8) is a two-dimensional set of algebraic equations, the solution of which, f T +1 M,N , for all v ⊥ -and v (cid:107) -indexes describes the VDF at time T + 1 based on f TM,N for all v ⊥ - and v (cid:107) -indexes.In order to transform Eq. (A8) into a single matrix equation with a tridiagonal matrix, we introduce the concept of adouble matrix. On both sides of Eq. (A8), we group the terms by the same v ⊥ -index in the VDF and rearrange thesegroups in increasing order in v ⊥ -index. In each group, we then rearrange terms in increasing order in their v (cid:107) -index inthe VDF. Then, we have − η ( µ ) (1) M,N f T +1 M − ,N − ξ ( µ ) (2) M,N f T +1 M − / ,N − / + ξ ( µ ) (1) M,N f T +1 M − / ,N +1 / − α ( µ ) (2) M,N f T +1 M,N − + α ( µ ) M,N f T +1 M,N − α ( µ ) (1) M,N f T +1 M,N +1 + ξ ( µ ) (4) M,N f T +1 M +1 / ,N − / − ξ ( µ ) (3) M,N f T +1 M +1 / ,N +1 / − η ( µ ) (2) M,N f T +1 M +1 ,N = − η ( − µ ) (1) M,N f TM − ,N − ξ ( − µ ) (2) M,N f TM − / ,N − / + ξ ( − µ ) (1) M,N f TM − / ,N +1 / − α ( − µ ) (2) M,N f TM,N − + α ( − µ ) M,N f TM,N − α ( − µ ) (1) M,N f TM,N +1 + ξ ( − µ ) (4) M,N f TM +1 / ,N − / − ξ ( − µ ) (3) M,N f TM +1 / ,N +1 / − η ( − µ ) (2) M,N f TM +1 ,N , (A12)4where α ( µ ) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:16) P nM,N P nM +1 / ,N (cid:17) (cid:101) D nM +1 / ,N + (cid:16) P nM,N P nM − / ,N (cid:17) (cid:101) D nM − / ,N + (cid:16) Q nN +1 / (cid:17) (cid:101) D nM,N +1 / + (cid:16) Q nN − / (cid:17) (cid:101) D nM,N − / (cid:21) , (A13) α ( µ ) (1) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:16) Q nN +1 / (cid:17) (cid:101) D nM,N +1 / (cid:21) , (A14) α ( µ ) (2) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:16) Q nN − / (cid:17) (cid:101) D nM,N − / (cid:21) , (A15) ξ ( µ ) (1) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:0) P nM,N Q nN (cid:1) (cid:101) D nM − / ,N + (cid:16) P nM,N +1 / Q nN +1 / (cid:17) (cid:101) D nM,N +1 / (cid:21) , (A16) ξ ( µ ) (2) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:0) P nM,N Q nN (cid:1) (cid:101) D nM − / ,N + (cid:16) P nM,N − / Q nN − / (cid:17) (cid:101) D nM,N − / (cid:21) , (A17) ξ ( µ ) (3) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:0) P nM,N Q nN (cid:1) (cid:101) D nM +1 / ,N + (cid:16) P nM,N +1 / Q nN +1 / (cid:17) (cid:101) D nM,N +1 / (cid:21) , (A18) ξ ( µ ) (4) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:0) P nM,N Q nN (cid:1) (cid:101) D nM +1 / ,N + (cid:16) P nM,N − / Q nN − / (cid:17) (cid:101) D nM,N − / (cid:21) , (A19) η ( µ ) (1) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:16) P nM,N P nM − / ,N (cid:17) (cid:101) D nM − / ,N (cid:21) , (A20)and η ( µ ) (2) M,N ≡ µ ∞ (cid:88) n = −∞ (cid:20) (cid:16) P nM,N P nM +1 / ,N (cid:17) (cid:101) D nM +1 / ,N (cid:21) . (A21)All terms on both sides of Eq. (A12) with a constant v ⊥ -index account for variations in v (cid:107) -space only. Therefore, theycan be grouped into a single system of one-dimensional algebraic equations.We transform all terms with v ⊥ -index M on both sides of Eq. (A12) into the tridiagonal matrices [ A ( µ ) M ][ F T +1 M ] and [ A ( − µ ) M ][ F TM ] , where F TM ≡ [ f TM, f TM, f TM, · · · f TM,N v f TM,N v + ] T × N v ( T represents the transpose of a matrix), and A ( µ ) M ≡ α ( µ ) M, − α ( µ ) (1) M, · · · α ( µ ) M, / − α ( µ ) (1) M, / · · · − α ( µ ) (2) M, α ( µ ) M, − α ( µ ) (1) M, · · · ... . . . ... · · · − α ( µ ) (2) M,N v α ( µ ) M,N v · · · − α ( µ ) (2) M,N v +1 / α ( µ ) M,N v +1 / N v × N v . (A22)We transform all terms with v ⊥ -index M − / on both sides of Eq. (A12) into the tridiagonal matrices [ B ( µ ) (1) M ][ F T +1 M − / ] and [ B ( − µ ) (1) M ][ F TM − / ] , where B ( µ ) (1) M ≡ ξ ( µ ) (1) M, · · · − ξ ( µ ) (2) M, / ξ ( µ ) (1) M, / · · · ... . . . ... · · · − ξ ( µ ) (2) M,N v ξ ( µ ) (1) M,N v · · · − ξ ( µ ) (2) M,N v +1 / N v × N v . (A23)5We transform all terms with v ⊥ -index M +1 / on both sides of Eq. (A12) into the tridiagonal matrices [ B ( µ ) (2) M ][ F T +1 M +1 / ] and [ B ( − µ ) (2) M ][ F TM +1 / ] , where B ( µ ) (2) M ≡ − ξ ( µ ) (3) M, · · · ξ ( µ ) (4) M, / − ξ ( µ ) (3) M, / · · · ... . . . ... · · · ξ ( µ ) (4) M,N v − ξ ( µ ) (3) M,N v · · · ξ ( µ ) (4) M,N v +1 / N v × N v . (A24)We transform all terms with v ⊥ -index M − on both sides of Eq. (A12) into the tridiagonal matrices [ C ( µ ) (1) M ][ F T +1 M − ] and [ C ( − µ ) (1) M ][ F TM − ] , where C ( µ ) (1) M ≡ − η ( µ ) (1) M, · · · − η ( µ ) (1) M, / · · · ... . . . ... · · · − η ( µ ) (1) M,N v +1 / N v × N v . (A25)Lastly, we transform all terms with v ⊥ -index M + 1 on both sides of Eq. (A12) into the tridiagonal matrices [ C ( µ ) (2) M ][ F T +1 M +1 ] and [ C ( − µ ) (2) M ][ F TM +1 ] , where C ( µ ) (2) M ≡ − η ( µ ) (2) M, · · · − η ( µ ) (2) M, / · · · ... . . . ... · · · − η ( µ ) (2) M,N v +1 / N v × N v . (A26)This strategy allows us to express Eq. (A12) as a single system of one-dimensional algebraic equations: [ C ( µ ) (1) M ][ F T +1 M − ] + [ B ( µ ) (1) M ][ F T +1 M − / ] + [ A ( µ ) M ][ F T +1 M ] + [ B ( µ ) (2) M ][ F T +1 M +1 / ] + [ C ( µ ) (2) M ][ F T +1 M +1 ]= [ C ( − µ ) (1) M ][ F TM − ] + [ B ( − µ ) (1) M ][ F TM − / ] + [ A ( − µ ) M ][ F TM ] + [ B ( − µ ) (2) M ][ F TM +1 / ] + [ C ( − µ ) (2) M ][ F TM +1 ] . (A27)Eq. (A27) only describes the VDF evolution in v ⊥ -space. However, each matrix term itself includes the VDF evolutionin v (cid:107) -space. We transform Eq. (A27) into a single tridiagonal matrix: E ( µ ) QLD F T +11 F T +13 / F T +12 ... F T +1 N v +1 / (2 N v ) × = E ( − µ ) QLD F T F T / F T ... F TN v +1 / (2 N v ) × , (A28)where E ( µ ) QLD ≡ A ( µ ) B ( µ ) (2)1 C ( µ ) (2)1 · · · B ( µ ) (1)3 / A ( µ ) / B ( µ ) (2)3 / C ( µ ) (2)3 / · · · C ( µ ) (1)2 B ( µ ) (1)2 A ( µ ) B ( µ ) (2)2 C ( µ ) (2)2 · · · ... . . . ... · · · C ( µ ) (1) N v B ( µ ) (1) N v A ( µ ) N v B ( µ ) (2) N v · · · C ( µ ) (1) N v +1 / B ( µ ) (1) N v +1 / A ( µ ) N v +1 / (2 N v ) × (2 N v ) . (A29)6Eq. (A28) is in the form of a double matrix, and E ( µ ) QLD in Eq. (A29) defines the evolution matrix. The innermatrices of E ( µ ) QLD evolve f TM,N in v (cid:107) -space while the outer matrices of E ( µ ) QLD evolve f TM,N in v ⊥ -space duringeach time step. By multiplying Eq. (A28) with the inverse of E ( µ ) QLD on both sides, Eq. (A28) provides the timeevolution of f TM,N in one time step simultaneously in the v ⊥ - and v (cid:107) -spaces. Therefore, it represents the numericalsolution of Eq. (10), which describes the quasi-linear diffusion of a VDF through all resonances. B. NUMERICAL ANALYSIS OF THE FOKKER-PLANCK EQUATIONIn this appendix, we present our numerical strategy to solve the Fokker-Planck equation for Coulomb collisions inEq. (40). Using the Crank-Nicolson scheme presented in Appendix A, we discretize Eq. (40) as f T +1 M,N − (cid:88) b Γ jb (cid:20) π (∆ v ) µ m j m b f b ( v ⊥ M , v (cid:107) N ) f T +1 M,N + µg (cid:107)⊥ M,N (cid:16) f T +1 M +1 / ,N +1 / − f T +1 M − / ,N +1 / − f T +1 M +1 / ,N − / + f T +1 M − / ,N − / (cid:17) + µg ⊥⊥ M,N (cid:16) f T +1 M +1 ,N − f T +1 M,N + f T +1 M − ,N (cid:17) + µg (cid:107)(cid:107) M,N (cid:16) f T +1 M,N +1 − f T +1 M,N + f T +1 M,N − (cid:17) + (∆ v ) µh ⊥ M,N (cid:16) f T +1 M +1 / ,N − f T +1 M − / ,N (cid:17) + (∆ v ) µh (cid:107) M,N (cid:16) f T +1 M,N +1 / − f T +1 M,N − / (cid:17)(cid:21) = f TM,N + (cid:88) b Γ jb (cid:20) π (∆ v ) µ m j m b f b ( v ⊥ M , v (cid:107) N ) f TM,N + µg (cid:107)⊥ M,N (cid:16) f TM +1 / ,N +1 / − f TM − / ,N +1 / − f TM +1 / ,N − / + f TM − / ,N − / (cid:17) + µg ⊥⊥ M,N (cid:16) f TM +1 ,N − f TM,N + f TM − ,N (cid:17) + µg (cid:107)(cid:107) M,N (cid:16) f TM,N +1 − f TM,N + f TM,N − (cid:17) + (∆ v ) µh ⊥ M,N (cid:16) f TM +1 / ,N − f TM − / ,N (cid:17) + (∆ v ) µh (cid:107) M,N (cid:16) f TM,N +1 / − f TM,N − / (cid:17)(cid:21) , (B30)where g ⊥⊥ M,N ≡ ∂ g/∂v ⊥ , g (cid:107)(cid:107) M,N ≡ ∂ g/∂v (cid:107) , g (cid:107)⊥ M,N ≡ ∂ g/∂v (cid:107) ∂v ⊥ , h ⊥ M,N ≡ ∂h/∂v ⊥ and h (cid:107) M,N ≡ ∂h/∂v (cid:107) , estimated at v ⊥ = v ⊥ M and v (cid:107) = v (cid:107) N .Eq. (B30) represents a system of two-dimensional algebraic equations. Therefore, we transform Eq. (B30) into asingle tridiagonal matrix using the same strategy for a double matrix as presented in Appendix A. E ( µ ) F F T +11 F T +13 / F T +12 ... F T +1 N v +1 / (2 N v ) × = E ( − µ ) F F T F T / F T ... F TN v +1 / (2 N v ) × , (B31)where E ( µ ) F ≡ X ( µ ) − Y ( µ ) Z ( µ ) · · · Y ( µ ) / X ( µ ) / − Y ( µ ) / Z ( µ ) / · · · Z ( µ ) Y ( µ ) X ( µ ) − Y ( µ ) Z ( µ ) · · · ... . . . ... · · · Z ( µ ) N v Y ( µ ) N v X ( µ ) N v − Y ( µ ) N v · · · Z ( µ ) N v +1 / Y ( µ ) N v +1 / X ( µ ) N v +1 / (2 N v ) × (2 N v ) , (B32)7 X ( µ ) M ≡ ε ( µ ) M, − ε ( µ ) (2) M, − ε ( µ ) (1) M, · · · ε ( µ ) (2) M, / ε ( µ ) M, / − ε ( µ ) (2) M, / − ε ( µ ) (1) M, / · · · − ε ( µ ) (1) M, ε ( µ ) (2) M, ε ( µ ) M, − ε ( µ ) (2) M, − ε ( µ ) (1) M, · · · ... . . . ... · · · − ε ( µ ) (1) M,N v ε ( µ ) (2) M,N v ε ( µ ) M,N v − ε ( µ ) (2) M,N v · · · − ε ( µ ) (1) M,N v +1 / ε ( µ ) (2) M,N v +1 / ε ( µ ) M,N v +1 / N v × N v , (B33) Y ( µ ) M ≡ (cid:37) ( µ ) (2) M, (cid:37) ( µ ) (1) M, · · · − (cid:37) ( µ ) (1) M, / (cid:37) ( µ ) (2) M, / (cid:37) ( µ ) (1) M, / · · · ... . . . ... · · · − (cid:37) ( µ ) (1) M,N v (cid:37) ( µ ) (2) M,N v (cid:37) ( µ ) (1) M,N v · · · − (cid:37) ( µ ) (1) M,N v +1 / (cid:37) ( µ ) (2) M,N v +1 / N v × N v , (B34) Z ( µ ) M ≡ − τ ( µ ) M, · · · − τ ( µ ) M, / · · · ... . . . ... · · · − τ ( µ ) M,N v +1 / N v × N v , (B35) ε ( µ ) M,N ≡ − µ (cid:88) b Γ jb (cid:20) π (∆ v ) m j m b f b ( v ⊥ M , v (cid:107) N ) − g ⊥⊥ M,N − g (cid:107)(cid:107) M,N (cid:21) , (B36) ε ( µ ) (1) M,N ≡ µ (cid:88) b Γ jb g (cid:107)(cid:107) M,N , (B37) ε ( µ ) (2) M,N ≡ µ (cid:88) b Γ jb h (cid:107) M,N (∆ v )2 , (B38) (cid:37) ( µ ) (1) M,N ≡ µ (cid:88) b Γ jb g (cid:107)⊥ M,N , (B39) (cid:37) ( µ ) (2) M,N ≡ µ (cid:88) b Γ jb h ⊥ M,N (∆ v )2 , (B40)and τ ( µ ) M,N ≡ µ (cid:88) b Γ jb g ⊥⊥ M,N . (B41)Like Eq. (A28), Eq. (B31) provides the time evolution of f TM,N in one time step simultaneously in the v ⊥ - and v (cid:107) -spaces.Therefore, it represents the numerical solution of Eq. (40), which describes the action of Coulomb collisions of particlesin f j with particles in f b .8 REFERENCES Albert, J. M. 2004, Space Weather, 2,doi: 10.1029/2004SW000069Brügmann, B., Tichy, W., & Jansen, N. 2004, PhysicalReview Letters, 92, doi: 10.1103/physrevlett.92.211101Feldman, W. C., Anderson, R. C., Bame, S. J., et al. 1983,Journal of Geophysical Research: Space Physics, 88, 9949,doi: 10.1029/JA088iA12p09949Fitzenreiter, R. J., Ogilvie, K. W., Bale, S. D., & Viñas,A. F. 2003, Journal of Geophysical Research: SpacePhysics, 108, doi: 10.1029/2003JA009865Gary, S. P., Saito, S., & Li, H. 2008, Geophysical ResearchLetters, 35, doi: 10.1029/2007GL032327Gendrin, R. 1968, Journal of Atmospheric and TerrestrialPhysics, 30, 1313, doi: 10.1016/S0021-9169(68)91158-6Gendrin, R. 1981, Reviews of Geophysics, 19, 171,doi: 10.1029/RG019i001p00171Gendrin, R., & Roux, A. 1980, Journal of GeophysicalResearch: Space Physics, 85, 4577,doi: 10.1029/JA085iA09p04577Gilat, A., & Subramaniam, V. V. 2011, Numerical methods :an introduction with applications using MATLAB (WileyPublishing)Glauert, S. A., & Horne, R. B. 2005, Journal of GeophysicalResearch: Space Physics, 110, doi: 10.1029/2004JA010851Graham, G. A., Rae, I. J., Owen, C. J., et al. 2017, Journalof Geophysical Research: Space Physics, 122, 3858,doi: 10.1002/2016JA023656Gurgiolo, C., & Goldstein, M. L. 2016, Annales Geophysicae,34, 1175, doi: 10.5194/angeo-34-1175-2016Gurgiolo, C., Goldstein, M. L., Viñas, A. F., & Fazakerley,A. N. 2012, Annales Geophysicae, 30, 163,doi: 10.5194/angeo-30-163-2012Gurnett, D. A., & Bhattacharjee, A. 2017, Introduction toPlasma Physics: With Space, Laboratory andAstrophysical Applications, Cambridge University Press,doi: 10.1017/9781139226059Heuer, M., & Marsch, E. 2007, Journal of GeophysicalResearch: Space Physics, 112, doi: 10.1029/2006JA011979Horaites, K., Boldyrev, S., & Medvedev, M. V. 2019,MNRAS, 484, 2474, doi: 10.1093/mnras/sty3504Horaites, K., Boldyrev, S., Wilson, Lynn B., I., Viñas, A. F.,& Merka, J. 2018, MNRAS, 474, 115,doi: 10.1093/mnras/stx2555Isenberg, P., & Vasquez, B. 2011, The Astrophysical Journal,731, 88, doi: 10.1088/0004-637X/731/2/88Isenberg, P. A., & Lee, M. A. 1996, Journal of GeophysicalResearch: Space Physics, 101, 11055,doi: 10.1029/96JA00293 Iserles, A. 2008, A First Course in the Numerical Analysis ofDifferential Equations, Cambridge University Press,doi: 10.1017/CBO9780511995569Kennel, C. F., & Engelmann, F. 1966, Physics of Fluids, 9,2377, doi: 10.1063/1.1761629Khazanov, G. V., Gamayunov, K. V., Jordanova, V. K., &Krivorutsky, E. N. 2002, Journal of Geophysical Research:Space Physics, 107, SMP 14, doi: 10.1029/2001JA000180Klein, K. G., & Chandran, B. D. G. 2016, ApJ, 820, 47,doi: 10.3847/0004-637X/820/1/47Lacombe, C., Alexandrova, O., Matteini, L., et al. 2014,ApJ, 796, 5, doi: 10.1088/0004-637X/796/1/5Lemons, D. S., & Feldman, W. C. 1983, J. Geophys. Res.,88, 6881, doi: 10.1029/JA088iA09p06881Ljepojevic, N. N., Burgess, A., & Moffatt, H. K. 1990,Proceedings of the Royal Society of London. A.Mathematical and Physical Sciences, 428, 71,doi: 10.1098/rspa.1990.0026Lyons, L. R. 1974, Journal of Plasma Physics, 12, 45–49,doi: 10.1017/S0022377800024910Lyons, L. R., Thorne, R. M., & Kennel, C. F. 1971, Journalof Plasma Physics, 6, 589–606,doi: 10.1017/S0022377800006310Lyutikov, M., & Gavriil, F. P. 2006, Monthly Notices of theRoyal Astronomical Society, 368, 690,doi: 10.1111/j.1365-2966.2006.10140.xMarsch, E. 2006, Living Reviews in Solar Physics, 3, 1,doi: 10.12942/lrsp-2006-1Marsch, E., & Bourouaine, S. 2011, Annales Geophysicae, 29,2089, doi: 10.5194/angeo-29-2089-2011Marsch, E., & Tu, C.-Y. 2001, Journal of GeophysicalResearch: Space Physics, 106, 8357,doi: 10.1029/2000JA000414Mourenas, D., Artemyev, A. V., Agapitov, O. V.,Krasnoselskikh, V., & Mozer, F. S. 2015, Journal ofGeophysical Research: Space Physics, 120, 3665,doi: 10.1002/2015JA021135Pagel, C., Gary, S. P., de Koning, C. A., Skoug, R. M., &Steinberg, J. T. 2007, Journal of Geophysical Research:Space Physics, 112, doi: 10.1029/2006JA011967Pilipp, W. G., Miggenrieder, H., Montgomery, M. D., et al.1987, Journal of Geophysical Research: Space Physics, 92,1075, doi: 10.1029/JA092iA02p01075Roberg-Clark, G. T., Drake, J. F., Reynolds, C. S., &Swisdak, M. 2016, The Astrophysical Journal, 830, L9,doi: 10.3847/2041-8205/830/1/l9Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957,Phys. Rev., 107, 1, doi: 10.1103/PhysRev.107.19