The Hamiltonian Mechanics of Stochastic Acceleration
aa r X i v : . [ m a t h - ph ] D ec The Hamiltonian Mechanics of Stochastic Acceleration
J. W. Burby
Princeton Plasma Physics Laboratory,Princeton, New Jersey 08543, USA
A. I. Zhmoginov
Department of Physics, University of California, Berkeley, California 94720, USA
H. Qin
Princeton Plasma Physics Laboratory,Princeton, New Jersey 08543, USA andDept. of Modern Physics, University of Science andTechnology of China, Hefei, Anhui 230026, China (Dated: September 28, 2018)
Abstract
We show how to find the physical Langevin equation describing the trajectories of particles un-dergoing collisionless stochastic acceleration. These stochastic differential equations retain not onlyone-, but two-particle statistics, and inherit the Hamiltonian nature of the underlying microscopicequations. This opens the door to using stochastic variational integrators to perform simulationsof stochastic interactions such as Fermi acceleration. We illustrate the theory by applying it totwo example problems. ntroduction . — The term “stochastic acceleration” refers to the chaotic motion of par-ticles subjected to a prescribed random force. Such motion occurs in myriad contexts; theturbulent electromagnetic fields present in the interstellar medium and the RF wave fieldsfound in magnetic fusion devices are just two examples. In the astrophysical context, it isthought to be partially responsible for the presence of cosmic rays in our solar system [1].In the magnetic fusion context, it might explain the presence of certain high-energy tailsobserved in the National Spherical Torus Experiment when neutral beams are fired intoRF-heated plasmas [2].Robust modeling of stochastic acceleration requires statistical approaches. The dominantapproach is to employ the Fokker-Planck equation [3–7] for the one-particle distributionfunction. However, when studying Richardson dispersion [8, 9], and more generally anyphenomenon governed by the two-particle distribution function [10], the one-particle Fokker-Planck equation is insufficient. This is because spatial correlations in the random force fieldprevent the two-particle distribution function from factoring as a product of one-particledistribution functions. A superior statistical model when multi-particle statistics are inquestion would be a Langevin equation for particle trajectories. A wisely-chosen Langevinequation could capture the physics of the one- and two-particle distribution functions whileproviding an attractive means to perform Monte Carlo simulations of stochastic acceleration.Currently, there are no satisfactory methods for finding such a Langevin equation.The purpose of this Letter is to describe, for the first time, a systematic procedure forpassing from a microscopic description of stochastic acceleration in terms of Hamiltonianequations of motion to the physically-correct Langevin equation for particle trajectories inthe long-time limit. We will also show that, aside from reproducing the correct multi-particlestatistics, this Langevin equation inherits the Hamiltonian structure of the microscopic dy-namics. Specifically, we will show that the Langevin equation is a Hamiltonian stochasticdifferential equation (SDE) [11]. Thus, this work proves that symmetries of the macroscopicphysical laws governing stochastic acceleration lead to conservation laws.We will focus our attention on stochastic acceleration problems similar to those studiedin [3–6]. These consist of a collection of non-interacting particles moving through a pre-scribed Hamiltonian force field. By assumption, the force will consist of a small-amplitudeperturbation superimposed over a time-independent background. The perturbed force feltby a particle will be assumed to have a correlation time much shorter than any bounce time2ssociated with the perturbation, zero mean, and temporally homogeneous statistics. Theseassumptions preclude treating Coulomb collisions because the polarization field producedby a particle cannot be modeled as a prescribed field; the polarization force depends on thehistory of a particle’s orbit. They also preclude the treatment of strong turbulence [12]. The main idea. — Mathematically, this type of problem can be described as follows.Each particle moves through a 2 n -dimensional single-particle phase space M according to adynamical law given by a time-dependent vector field X t ; if z t ∈ M denotes the trajectoryof a particle in M , then ˙ z t = X t ( z t ) . (1)Because the only forces present are Hamiltonian, X t must be Hamiltonian in the sensethat there is some Poisson bracket {· , ·} and some time-dependent Hamiltonian, H t , suchthat ˙ z i = { z i , H t } , where z i denotes an arbitrary coordinate system on M [13]. By standardmathematical convention, this is written X t = X H t [14]. The presumed form of the force thenimplies H t = H + ǫh t , where ǫ ≪ H describes the mean time-independent background,and h t describes the small-amplitude random perturbation. Moreover, X h t evaluated ona particle trajectory must have a correlation time τ ac much shorter than some constant τ ,which, in turn, is much shorter than any bounce time associated with the perturbation τ b , τ ac ≪ τ ≪ τ b .Our goal in this Letter is to find the correct coarse-grained version of the microscopicequations of motion, X H t . Specifically, we seek a Langevin equation in the form δz t = X ( z t ) d t + X k ≥ X k ( z t ) δW kt (2)whose solutions correctly reproduce the late-time statistical behavior of solutions to the mi-croscopic equations of motion. Here X k are vector fields on M that must be determined, W k are independent ordinary Wiener processes, and δ denotes the Stratonovich differen-tial [15] (sometimes also written ◦ d). We will identify the X k by demanding that Eq. (2)possess two properties: it must generate the Fokker-Planck equations for the one- and two-particle distribution functions, f t ( z ) and g t ( z , z ). The two-particle distribution functionis defined such that the probability particle 1 is in the region U ⊂ M and particle 2 is inthe region U ⊂ M at time t is given by R U R U g t d z d z , where d z denotes the Liouvillemeasure [14]. Baxendale [16] has proven that a Langevin equation is uniquely determined3y its one- and two- particle Fokker-Planck equations. Therefore, these conditions uniquelyspecify the Langevin equation we seek. In particular, the requirement that two-particlestatistics be accurately reproduced is critical; Baxendale’s work implies that constrainingthe Langevin equation only to be consistent with the one-particle Fokker-Planck equationwould not identify it uniquely.Physically, the reason that the two-particle Fokker-Planck equation contains more infor-mation than the one-particle Fokker-Planck equation can be understood as follows. After ashort amount of time ∆ t , the displacement of a particle initially located at z at time t isgiven approximately by ∆ t X t ( z ). Similarly, the displacement of a particle initially locatedat z is nearly ∆ t X t ( z ). Because the random force field generally has spatial correlations, X t ( z ) and X t ( z ) are not statistically independent. Thus, the probability distribution of( z ′ , z ′ ), where z ′ i ≈ z i + ∆ t X t ( z i ), will not be given by the product of the distribution of z ′ with that of z ′ . This failure-to-factor precludes determining the two-particle distributionfunction from the mere knowledge of the one-particle distribution function. Note that thisis true in spite of the fact that these particles do not interact ; because the random force isassumed to be prescribed, the time-evolution of z is decoupled from the time-evolution of z . Identifying the Langevin equation. — The one-particle Fokker-Planck equation associatedwith Eq. (2) is given by [15, 16] ∂f t ∂t = − div( f t X ) + 12 X k ≥ div(div( f t X k ) X k )= A f t , (3)while the two-particle Fokker-Planck equation [16–18] is given by ∂g t ∂t = A (1)1 g t + A (2)1 g t + X k ≥ div (1) div (2) : g t X k ( z ) ⊗ X k ( z ) . (4)The divergence operators in these expressions are defined relative to the Liouville volumeform and the colon indicates the full contraction of second-rank tensors, a : b ≡ a ij b ij .Because these equations follow from Eq. (2) via rigorous mathematics, we will refer to themas the mathematical Fokker-Planck equations.On the other hand, under our assumption that the correlation time of the perturbedforce is much shorter than a bounce time, standard coarse-graining procedures [19, 20]4ogether with a decomposition theorem for time-ordered exponentials [21] lead to the late-time evolution laws for the one- and two-particle distribution functions associated with themicroscopic equations of motion, Eq. (1). The physical one-particle Fokker-Planck equationis given by ∂f t ∂t = − (cid:26) f t , H + ǫ τ E [ s ] (cid:27) + ǫ τ E [ {{ f t , s } , s } ]= A f t , (5)while the physical two-particle Fokker-Planck equation (see the supplementary material fora derivation) is given by ∂g t ∂t = A (1)1 g t + A (2)1 g t + ǫ τ E [ α : d (1) d (2) g t ] . (6)The notation introduced in these two equations is defined as follows: E denotes an expecta-tion value; the functions s , s are defined by s = Z τ exp( λX H ) ∗ h τ − λ d λ (7a) s = 12 Z τ Z a { exp( bX H ) ∗ h τ − b , exp( aX H ) ∗ h τ − a } d b d a ; (7b)exp( Y ) : M → M denotes the time-one advance map of the dynamical system definedby the vector field Y ; (exp( Y ) ∗ h )( z ) ≡ h (exp( − Y )( z )); the superscripts indicate whichargument of g t that A and the exterior derivative d should be applied to; and α ( z , z ) ≡ E [ X s ( z ) ⊗ X s ( z )] is the two-point covariance tensor.The X k must be chosen so that the mathematical Fokker-Planck equations, Eqs. (3) and(4), are equivalent to the physical Fokker-Planck equations, Eqs. (5) and (6). However, adirect comparison of these two pairs of equations is difficult with Eqs. (5) and (6) in theircurrent form. To eliminate this issue, we will obtain a special decomposition of the two-pointcovariance tensor α ( z , z ).As a first step, notice that if we fix a one-form ξ ∈ T ∗ z M , then we can define a vectorfield Y ξ on M by contracting ξ with α on the left according to Y ξ ( z ) = α ( z , z )( ξ, · )= E [ ξ ( X s ( z )) X s ( z )] . (8)5y forming all possible linear combinations of vector fields of this form, we can construct a(potentially infinite dimensional) linear space of vector fields [22, 23], which we will denote H , H = { linear combinations of Y ξ , ξ ∈ T ∗ M } . (9)Because each Y ξ is of the form Y ξ ( z ) = X ¯ H ( z ) with ¯ H ( z ) = E [ ξ ( X s ( z o )) s ( z )], and the sumof Hamiltonian vector fields is again Hamiltonian, H consists entirely of Hamiltonian vectorfields. Moreover, following Baxendale [16, 23], we see that H is a real Hilbert space whoseinner product is defined by the formula h Y ξ , Y η i H = α ( z , z )( ξ, η )= E [ ξ ( X s ( z )) η ( X s ( z ))] , (10)where ξ ∈ T ∗ z M and η ∈ T ∗ z M . Therefore we may choose an orthonormal basis { e k } k ≥ for H , where each e k must be of the form e k = X H k . A simple calculation then leads to thedesired decomposition of α : α ( z , z ) = X k ≥ X H k ( z ) ⊗ X H k ( z ) . (11)Using this decomposition of the two-point covariance tensor, it is straightforward tomanipulate Eqs. (5) and (6) into the same form as Eqs. (3) and (4). After doing so, it istrivial to identify the correct X k . Indeed, we have found that the physical Langevin equationis given by δz t = X ˜ H ( z t ) d t + X k ≥ X ˜ H k ( z t ) δW kt , (12)where ˜ H = H + ǫ τ E [ s ] , ˜ H k = ǫ √ τ H k (13)Recall that the X H k are defined to be an orthonormal basis of the Hilbert space H definedin Eq. (9). Also recall that all of the above manipulations have been performed under theassumption that the correlation time of the perturbed force felt by a particle is much shorterthan any bounce time associated with the perturbation.Because the coefficients in the Langevin equation for stochastic acceleration, Eq. (12),are all Hamiltonian vector fields, this equation is an example of a stochastic Hamiltonian ystem , the foundations of which are developed in [11]. It is in this sense that the Langevinequation for stochastic acceleration inherits the Hamiltonian structure of the microscopicequations. In particular, SDEs of this type are known to arise from a stochastic variationalprinciple for which Noether’s theorem applies. Thus, even at the dissipative macroscopiclevel, symmetries imply the presence of conservation laws. Example 1. — We will find the physical Langevin equation for two example stochasticacceleration problems. Generally speaking, finding the coefficients of the physical Langevinequation involves finding an orthonormal basis for the space H , a task which may be ana-lytically intractable. But, by Mercer’s theorem [24], this task can be cast as an eigenvalueproblem for which there are existing numerical solution methods. In any case, in theseexamples, the analytical route is tractable.First, consider a single-species, unmagnetized plasma subjected to a random weak elec-trostatic pulse at τ -second intervals. Assume that the pulses are uniform in space andconstant in magnitude, but uniformly and independently distributed in direction. Thus, the k ’th pulse is generated by a potential of the form φ k ( x , t ) = ( z k · x ) φ o u ( t − kτ ), where z k is arandom vector uniformly distributed over the unit sphere and u ( t ) is a temporal windowingfunction localized at t = τ / τ , we must (a) calculate s and s using Eqs. (7a) and (7b), (b) find an orthonor-mal basis { X H k } k ≥ for the space H defined in Eq. (9), and (c) write down Eq. (12) with ˜ H and ˜ H k calculated using Eq. (13). The results of these three steps are as follows.(a) A quick calculation shows that s = m o z · x − m z · v (14a) s = const (14b)where m o = ( q/m ) φ o R τ u ( s )d s , m = ( q/m ) φ o R τ ( τ − s ) u ( s )d s , and q/m is the charge-to-mass ratio.(b) Each Y ξ must be of the form Y ξ = X g βγ , where g βγ ( x , v ) = 13 ( m β + m o γ ) · ( m v − m o x ) , (15)and β , γ are arbitrary constant 3-component vectors. Using this expression, it is simple to7nd an orthonormal basis for H . One is given by { X ¯ H k } k =1 .. , with H i ( x , v ) = 1 √ e i · ( m v − m o x ) , (16)where { e i } i =1 .. is the standard basis for R .(c) Finally, the physical Langevin equation is given by δx i = v i d t + 1 √ τ m δW i (17a) δv i = 1 √ τ m o δW i , (17b)where i = 1 , , ∂f t ∂t + v · ∇ f t = 16 τ ( m ∇ f t + m o m ∇ · ∇ v f t + m o m ∇ v · ∇ f t + m o ∇ v f t ) . (18)On the other hand, given an arbitrary function φ ( x , v ), the SDE δx i = v i dt + m √ τ (cid:0) cos( φ ) δW ,i − sin( φ ) δW ,i (cid:1) (19a) δv i = m o √ τ (cid:0) cos( φ ) δW ,i − sin( φ ) δW ,i (cid:1) , (19b)where the W ,i , W ,j are six independent ordinary Wiener processes, will also generateEq. (18). However, when φ is not constant, the two-particle Fokker-Planck equation gener-ated by Eq. (19) will differ from the two-point Fokker-Planck equation generated by Eq. (17).This can be verified using Eq. (6). The procedure identified here selects φ = 0 as the phys-ical choice. In particular, it shows that a Langevin equation with the correct one-particleFokker-Planck equation may still incorrectly reproduce the two-particle distribution func-tion.The inadequacy of Eq. (19) can also be understood intuitively as follows. Chaotic motionsof any two particles experiencing the electrostatic pulses are “synchronized” since the pulsesare independent of x and v . The Langevin equation (19), on the other hand, desynchronizesparticle trajectories by involving additional Wiener processes, in spite of giving the correctone-particle Fokker-Planck equation. Example 2. — Next, consider a minority population of magnetized fast ions movingthrough a plane lower-hybrid wave that propagates perpendicular to the magnetic field.8ssume the wave has a high harmonic number and a wavelength small compared to a typi-cal ion gyroradius. Karney [25] has shown that the dynamics of the perpendicular velocityof these ions are governed by a canonical time-dependent Hamiltonian system with Hamil-tonian H t = I − ǫ sin( √ I sin θ − νt ) , (20)where I is the normalized magnetic moment, t the time normalized by the gyroperiod, θ the gyrophase, ν the harmonic number, and ǫ the normalized wave amplitude. Moreover,when ǫ exceeds a threshold value, an ion’s motion becomes chaotic. This chaotic motioncomes as the result of the effective randomization of the wave phase felt by an ion after agyroperiod. Thus, above the threshold for chaos, we can model the wave phase as beingrandomized every gyroperiod by a random variable η . That is, we can replace the exact chaotic ion motion with a stochastic approximation; see Ref. [26] for Chirikov’s applicationof the same modeling approach to the standard map. This allows us to apply the formalismdeveloped in this Letter to find the physical Langevin equation describing the stochasticparticle trajectories at times much longer than the gyroperiod.As in the previous example, the first step is to calculate s and s . Set τ = 2 π and adoptthe rough approximation ∞ X n = −∞ J n ν − n exp( inθ ) ≈ J n o δ exp( in o θ ) , (21)where ν = n o + δ , | δ | < , and J n = J n ( √ I ) denotes the Bessel function of the first kind[27]. This approximation amounts to selecting the most slowly varying term in the sum inEq. (21). Then, upon directly evaluating the integrals in Eqs. (7a) and (7b), the resultingexpressions for s and E [ s ] are s = 2 π sinc( πδ ) J n o sin( n o θ + η ) (22a) E [ s ] = π ∞ X m = −∞ J m +1 − J m − m − ν + π πδ ) J n o +1 − J n o − δ , (22b)where η is a random variable uniformly distributed over the interval [0 , π ] and sinc( x ) =sin( x ) /x . 9ext, the space H can be constructed using the above expression for s . In this case, H is two-dimensional and has a basis { X H , X H } , where H ( I, θ ) = √ π sinc( πδ ) J n o ( √ I ) cos( n o θ ) (23a) H ( I, θ ) = √ π sinc( πδ ) J n o ( √ I ) sin( n o θ ) . (23b)Finally, the coefficients for the Langevin equation, Eq. (12), can be derived using Eq. (13).The result is δI = ǫ √ π sinc( πδ ) n o J n o ( √ I ) × (cid:0) sin( n o θ ) δW − cos( n o θ ) δW (cid:1) (24a) δθ = (cid:18) ǫ π ∂∂I E [ s ] (cid:19) d t + (cid:18) ǫ r π I sinc( πδ ) J ′ n o ( √ I ) × (cid:0) cos( n o θ ) δW + sin( n o θ ) δW (cid:1) (cid:19) . (24b)The diffusion of the magnetic moment I predicted by Eq. (24) has already been studied byKarney [25]. However, Eq. (24) extends and compliments Karney’s results by predicting theappropriate diffusion in gyrophase, as well as the correct two-particle statistics. Concluding remarks. — We have shown how to derive the physical Langevin equationfor particle trajectories undergoing stochastic acceleration. This SDE correctly generatesthe correct one- and two-particle Fokker-Planck equations and inherits the Hamiltonianstructure of the microscopic equations of motion. This inheritance is theoretically satisfyingbecause it is a direct consequence of demanding consistency with the physical one- and two-particle Fokker-Planck equations. It also implies that symmetries of the macroscopic physicallaws governing stochastic acceleration imply the presence of conservation laws. While thisrelationship is well known at the microscopic level, it is a pleasant surprise that it remainsintact upon passing to dissipative macroscopic equations.A Hamiltonian Langevin equation [11] is a Stratonovich SDE of the form given in Eq. (12).If a loop of initial conditions for this SDE evolves under a given realization of the noise,then the action of that loop is constant in time. In addition, these equations arise froma stochastic action principle [11] for which Noether’s theorem applies. Thus, by showingthe physical Langevin equation is Hamiltonian, we have also identified potentially powerful10ools for the analysis of stochastic acceleration. In particular, using the methods of Bou-Rabee [28], the stochastic action principle can be used to develop variational integratorsfor Eq. (12). Because these integrators are known to possess superior long-term statisticalfidelity [29], this approach may prove to be useful in Monte Carlo simulations of stochasticacceleration.The authors would like to express their appreciation to I. Dodin, J. A. Krommes, J.Parker, and G. W. Hammett. This work was supported by DOE contracts DE-AC02-09CH11466 and DE-FG02-04ER41289. [1] E. Fermi, Phys. Rev. , 1169 (1949).[2] D. Liu, W. W. Heidbrink, M. Podest`a, R. E. Bell, E. D. Fredrickson, S. S. Medley, R. W.Harvey, and E. Ruskov, Plasma Phys. Controlled Fusion , 025006 (2009).[3] P. A. Sturrock, Phys. Rev. , 186 (1966).[4] D. E. Hall and P. A. Sturrock, Phys. Fluids , 2620 (1967).[5] D. D. Barbosa, Astrophys. J. , 383 (1979).[6] V. Petrosian and S. Liu, Astrophys. J. , 550 (2004).[7] R. J. Hamilton and V. Petrosian, Astrophys. J. , 350 (1992).[8] L. F. Richardson, Proc. R. Soc. London, Ser. A , 709 (1926).[9] M.-C. Jullien, J. Paret, and P. Tabeling, Phys. Rev. Lett. , 2872 (1999).[10] A. K. Mukhopadhyay and J. Goree, Phys. Rev. Lett. , 165003 (2012).[11] J. L´azaro-Cam´ı and J. Ortega, Rep. Math. Phys. , 65 (2008).[12] D. F. DuBois and M. Espedal, Plasma Phys. , 1209 (1978).[13] C. Grebogi, A. N. Kaufman, and R. G. Littlejohn, Phys. Rev. Lett. , 1668 (1979).[14] R. Abraham and J. Marsden, Foundations of Mechanics , AMS Chelsea publishing (AmericanMathematical Soc., 2008).[15] C. Gardiner,
Stochastic Methods: A Handbook for the Natural and Social Sciences , SpringerSeries in Synergetics (Springer, 2009).[16] P. Baxendale, Compositio Math. , 19 (1984).[17] B. Schmalfuss, Dyn .Syst. , 303 (2001).[18] H. Kunita, Lectures on stochastic flows and applications , Tata Institute Lectures on Mathe- atics and Physics (Springer, 1987).[19] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications , LectureNotes in Mathematics (Springer-Verlag, 1996).[20] M. Bazant, “18.366 random walks and diffusion, fall 2006,” (MIT OpenCourseWare:Massachusetts Institute of Technology), http://ocw.mit.edu/courses/mathematics/18-366-random-walks-and-diffusion-fall-2006 (see Lecture 8 under Study Materials).[21] C. S. Lam, J. Math. Phys. , 5543 (1998).[22] N. Aronszajn, Trans. Amer. Math. Soc. , 337 (1950).[23] P. Baxendale, Amer. J. Math. , 891 (1976).[24] J. Mercer, Philos. Trans. Roy. Soc. London Ser. A , 415 (1909).[25] C. F. F. Karney, Phys. Fluids , 2188 (1979).[26] B. V. Chirikov, Phys. Rep. , 263 (1979).[27] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions , Applied mathematicsseries (Dover Publications, Incorporated, 1964).[28] N. Bou-Rabee and H. Owhadi, IMA J. Numer. Anal. , 421 (2009).[29] N. Bou-Rabee and H. Owhadi, SIAM J. Numer. Anal. , 278 (2010)., 278 (2010).