- Email: [email protected]
The coupling of quasi-linear pitch angle and energy diffusion
ARTICLE IN PRESS Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
Contents lists available at ScienceDirect
Journal of Atmospheric and Solar-Terrestrial Physics journal homepage: www.elsevier.com/locate/jastp
The coupling of quasi-linear pitch angle and energy diffusion Jay M. Albert Air Force Research Laboratory, Space Vehicles Directorate, 29 Randolph Road, Hanscom AFB, MA 01731-3010, USA
a r t i c l e in f o
a b s t r a c t
Article history: Accepted 28 November 2008 Available online 24 December 2008
Wave–particle interactions are described in the quasi-linear formulation by bounce-averaged diffusion coefficients for equatorial pitch angle a0 and for momentum p, along with a mixed or cross diffusion coefficient Da0 p . Because the cross terms complicate the associated time-dependent diffusion equation, and because Da0 p has a somewhat unfamiliar character (e.g., it may be negative), it has frequently been omitted in numerical simulations of multidimensional diffusion. Generally, Da0 p becomes increasingly important for small values of f pe =f ce , as does Dpp , and for small widths of the wave frequency and wavenormal angle distributions. Here we use very simple models of diffusion to investigate the effect of the cross terms, and to demonstrate numerical problems associated with them. Published by Elsevier Ltd.
Keywords: Cross diffusion Diffusion coefficients Wave–particle interactions
1. Introduction Cyclotron-resonant wave–particle interactions are fundamental to the behavior of radiation belt electrons, having long been identified as the dominant loss term for inner zone electrons through pitch angle diffusion into the loss cone (Lyons et al., 1972; Abel and Thorne, 1998). More recently, they have been implicated in the dynamics of outer zone electrons both as a loss mechanism (Lorentzen et al., 2001; Thorne et al., 2005) and as a source via energy diffusion from a less energetic seed population (Summers et al., 1998), particularly by whistler mode chorus waves (Horne et al., 2003, 2005). In a strong background magnetic field, such as that of the Earth, the conventional formulation of quasi-linear theory uses spherical coordinates for the particle momentum, which yields a diffusion equation which separates into terms for pitch angle a, momentum p, and ‘‘cross’’ or ‘‘mixed’’ diffusion. All three diffusion coefficients arise from the same cyclotron-resonant process, and have definite relationships which, after bounce averaging, typically lead to the ordering Da0 a0 4jDa0 p j=p4Dpp =p2 . Although Da0 p is no more difficult to calculate than the other two, cross diffusion greatly complicates numerical solution of the diffusion equation (which would otherwise be straightforward). Because of the ordering, it seems questionable that cross diffusion may be neglected if one is interested in the effects of energy (or equivalently, p) diffusion. Recently, Albert and Young (2005) presented a technique for constructing new variables in which the cross diffusion term rigorously vanishes; the physical effects are retained by the coordinate transformations. Diffusion coefficients for electrons at L ¼ 4:5 interacting with a realistic model of storm-time chorus waves were used to simulate the production of MeV electrons on a
one day time scale, confirming earlier estimates. The associated, time-developing pitch angle distribution was also naturally provided by the numerical solution for f ða0 ; p; tÞ. Fully three-dimensional simulations, including the effects of radial diffusion, are still under development, although useful estimates without accounting for cross diffusion have been made (Varotsou et al., 2005). This paper discusses the relationship between cross diffusion and ‘‘pure’’ diffusion, and the dependence on corresponding physical parameters. Connection is made to resonant ellipses and to single-wave characteristics (which are also ellipses) and their generalization, resonant diffusion curves. Then, a highly simplified two-dimensional model problem is introduced. The well-known Green function solutions for an infinite domain (whose contours are also ellipses) are examined as the relative strength of cross diffusion is varied. Simple numerical methods are used on a finite-domain version of the model problem, demonstrating both the effects of cross diffusion and the numerical problems that arise. The origin of these problems for the simple finite difference algorithm is isolated in the cross terms and expected gradients of f and shown to be inevitable, motivating the use of more sophisticated approaches.
2. Formulation The two-dimensional bounce-averaged diffusion equation may be written as either #" # " qf =qJ1 D11 D12 qf q q (1) ¼ qt qJ 1 qJ2 D12 D22 qf =qJ2 or
E-mail address: [email protected] 1364-6826/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.jastp.2008.11.014
"
Da0 a0 qf 1 q q ¼ G Da 0 p qt G qa0 qp
Da0 p Dpp
#"
#
qf =qa0 . qf =qp
(2)
ARTICLE IN PRESS J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
Here G ¼ p2 Tða0 Þ sin a0 cos a0 , and Tða0 Þ is the normalized bounce period ðc=4LRe Þtb (Schulz, 1991). It is important to note that in this notation Da0 a0 has dimensions of a20 =t, while Dpp has dimensions of p2 =t. We also consider the simple model equation # #" " qf q q Dxx Dxy qf =qx , (3) ¼ qt qx qy Dxy Dyy qf =qy where the matrix elements are constant and satisfy Dxx Dyy D2xy 40. 2.1. Normalizations In a dipole magnetic field, the three adiabatic invariants may be written as (Schulz and Lanzerotti, 1974, chapter 1) ^J ¼ L3 p^ 2 sin2 a0 ; 1
^J ¼ LpYðsin ^ a0 Þ; 2
^J ¼ 1 , 3 L
(4)
where p^ ¼ p=mc and ^J ¼ J1 ; 1 J 10 J 10 ¼ p
^J ¼ J 2 ; 2 J20
mc2
O0
diffusion, cross diffusion in ð^J 1 ; ^J 2 Þ is essentially as large as the ‘‘diagonal’’ terms in ^J 1 and ^J2 and cannot be neglected.
3. Relationships among diffusion coefficients We investigate conditions leading to large values of D2a0 p =ðDa0 a0 Dpp Þ. The standard quasi-linear expressions for the diffusion coefficients are bounce averages of local diffusion coefficients, which are sums (over resonant harmonic number n) of integrals over the resonant waves, usually parameterized by wavenormal angle y. At each local, single, resonance, specified by o kk vk ¼ nOe =g, the integrands ðDnaay ; Dnayp ; Dnppy Þ obey ! sin a cos a ny Dnaay , Dap =p ¼ sin2 a þ nOe =og !2 sin a cos a ny 2 Dnaay , (10) Dpp =p ¼ sin2 a þ nOe =og where Dnaay is a complicated but positive quantity given by Lyons (1974). The same expressions hold for the normalized quantities ^ ny =p^ 2 . ^ ny =p^ and D D
^J ¼ J 3 , 3 J 30
ap
;
J20 ¼ 2mcRE ;
1665
J 30 ¼ 2pmO0 R2E .
(5)
Here O0 ¼ eB0 =mc, B0 0:31 G is the Earth’s surface magnetic field strength at the equator, and Yðsin a0 Þ is the usual dipole function, which ranges from about 1:74 to 0 as discussed by Schulz and Lanzerotti (1974, p. 20). Thus J10 ; J 20 ; J 30 stand approximately in the ratios 1 : c=O0 RE : ðc=O0 RE Þ2 , with c=O0 RE 105 , but all the normalized variables are considered to be 1. Normalization by constants does not change the form of the diffusion equation. For example,
q qf q ^ qf D ¼ , D qJ 1 12 qJ 2 q^J 1 12 q^J 2
(6)
pp
For whistler mode waves, ooOe (and possibly 5Oe ). Thus negative values of n, which are often dominant, typically give positive values of Dnayp while nX0 gives negative values. Large values of Dnayp =Dnaay are favored by small values (1) of ope =Oe , as pointed out by Horne et al. (2003). Albert (2005) explained this in terms of the o dependence of ðkc=oÞ2 , which typically leads to larger resonant values of o=Oe for smaller ope =Oe values. Note that ðDnayp Þ2 =ðDnaay Dnppy Þ is exactly 1, but D2a0 p =ðDa0 a0 Dpp Þ is typically smaller; reduction occurs through integration and summation over latitude, n, and wavenormal angle, in accordance with the Cauchy–Schwarz inequality (Albert, 2004). The reduction will thus be minimal if resonant waves only occur in narrow ranges of y, latitude, and n. In such cases, D2a0 p will be close to its largest possible value, Da0 a0 Dpp .
^ 12 ¼ D12 =ðJ J Þ, and where D 10 20
q qf q ^ qf ¼ D D , qa0 a0 p qp qa0 a0 p qp^
(7)
^ a p ¼ Da p =mc. Although radial diffusion is being ignored where D 0 0 here, we also note that
q qf q DLL qf . D ¼ L2 qJ 3 33 qJ 3 qL L2 qL
(8)
These simple relations allow the diffusion equation to be written in the form of Eq. (1) or (2) with all the matrix elements nondimensional and comparable in numerical value. 2.2. Diffusion in ðJ1 ; J2 Þ Taking into account the normalization discussion above, we ^ relate the diffusion coefficients in ð^J1 ; ^J 2 Þ to those in ða0 ; pÞ: ^ 11 ¼ D
q^J1 qa0
!2 Da0 a0 þ 2
^ ^ ^ 12 ¼ qJ 1 qJ 2 Da a þ D 0 0
qa0 qa0
^ 22 ¼ D
q^J2 qa0
!2 Da0 a0 þ 2
q^J 1 q^J 1 ^ q^J1 Da 0 p þ ^ qa0 qp qp^
!2 ^ pp , D
!
q^J 1 q^J 2 q^J 1 q^J 2 ^ q^J q^J ^ D þ 1 2D , þ qa0 qp^ qp^ qa0 a0 p qp^ qp^ pp q^J 2 q^J 2 ^ q^J2 Da 0 p þ ^ qa0 qp qp^
!2 ^ pp , D
(9)
From Eq. (4), the partial derivatives are all 1 for radiation belt particles. Thus even for the limiting case of pure pitch angle
3.1. Resonant diffusion curves Eqs. (10) derive from the more fundamental single wave characteristics, which express the constancy of particle energy in the frame of a single resonant wave, as reviewed by Walker (1993). They are ellipses in ðv? ; vk Þ (not to be confused with the resonant ellipses, which simply express the resonance condition). In the presence of a continuous band of waves, they are still valid in differential form, for the instantaneously resonant individual wave. Summers et al. (1998) gave this differential equation in terms of ðv? ; vk Þ; transforming to ðp; aÞ gives a relation between Da and Dp for short times Dt. It thus gives relations between Dnaay ¼ hðDaÞ2 =2Dti, Dnayp ¼ hDaDp=2Dti, and Dnppy ¼ hðDpÞ2 =2Dti, which are exactly Eqs. (10). Integrating the expression for dp=da over the frequency distribution (Summers et al., 1998) gives resonant diffusion curves, which Horne et al. (2003) used to qualitatively analyze electron energization and loss by a band of chorus waves, and the corresponding effect on the waves. The effects associated with cross diffusion are therefore deeply connected to self-consistent wave-particle interactions. However, since these resonant diffusion curves are constructed for single values ofn, y, and l (parallelpropagating waves, with n ¼ 1, at fixed latitude), the diffusion obeys Daa Dpp D2ap ¼ 0 and is one dimensional rather than two dimensional (Albert, 2004). Thus, in the analysis of Horne et al. (2003), phase space density spreads out along one-dimensional paths rather than in two dimensions as described by Eq. (1) or (2).
ARTICLE IN PRESS 1666
J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
1.0
4. Green functions on an infinite domain On an infinite domain, the Green function for Eq. (3), for a distribution initially concentrated at ðx0 ; y0 Þ, is Gðx; y; tÞ ¼ exp½ðDyy ðdxÞ2 2Dxy dxdy pffiffiffiffi þ Dxx ðdyÞ2 Þ=4t L=4pt L,
(11) D2xy .
Contours of where dx ¼ x x0 , dy ¼ y y0 , and L ¼ Dxx Dyy constant G are given by ellipses; it is intuitively clear that the axes give the directions of fastest and slowest diffusion. The major and minor radii of the ellipse are proportional to the eigenvalues of the diffusion matrix, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dxx þ Dyy Dxx Dyy 2 (12) l ¼ þ D2xy , 2 2 and the ellipse is rotated relative to the ðx; yÞ axes by angle tan1 ½Dxy =ðlþ Dxx Þ, to be aligned with the corresponding eigenvectors. Thus Dxy affects the size of the ellipse, its shape (eccentricity), and its orientation. Fig. 1 shows the Green function Gðx; yÞ centered at ð0:5; 0:5Þ, with Dxx t ¼ 2 103 , Dyy t ¼ 1 103 , and several values of Dxy t parameterized by ¼ Dxy =ðDxx Dyy Þ1=2, with 1oo1. The plots, normalized here to have unit maxima, show the increasing eccentricity and rotation of the G contours with increasing magnitude of Dxy . Fig. 2 shows contours of the Green function both with and without cross diffusion. First, G was evaluated with Dxx t ¼ 102 , Dyy t ¼ 0:2 102 , and Dxy t ¼ 0, and its maximum value G0 was noted. The dashed curves show contours with levels ðG=G0 Þ2 ¼ ð0:1; 0:01; 0:001Þ (red, green, blue, respectively). The solid curves differ only in that G was evaluated with Dxy t corresponding to ¼ 0:8. At points of red–green or green–blue intersection the values p ffiffiffiffiffiffi of G with and without cross diffusion differ by a factor of 10, and at red–blue intersections they differ by a factor of 10.
0.5
0.0 0.0
0.5
1.0
Fig. 2. Contours of the diffusion Green function Gðx; yÞ, both with and without cross diffusion. The dashed curves are for ðDxx t; Dyy t; Dxy tÞ ¼ ð1; 0:5; 0Þ 102, showing levels G=G0 ¼ ð0:316; 0:10; 0:0316Þ (red, green, blue, respectively), where G0 is the maximum value of G. The solid curves show the effect of changing Dxy to ðDxx Dyy Þ1=2 with ¼ 0:8.
5. Finite boundaries: numerical solutions The diffusion Green functions considered above are strictly valid only for infinite domains, so we now consider numerical methods on a finite domain, 0px; yp1. As a first experiment, we take as initial conditions a centered Gaussian, namely f ðx; yÞ ¼ exp½ðx 0:5Þ2 =0:01 ðy 0:5Þ2 =0:01. Boundary conditions are
Fig. 3. Numerical solution of a simple diffusion equation on 0px; yp1, with constant diffusion coefficients given in the text. The boundary conditions are f ðx; 0Þ ¼ 0, f ðx; 1Þ ¼ 0, f ð0; yÞ ¼ 0, qf =qxð1; yÞ ¼ 0, and the initial conditions given by a Gaussian peaked at the center. Numerically negative values are indicated by black plus symbols.
Fig. 1. The Green function Gðx; yÞ with Dxx t ¼ 2 103 , Dyy t ¼ 1 103 , and several values of Dxy t given by Dxy ¼ ðDxx Dyy Þ1=2.
chosen to be f ðx; 0Þ ¼ 0, f ðx; 1Þ ¼ 0, f ð0; yÞ ¼ 0, and qf =qxð1; yÞ ¼ 0. The condition at x ¼ 1 is modeled after the realistic boundary condition qf =qa0 ¼ 0 at a0 ¼ 90 . The value of f ð0:5; 0:5Þ is not held fixed but is allowed to decay with time. We use the straightforward spatial discretization and explicit time differencing discussed in Section 6, with 51 grid points, and with the timestep set to 0.01 times the maximum value allowed by linear stability; that is, Dt ¼ 106 . Fig. 3 shows numerical results at a few different times, for cases with Dxx ¼ Dyy ¼ 1 and three different values of Dxy , with ¼ Dxy =ðDxx Dyy Þ1=2 as above. At t ¼ 0:01, f still resembles the infinite-domain results, but by t ¼ 0:05 the effects of the boundaries are apparent. It is also clear that the effect of cross diffusion is significant, particularly at the intermediate times. Note that even
ARTICLE IN PRESS J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
1667
here. The time levels at which f is evaluated depend on the n specific numerical algorithm. Using f everywhere gives a fully explicit scheme, which for constant diffusion coefficients takes the form nþ1
f ij
n
n
n
n
n
¼ ½1 ð2A þ 2BÞf ij þ Af iþ1j þ Af i1j þ Bf ijþ1 þ Bf ij1 þ
n 2Cf iþ1jþ1
n 2Cf i1jþ1
2
n 2Cf iþ1j1
þ
n 2Cf i1j1 .
(15)
2
Here A ¼ Dxx Dt=ðDxÞ , B ¼ Dyy Dt=ðDyÞ , and C ¼ Dxy Dt=ð4DxDyÞ, evaluated at the appropriate grid points. Although this is probably the simplest reasonable finite difference algorithm imaginable, the spatial discretizations are each locally (hence globally) conservative. Furthermore, the usual local, linear analysis indicates that the explicit scheme is stable for sufficiently small timestep, given by ½Dxx =ðDxÞ2 þ Dyy =ðDyÞ2 Dto1=2
Fig. 4. Same as Fig. 3 except that the initial conditions are peaked for large x and small y, as given in the text and shown in the top row, and the boundary values f ðx; 0Þ are held fixed.
this simple simulation, with constant diffusion coefficients and a quite small timestep, leads to negative values of f for ¼ 0:8, as indicated by black plus symbols in Fig. 3. To make the run more closely resemble a realistic simulation in ða0 ; pÞ, the initial conditions were changed to f ðx; yÞ ¼ exp ðy=0:02Þ sinðpx=2Þ, modeling a constant source at low p or energy, and the values at y ¼ 0 were held fixed in time. Fig. 4 shows the results, in the same format as Fig. 3. The f ðx; yÞ results are necessarily identical on the y ¼ 0 boundary, but evolve quite differently at larger y values, especially at low x values, where the large- values of f are seen to be lower by more than a factor of 10. This is consistent with the detailed, more realistic simulations in ða0 ; pÞ performed by Albert and Young (2005), who also found that cross diffusion acts to reduce the rate of diffusion to high energy. As in Fig. 3, negative values of f occur and are indicated by the black plus symbols. Simulations using the same numerical algorithm and physically realistic, widely varying diffusion coefficients show much more drastic development of negative values, and not just near the boundaries.
6. Numerical problems To see how negative values can arise numerically, consider Eq. (3) with constant diffusion coefficients. The left-hand side is nþ1 n most simply discretized in time by qf =qt ! ðf ij f ij Þ=Dt, where, as usual, the superscript refers to the timestep and the subscripts index the ðx; yÞ grid. On the right-hand side of Eq. (3), the term involving Dxx is discretized according to
q=qxðDqf =qxÞ ! ½Diþ1j ðf iþ1j f ij Þ=Dx 2
D
1 ðf i2j ij
f i1j Þ=Dx=Dx,
(13)
and similarly for the term with Dyy , while the terms involving Dxy are discretized by
q=qyðDqf =qxÞ ! ½Dijþ1 ðf iþ1jþ1 f i1jþ1 Þ=2Dx Dij1 ðf iþ1j1 f i1j1 Þ=2Dx=2Dy,
(14)
and similarly for q=qxðDqf =qyÞ. This is preferable to separating Eq. (3) into diffusive and advective pieces before discretizing, since the treatment of advection introduces a whole new class of numerical issues while failing to solve the problems discussed
(16)
(Richtmyer and Morton, 1967). In the words of Press et al. (1992) (in a different context), this is ‘‘a fine example of an algorithm that is easy to derive, takes little storage, and executes quickly. Too bad it doesn’t work!’’ In particular, this algorithm does not guarantee n positivity. For small enough timestep the factors multiplying f ij , n n f i1j , and f ij1 are necessarily positive, but some of the factors n multiplying f i1j1 are not. Thus in the presence of cross diffusion, n positive but unfavorable values of f can yield negative values of nþ1 even with an arbitrarily small timestep. f In the radiation belts, f is usually largest for large a0 and small energy, and the cross diffusion coefficient is typically positive, ^ or ða0 ; pÞ ^ grid, 2Cf iþ1j1 is likely to be negative so on an ða0 ; EÞ and can be dominant. Realistic diffusion coefficients range over several orders of magnitude and Da0 p can abruptly change sign (Albert, 2005), which seems unlikely to help the situation. 6.1. Other algorithms nþ1
n
Using f instead of f in the discretized diffusion operators gives a fully implicit scheme, which (with constant D’s) is linearly stable for any size timestep. Though this is very inefficient to implement except in one dimension, so-called alternating direction implicit (ADI) schemes are much more practical and retain much of the theoretical stability. Variants of ADI which treat cross diffusion have been formulated by Beam and Warming (1980) and Gourlay and McKee (1977), but they have also been found vulnerable to unphysical, negative values. The explicit scheme could be ‘‘fixed’’ with a Lax approach, n nþ1 n which replaces f ij in the time difference ðf ij f ij Þ=Dt with its average over neighboring grid points. This has the effect of deliberately adding grid diffusion to the physical diffusion (Press et al., 1992). This may be reasonable in a problem dominated by other processes, but it is not appropriate when the physical diffusion is the focus of interest. (Since the grid diffusion would have to be comparable to the maximum value of the cross diffusion, smaller values would be completely masked.) Other promising numerical approaches include the supportoperators method of Shashkov and Steinberg (1996), which is optimized for rough diffusion coefficients; the method of lines (Sommeijer et al., 1997) which discretizes the diffusion equation in space but not in time; and a wide variety of finite volume and finite element methods. Tao et al. (2008) reformulated Eq. (2) as a Fokker–Planck equation, which they treated with random walk techniques; this is straightforward, but relatively inefficient. Albert and Young (2005) attacked the problematic cross terms directly, by constructing new variables such that cross terms are absent from the transformed equation. The resulting equation is essentially of the same form as Eq. (1) or (2) with the cross terms simply omitted, but the new variables retain the full physical
ARTICLE IN PRESS 1668
J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
effects. This equation can be treated with simple, well-known finite difference algorithms that are robust and reliable.
7. Summary The significance of pitch angle-energy cross diffusion has been discussed from several points of view. Normalizing the diffusion equation, written in either ða0 ; pÞ or ðJ1 ; J 2 Þ, shows that the diffusion operators for all the variables are present on an equal footing. For single resonant waves, the diffusion coefficients themselves have well-defined ratios Dap =Daa and D2ap =ðDaa Dpp Þ, which are inseparable from the fundamental theoretical concepts of single wave characteristics and resonant diffusion curves. Integrating over the wave distribution, summing over resonances, and bounce averaging modify these values, ratios, and the concept itself of a diffusion path in momentum space (it becomes two-dimensional rather than one-dimensional). A simple two-dimensional model diffusion equation with constant coefficients was written, and its Green function given and illustrated. Its contours are ellipses, whose size, eccentricty, and orientation depend on all three diffusion coefficients. It was demonstrated that the cross term can change the value at a given location by a factor of 10. This equation was then combined with initial and boundary conditions modeled on actual pitch angle and energy diffusion, and solved numerically with a simple, standard algorithm. The solutions again show the effect of the cross terms, but also yield negative values of f , hinting at the bad behavior found when the same algorithm is applied to realistic diffusion coefficients. The problem was traced to the finite differencing of the cross terms; the textbook approach (Richtmyer and Morton, 1967), while numerically stable (for constant coefficients), is inherently unsatisfactory. However, as argued above, these terms are physically meaningful and should be included somehow. Fortunately, there are now several alternative means available to do so.
Acknowledgments This work was supported by the Space Vehicles Directorate of the Air Force Research Laboratory.
References Abel, B., Thorne, R.M., 1998. Electron scattering loss in Earth’s inner magnetosphere, 1, dominant physical processes. Journal of Geophysical Research 103, 2385. Albert, J.M., 2004. Using quasi-linear diffusion to model acceleration and loss from wave–particle interactions. Space Weather 2, S09S03. Albert, J.M., 2005. Evaluation of quasi-linear diffusion coefficients for whistler mode waves in a low density plasma. Journal of Geophysical Research 110, A03218. Albert, J.M., Young, S.L., 2005. Multidimensional quasi-linear diffusion of radiation belt electrons. Geophysical Research Letters 32, L14110. Beam, R.M., Warming, R.F., 1980. Alternating direct implicit methods for parabolic equations with a mixed derivative. SIAM Journal on Scientific and Statistical Computing 1, 131. Gourlay, A.R., McKee, S., 1977. The construction of hopscotch methods for parabolic and elliptic equations in two space dimensions with a mixed derivative. Journal of Computational and Applied Mathematics 3, 201. Horne, R.B., Glauert, S.A., Thorne, R.M., 2003. Resonant diffusion of radiation belt electrons by whistler-mode chorus. Geophysical Research Letters 30 (9), 1493. Horne, R.B., Thorne, R.M., Glauert, S.A., Albert, J.M., Meredith, N.P., Anderson, R.R., 2005. Timescale for radiation belt electron acceleration by whistler mode chorus waves. Journal of Geophysical Research 110, A03225. Lorentzen, K.R., Blake, J.B., Inan, U.S., Bortnik, J., 2001. Observations of relativistic electron microbursts in association with VLF chorus. Journal of Geophysical Research 106, 6017. Lyons, L.R., 1974. Pitch angle and energy diffusion coefficients from resonant interactions with ion-cyclotron and whistler waves. Journal of Plasma Physics 12, 417. Lyons, L.R., Thorne, R.M., Kennel, C.F., 1972. Pitch-angle diffusion of radiation belt electrons within the plasmasphere. Journal of Geophysical Research 77, 3455. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes. Cambridge University Press, New York. Richtmyer, R.D., Morton, K.W., 1967. Difference Methods for Initial-Value Problems. Wiley, New York. Schulz, M., 1991. The magnetosphere. In: Jacobs, J.A. (Ed.), Geomagnetism, vol. 4. Academic, San Diego, CA, pp. 87–293. Schulz, M., Lanzerotti, L.J., 1974. Particle Diffusion in the Radiation Belts. Springer, New York. Shashkov, M., Steinberg, S., 1996. Solving diffusion equations with rough coefficients in rough grids. Journal of Computational Physics 129, 383. Sommeijer, B.P., Shampine, L.F., Verwer, J.G., 1997. RKC: an explicit solver for parabolic PDEs. Journal of Computational and Applied Mathematics 88, 315. Summers, D., Thorne, R.M., Xiao, F., 1998. Relativistic theory of wave-particle resonant diffusion with application to electron acceleration in the magnetosphere. Journal of Geophysical Research 103(A9) (20), 487. Tao, X., Chan, A.A., Albert, J.M., Miller, J.A., 2008 Stochastic modeling of multidimensional diffusion in the radiation belts. Journal of Geophysical Research 113, A07212. Thorne, R.M., O’Brien, T.P., Shprits, Y.Y., Summers, D., Horne, R.B., 2005. Timescale for MeV electron microburst loss during geomagnetic storms. Journal of Geophysical Research 110, A09202. Varotsou, A., Boscher, D., Bourdarie, S., Horne, R.B., Glauert, S.A., Meredith, N.P., 2005. Simulation of the outer radiation belt electrons near geosynchronous orbit including both radial diffusion and resonant interaction with whistlermode chorus waves. Geophysical Research Letters 32, L19106. Walker, A.D.M., 1993. Plasma Waves in the Magnetosphere. Springer, New York.
Contents lists available at ScienceDirect
Journal of Atmospheric and Solar-Terrestrial Physics journal homepage: www.elsevier.com/locate/jastp
The coupling of quasi-linear pitch angle and energy diffusion Jay M. Albert Air Force Research Laboratory, Space Vehicles Directorate, 29 Randolph Road, Hanscom AFB, MA 01731-3010, USA
a r t i c l e in f o
a b s t r a c t
Article history: Accepted 28 November 2008 Available online 24 December 2008
Wave–particle interactions are described in the quasi-linear formulation by bounce-averaged diffusion coefficients for equatorial pitch angle a0 and for momentum p, along with a mixed or cross diffusion coefficient Da0 p . Because the cross terms complicate the associated time-dependent diffusion equation, and because Da0 p has a somewhat unfamiliar character (e.g., it may be negative), it has frequently been omitted in numerical simulations of multidimensional diffusion. Generally, Da0 p becomes increasingly important for small values of f pe =f ce , as does Dpp , and for small widths of the wave frequency and wavenormal angle distributions. Here we use very simple models of diffusion to investigate the effect of the cross terms, and to demonstrate numerical problems associated with them. Published by Elsevier Ltd.
Keywords: Cross diffusion Diffusion coefficients Wave–particle interactions
1. Introduction Cyclotron-resonant wave–particle interactions are fundamental to the behavior of radiation belt electrons, having long been identified as the dominant loss term for inner zone electrons through pitch angle diffusion into the loss cone (Lyons et al., 1972; Abel and Thorne, 1998). More recently, they have been implicated in the dynamics of outer zone electrons both as a loss mechanism (Lorentzen et al., 2001; Thorne et al., 2005) and as a source via energy diffusion from a less energetic seed population (Summers et al., 1998), particularly by whistler mode chorus waves (Horne et al., 2003, 2005). In a strong background magnetic field, such as that of the Earth, the conventional formulation of quasi-linear theory uses spherical coordinates for the particle momentum, which yields a diffusion equation which separates into terms for pitch angle a, momentum p, and ‘‘cross’’ or ‘‘mixed’’ diffusion. All three diffusion coefficients arise from the same cyclotron-resonant process, and have definite relationships which, after bounce averaging, typically lead to the ordering Da0 a0 4jDa0 p j=p4Dpp =p2 . Although Da0 p is no more difficult to calculate than the other two, cross diffusion greatly complicates numerical solution of the diffusion equation (which would otherwise be straightforward). Because of the ordering, it seems questionable that cross diffusion may be neglected if one is interested in the effects of energy (or equivalently, p) diffusion. Recently, Albert and Young (2005) presented a technique for constructing new variables in which the cross diffusion term rigorously vanishes; the physical effects are retained by the coordinate transformations. Diffusion coefficients for electrons at L ¼ 4:5 interacting with a realistic model of storm-time chorus waves were used to simulate the production of MeV electrons on a
one day time scale, confirming earlier estimates. The associated, time-developing pitch angle distribution was also naturally provided by the numerical solution for f ða0 ; p; tÞ. Fully three-dimensional simulations, including the effects of radial diffusion, are still under development, although useful estimates without accounting for cross diffusion have been made (Varotsou et al., 2005). This paper discusses the relationship between cross diffusion and ‘‘pure’’ diffusion, and the dependence on corresponding physical parameters. Connection is made to resonant ellipses and to single-wave characteristics (which are also ellipses) and their generalization, resonant diffusion curves. Then, a highly simplified two-dimensional model problem is introduced. The well-known Green function solutions for an infinite domain (whose contours are also ellipses) are examined as the relative strength of cross diffusion is varied. Simple numerical methods are used on a finite-domain version of the model problem, demonstrating both the effects of cross diffusion and the numerical problems that arise. The origin of these problems for the simple finite difference algorithm is isolated in the cross terms and expected gradients of f and shown to be inevitable, motivating the use of more sophisticated approaches.
2. Formulation The two-dimensional bounce-averaged diffusion equation may be written as either #" # " qf =qJ1 D11 D12 qf q q (1) ¼ qt qJ 1 qJ2 D12 D22 qf =qJ2 or
E-mail address: [email protected] 1364-6826/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.jastp.2008.11.014
"
Da0 a0 qf 1 q q ¼ G Da 0 p qt G qa0 qp
Da0 p Dpp
#"
#
qf =qa0 . qf =qp
(2)
ARTICLE IN PRESS J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
Here G ¼ p2 Tða0 Þ sin a0 cos a0 , and Tða0 Þ is the normalized bounce period ðc=4LRe Þtb (Schulz, 1991). It is important to note that in this notation Da0 a0 has dimensions of a20 =t, while Dpp has dimensions of p2 =t. We also consider the simple model equation # #" " qf q q Dxx Dxy qf =qx , (3) ¼ qt qx qy Dxy Dyy qf =qy where the matrix elements are constant and satisfy Dxx Dyy D2xy 40. 2.1. Normalizations In a dipole magnetic field, the three adiabatic invariants may be written as (Schulz and Lanzerotti, 1974, chapter 1) ^J ¼ L3 p^ 2 sin2 a0 ; 1
^J ¼ LpYðsin ^ a0 Þ; 2
^J ¼ 1 , 3 L
(4)
where p^ ¼ p=mc and ^J ¼ J1 ; 1 J 10 J 10 ¼ p
^J ¼ J 2 ; 2 J20
mc2
O0
diffusion, cross diffusion in ð^J 1 ; ^J 2 Þ is essentially as large as the ‘‘diagonal’’ terms in ^J 1 and ^J2 and cannot be neglected.
3. Relationships among diffusion coefficients We investigate conditions leading to large values of D2a0 p =ðDa0 a0 Dpp Þ. The standard quasi-linear expressions for the diffusion coefficients are bounce averages of local diffusion coefficients, which are sums (over resonant harmonic number n) of integrals over the resonant waves, usually parameterized by wavenormal angle y. At each local, single, resonance, specified by o kk vk ¼ nOe =g, the integrands ðDnaay ; Dnayp ; Dnppy Þ obey ! sin a cos a ny Dnaay , Dap =p ¼ sin2 a þ nOe =og !2 sin a cos a ny 2 Dnaay , (10) Dpp =p ¼ sin2 a þ nOe =og where Dnaay is a complicated but positive quantity given by Lyons (1974). The same expressions hold for the normalized quantities ^ ny =p^ 2 . ^ ny =p^ and D D
^J ¼ J 3 , 3 J 30
ap
;
J20 ¼ 2mcRE ;
1665
J 30 ¼ 2pmO0 R2E .
(5)
Here O0 ¼ eB0 =mc, B0 0:31 G is the Earth’s surface magnetic field strength at the equator, and Yðsin a0 Þ is the usual dipole function, which ranges from about 1:74 to 0 as discussed by Schulz and Lanzerotti (1974, p. 20). Thus J10 ; J 20 ; J 30 stand approximately in the ratios 1 : c=O0 RE : ðc=O0 RE Þ2 , with c=O0 RE 105 , but all the normalized variables are considered to be 1. Normalization by constants does not change the form of the diffusion equation. For example,
q qf q ^ qf D ¼ , D qJ 1 12 qJ 2 q^J 1 12 q^J 2
(6)
pp
For whistler mode waves, ooOe (and possibly 5Oe ). Thus negative values of n, which are often dominant, typically give positive values of Dnayp while nX0 gives negative values. Large values of Dnayp =Dnaay are favored by small values (1) of ope =Oe , as pointed out by Horne et al. (2003). Albert (2005) explained this in terms of the o dependence of ðkc=oÞ2 , which typically leads to larger resonant values of o=Oe for smaller ope =Oe values. Note that ðDnayp Þ2 =ðDnaay Dnppy Þ is exactly 1, but D2a0 p =ðDa0 a0 Dpp Þ is typically smaller; reduction occurs through integration and summation over latitude, n, and wavenormal angle, in accordance with the Cauchy–Schwarz inequality (Albert, 2004). The reduction will thus be minimal if resonant waves only occur in narrow ranges of y, latitude, and n. In such cases, D2a0 p will be close to its largest possible value, Da0 a0 Dpp .
^ 12 ¼ D12 =ðJ J Þ, and where D 10 20
q qf q ^ qf ¼ D D , qa0 a0 p qp qa0 a0 p qp^
(7)
^ a p ¼ Da p =mc. Although radial diffusion is being ignored where D 0 0 here, we also note that
q qf q DLL qf . D ¼ L2 qJ 3 33 qJ 3 qL L2 qL
(8)
These simple relations allow the diffusion equation to be written in the form of Eq. (1) or (2) with all the matrix elements nondimensional and comparable in numerical value. 2.2. Diffusion in ðJ1 ; J2 Þ Taking into account the normalization discussion above, we ^ relate the diffusion coefficients in ð^J1 ; ^J 2 Þ to those in ða0 ; pÞ: ^ 11 ¼ D
q^J1 qa0
!2 Da0 a0 þ 2
^ ^ ^ 12 ¼ qJ 1 qJ 2 Da a þ D 0 0
qa0 qa0
^ 22 ¼ D
q^J2 qa0
!2 Da0 a0 þ 2
q^J 1 q^J 1 ^ q^J1 Da 0 p þ ^ qa0 qp qp^
!2 ^ pp , D
!
q^J 1 q^J 2 q^J 1 q^J 2 ^ q^J q^J ^ D þ 1 2D , þ qa0 qp^ qp^ qa0 a0 p qp^ qp^ pp q^J 2 q^J 2 ^ q^J2 Da 0 p þ ^ qa0 qp qp^
!2 ^ pp , D
(9)
From Eq. (4), the partial derivatives are all 1 for radiation belt particles. Thus even for the limiting case of pure pitch angle
3.1. Resonant diffusion curves Eqs. (10) derive from the more fundamental single wave characteristics, which express the constancy of particle energy in the frame of a single resonant wave, as reviewed by Walker (1993). They are ellipses in ðv? ; vk Þ (not to be confused with the resonant ellipses, which simply express the resonance condition). In the presence of a continuous band of waves, they are still valid in differential form, for the instantaneously resonant individual wave. Summers et al. (1998) gave this differential equation in terms of ðv? ; vk Þ; transforming to ðp; aÞ gives a relation between Da and Dp for short times Dt. It thus gives relations between Dnaay ¼ hðDaÞ2 =2Dti, Dnayp ¼ hDaDp=2Dti, and Dnppy ¼ hðDpÞ2 =2Dti, which are exactly Eqs. (10). Integrating the expression for dp=da over the frequency distribution (Summers et al., 1998) gives resonant diffusion curves, which Horne et al. (2003) used to qualitatively analyze electron energization and loss by a band of chorus waves, and the corresponding effect on the waves. The effects associated with cross diffusion are therefore deeply connected to self-consistent wave-particle interactions. However, since these resonant diffusion curves are constructed for single values ofn, y, and l (parallelpropagating waves, with n ¼ 1, at fixed latitude), the diffusion obeys Daa Dpp D2ap ¼ 0 and is one dimensional rather than two dimensional (Albert, 2004). Thus, in the analysis of Horne et al. (2003), phase space density spreads out along one-dimensional paths rather than in two dimensions as described by Eq. (1) or (2).
ARTICLE IN PRESS 1666
J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
1.0
4. Green functions on an infinite domain On an infinite domain, the Green function for Eq. (3), for a distribution initially concentrated at ðx0 ; y0 Þ, is Gðx; y; tÞ ¼ exp½ðDyy ðdxÞ2 2Dxy dxdy pffiffiffiffi þ Dxx ðdyÞ2 Þ=4t L=4pt L,
(11) D2xy .
Contours of where dx ¼ x x0 , dy ¼ y y0 , and L ¼ Dxx Dyy constant G are given by ellipses; it is intuitively clear that the axes give the directions of fastest and slowest diffusion. The major and minor radii of the ellipse are proportional to the eigenvalues of the diffusion matrix, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dxx þ Dyy Dxx Dyy 2 (12) l ¼ þ D2xy , 2 2 and the ellipse is rotated relative to the ðx; yÞ axes by angle tan1 ½Dxy =ðlþ Dxx Þ, to be aligned with the corresponding eigenvectors. Thus Dxy affects the size of the ellipse, its shape (eccentricity), and its orientation. Fig. 1 shows the Green function Gðx; yÞ centered at ð0:5; 0:5Þ, with Dxx t ¼ 2 103 , Dyy t ¼ 1 103 , and several values of Dxy t parameterized by ¼ Dxy =ðDxx Dyy Þ1=2, with 1oo1. The plots, normalized here to have unit maxima, show the increasing eccentricity and rotation of the G contours with increasing magnitude of Dxy . Fig. 2 shows contours of the Green function both with and without cross diffusion. First, G was evaluated with Dxx t ¼ 102 , Dyy t ¼ 0:2 102 , and Dxy t ¼ 0, and its maximum value G0 was noted. The dashed curves show contours with levels ðG=G0 Þ2 ¼ ð0:1; 0:01; 0:001Þ (red, green, blue, respectively). The solid curves differ only in that G was evaluated with Dxy t corresponding to ¼ 0:8. At points of red–green or green–blue intersection the values p ffiffiffiffiffiffi of G with and without cross diffusion differ by a factor of 10, and at red–blue intersections they differ by a factor of 10.
0.5
0.0 0.0
0.5
1.0
Fig. 2. Contours of the diffusion Green function Gðx; yÞ, both with and without cross diffusion. The dashed curves are for ðDxx t; Dyy t; Dxy tÞ ¼ ð1; 0:5; 0Þ 102, showing levels G=G0 ¼ ð0:316; 0:10; 0:0316Þ (red, green, blue, respectively), where G0 is the maximum value of G. The solid curves show the effect of changing Dxy to ðDxx Dyy Þ1=2 with ¼ 0:8.
5. Finite boundaries: numerical solutions The diffusion Green functions considered above are strictly valid only for infinite domains, so we now consider numerical methods on a finite domain, 0px; yp1. As a first experiment, we take as initial conditions a centered Gaussian, namely f ðx; yÞ ¼ exp½ðx 0:5Þ2 =0:01 ðy 0:5Þ2 =0:01. Boundary conditions are
Fig. 3. Numerical solution of a simple diffusion equation on 0px; yp1, with constant diffusion coefficients given in the text. The boundary conditions are f ðx; 0Þ ¼ 0, f ðx; 1Þ ¼ 0, f ð0; yÞ ¼ 0, qf =qxð1; yÞ ¼ 0, and the initial conditions given by a Gaussian peaked at the center. Numerically negative values are indicated by black plus symbols.
Fig. 1. The Green function Gðx; yÞ with Dxx t ¼ 2 103 , Dyy t ¼ 1 103 , and several values of Dxy t given by Dxy ¼ ðDxx Dyy Þ1=2.
chosen to be f ðx; 0Þ ¼ 0, f ðx; 1Þ ¼ 0, f ð0; yÞ ¼ 0, and qf =qxð1; yÞ ¼ 0. The condition at x ¼ 1 is modeled after the realistic boundary condition qf =qa0 ¼ 0 at a0 ¼ 90 . The value of f ð0:5; 0:5Þ is not held fixed but is allowed to decay with time. We use the straightforward spatial discretization and explicit time differencing discussed in Section 6, with 51 grid points, and with the timestep set to 0.01 times the maximum value allowed by linear stability; that is, Dt ¼ 106 . Fig. 3 shows numerical results at a few different times, for cases with Dxx ¼ Dyy ¼ 1 and three different values of Dxy , with ¼ Dxy =ðDxx Dyy Þ1=2 as above. At t ¼ 0:01, f still resembles the infinite-domain results, but by t ¼ 0:05 the effects of the boundaries are apparent. It is also clear that the effect of cross diffusion is significant, particularly at the intermediate times. Note that even
ARTICLE IN PRESS J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
1667
here. The time levels at which f is evaluated depend on the n specific numerical algorithm. Using f everywhere gives a fully explicit scheme, which for constant diffusion coefficients takes the form nþ1
f ij
n
n
n
n
n
¼ ½1 ð2A þ 2BÞf ij þ Af iþ1j þ Af i1j þ Bf ijþ1 þ Bf ij1 þ
n 2Cf iþ1jþ1
n 2Cf i1jþ1
2
n 2Cf iþ1j1
þ
n 2Cf i1j1 .
(15)
2
Here A ¼ Dxx Dt=ðDxÞ , B ¼ Dyy Dt=ðDyÞ , and C ¼ Dxy Dt=ð4DxDyÞ, evaluated at the appropriate grid points. Although this is probably the simplest reasonable finite difference algorithm imaginable, the spatial discretizations are each locally (hence globally) conservative. Furthermore, the usual local, linear analysis indicates that the explicit scheme is stable for sufficiently small timestep, given by ½Dxx =ðDxÞ2 þ Dyy =ðDyÞ2 Dto1=2
Fig. 4. Same as Fig. 3 except that the initial conditions are peaked for large x and small y, as given in the text and shown in the top row, and the boundary values f ðx; 0Þ are held fixed.
this simple simulation, with constant diffusion coefficients and a quite small timestep, leads to negative values of f for ¼ 0:8, as indicated by black plus symbols in Fig. 3. To make the run more closely resemble a realistic simulation in ða0 ; pÞ, the initial conditions were changed to f ðx; yÞ ¼ exp ðy=0:02Þ sinðpx=2Þ, modeling a constant source at low p or energy, and the values at y ¼ 0 were held fixed in time. Fig. 4 shows the results, in the same format as Fig. 3. The f ðx; yÞ results are necessarily identical on the y ¼ 0 boundary, but evolve quite differently at larger y values, especially at low x values, where the large- values of f are seen to be lower by more than a factor of 10. This is consistent with the detailed, more realistic simulations in ða0 ; pÞ performed by Albert and Young (2005), who also found that cross diffusion acts to reduce the rate of diffusion to high energy. As in Fig. 3, negative values of f occur and are indicated by the black plus symbols. Simulations using the same numerical algorithm and physically realistic, widely varying diffusion coefficients show much more drastic development of negative values, and not just near the boundaries.
6. Numerical problems To see how negative values can arise numerically, consider Eq. (3) with constant diffusion coefficients. The left-hand side is nþ1 n most simply discretized in time by qf =qt ! ðf ij f ij Þ=Dt, where, as usual, the superscript refers to the timestep and the subscripts index the ðx; yÞ grid. On the right-hand side of Eq. (3), the term involving Dxx is discretized according to
q=qxðDqf =qxÞ ! ½Diþ1j ðf iþ1j f ij Þ=Dx 2
D
1 ðf i2j ij
f i1j Þ=Dx=Dx,
(13)
and similarly for the term with Dyy , while the terms involving Dxy are discretized by
q=qyðDqf =qxÞ ! ½Dijþ1 ðf iþ1jþ1 f i1jþ1 Þ=2Dx Dij1 ðf iþ1j1 f i1j1 Þ=2Dx=2Dy,
(14)
and similarly for q=qxðDqf =qyÞ. This is preferable to separating Eq. (3) into diffusive and advective pieces before discretizing, since the treatment of advection introduces a whole new class of numerical issues while failing to solve the problems discussed
(16)
(Richtmyer and Morton, 1967). In the words of Press et al. (1992) (in a different context), this is ‘‘a fine example of an algorithm that is easy to derive, takes little storage, and executes quickly. Too bad it doesn’t work!’’ In particular, this algorithm does not guarantee n positivity. For small enough timestep the factors multiplying f ij , n n f i1j , and f ij1 are necessarily positive, but some of the factors n multiplying f i1j1 are not. Thus in the presence of cross diffusion, n positive but unfavorable values of f can yield negative values of nþ1 even with an arbitrarily small timestep. f In the radiation belts, f is usually largest for large a0 and small energy, and the cross diffusion coefficient is typically positive, ^ or ða0 ; pÞ ^ grid, 2Cf iþ1j1 is likely to be negative so on an ða0 ; EÞ and can be dominant. Realistic diffusion coefficients range over several orders of magnitude and Da0 p can abruptly change sign (Albert, 2005), which seems unlikely to help the situation. 6.1. Other algorithms nþ1
n
Using f instead of f in the discretized diffusion operators gives a fully implicit scheme, which (with constant D’s) is linearly stable for any size timestep. Though this is very inefficient to implement except in one dimension, so-called alternating direction implicit (ADI) schemes are much more practical and retain much of the theoretical stability. Variants of ADI which treat cross diffusion have been formulated by Beam and Warming (1980) and Gourlay and McKee (1977), but they have also been found vulnerable to unphysical, negative values. The explicit scheme could be ‘‘fixed’’ with a Lax approach, n nþ1 n which replaces f ij in the time difference ðf ij f ij Þ=Dt with its average over neighboring grid points. This has the effect of deliberately adding grid diffusion to the physical diffusion (Press et al., 1992). This may be reasonable in a problem dominated by other processes, but it is not appropriate when the physical diffusion is the focus of interest. (Since the grid diffusion would have to be comparable to the maximum value of the cross diffusion, smaller values would be completely masked.) Other promising numerical approaches include the supportoperators method of Shashkov and Steinberg (1996), which is optimized for rough diffusion coefficients; the method of lines (Sommeijer et al., 1997) which discretizes the diffusion equation in space but not in time; and a wide variety of finite volume and finite element methods. Tao et al. (2008) reformulated Eq. (2) as a Fokker–Planck equation, which they treated with random walk techniques; this is straightforward, but relatively inefficient. Albert and Young (2005) attacked the problematic cross terms directly, by constructing new variables such that cross terms are absent from the transformed equation. The resulting equation is essentially of the same form as Eq. (1) or (2) with the cross terms simply omitted, but the new variables retain the full physical
ARTICLE IN PRESS 1668
J.M. Albert / Journal of Atmospheric and Solar-Terrestrial Physics 71 (2009) 1664–1668
effects. This equation can be treated with simple, well-known finite difference algorithms that are robust and reliable.
7. Summary The significance of pitch angle-energy cross diffusion has been discussed from several points of view. Normalizing the diffusion equation, written in either ða0 ; pÞ or ðJ1 ; J 2 Þ, shows that the diffusion operators for all the variables are present on an equal footing. For single resonant waves, the diffusion coefficients themselves have well-defined ratios Dap =Daa and D2ap =ðDaa Dpp Þ, which are inseparable from the fundamental theoretical concepts of single wave characteristics and resonant diffusion curves. Integrating over the wave distribution, summing over resonances, and bounce averaging modify these values, ratios, and the concept itself of a diffusion path in momentum space (it becomes two-dimensional rather than one-dimensional). A simple two-dimensional model diffusion equation with constant coefficients was written, and its Green function given and illustrated. Its contours are ellipses, whose size, eccentricty, and orientation depend on all three diffusion coefficients. It was demonstrated that the cross term can change the value at a given location by a factor of 10. This equation was then combined with initial and boundary conditions modeled on actual pitch angle and energy diffusion, and solved numerically with a simple, standard algorithm. The solutions again show the effect of the cross terms, but also yield negative values of f , hinting at the bad behavior found when the same algorithm is applied to realistic diffusion coefficients. The problem was traced to the finite differencing of the cross terms; the textbook approach (Richtmyer and Morton, 1967), while numerically stable (for constant coefficients), is inherently unsatisfactory. However, as argued above, these terms are physically meaningful and should be included somehow. Fortunately, there are now several alternative means available to do so.
Acknowledgments This work was supported by the Space Vehicles Directorate of the Air Force Research Laboratory.
References Abel, B., Thorne, R.M., 1998. Electron scattering loss in Earth’s inner magnetosphere, 1, dominant physical processes. Journal of Geophysical Research 103, 2385. Albert, J.M., 2004. Using quasi-linear diffusion to model acceleration and loss from wave–particle interactions. Space Weather 2, S09S03. Albert, J.M., 2005. Evaluation of quasi-linear diffusion coefficients for whistler mode waves in a low density plasma. Journal of Geophysical Research 110, A03218. Albert, J.M., Young, S.L., 2005. Multidimensional quasi-linear diffusion of radiation belt electrons. Geophysical Research Letters 32, L14110. Beam, R.M., Warming, R.F., 1980. Alternating direct implicit methods for parabolic equations with a mixed derivative. SIAM Journal on Scientific and Statistical Computing 1, 131. Gourlay, A.R., McKee, S., 1977. The construction of hopscotch methods for parabolic and elliptic equations in two space dimensions with a mixed derivative. Journal of Computational and Applied Mathematics 3, 201. Horne, R.B., Glauert, S.A., Thorne, R.M., 2003. Resonant diffusion of radiation belt electrons by whistler-mode chorus. Geophysical Research Letters 30 (9), 1493. Horne, R.B., Thorne, R.M., Glauert, S.A., Albert, J.M., Meredith, N.P., Anderson, R.R., 2005. Timescale for radiation belt electron acceleration by whistler mode chorus waves. Journal of Geophysical Research 110, A03225. Lorentzen, K.R., Blake, J.B., Inan, U.S., Bortnik, J., 2001. Observations of relativistic electron microbursts in association with VLF chorus. Journal of Geophysical Research 106, 6017. Lyons, L.R., 1974. Pitch angle and energy diffusion coefficients from resonant interactions with ion-cyclotron and whistler waves. Journal of Plasma Physics 12, 417. Lyons, L.R., Thorne, R.M., Kennel, C.F., 1972. Pitch-angle diffusion of radiation belt electrons within the plasmasphere. Journal of Geophysical Research 77, 3455. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes. Cambridge University Press, New York. Richtmyer, R.D., Morton, K.W., 1967. Difference Methods for Initial-Value Problems. Wiley, New York. Schulz, M., 1991. The magnetosphere. In: Jacobs, J.A. (Ed.), Geomagnetism, vol. 4. Academic, San Diego, CA, pp. 87–293. Schulz, M., Lanzerotti, L.J., 1974. Particle Diffusion in the Radiation Belts. Springer, New York. Shashkov, M., Steinberg, S., 1996. Solving diffusion equations with rough coefficients in rough grids. Journal of Computational Physics 129, 383. Sommeijer, B.P., Shampine, L.F., Verwer, J.G., 1997. RKC: an explicit solver for parabolic PDEs. Journal of Computational and Applied Mathematics 88, 315. Summers, D., Thorne, R.M., Xiao, F., 1998. Relativistic theory of wave-particle resonant diffusion with application to electron acceleration in the magnetosphere. Journal of Geophysical Research 103(A9) (20), 487. Tao, X., Chan, A.A., Albert, J.M., Miller, J.A., 2008 Stochastic modeling of multidimensional diffusion in the radiation belts. Journal of Geophysical Research 113, A07212. Thorne, R.M., O’Brien, T.P., Shprits, Y.Y., Summers, D., Horne, R.B., 2005. Timescale for MeV electron microburst loss during geomagnetic storms. Journal of Geophysical Research 110, A09202. Varotsou, A., Boscher, D., Bourdarie, S., Horne, R.B., Glauert, S.A., Meredith, N.P., 2005. Simulation of the outer radiation belt electrons near geosynchronous orbit including both radial diffusion and resonant interaction with whistlermode chorus waves. Geophysical Research Letters 32, L19106. Walker, A.D.M., 1993. Plasma Waves in the Magnetosphere. Springer, New York.