Adaptive time-stepping Monte Carlo integration of Coulomb collisions

0.5 MB·

by  Konsta Särkimäki

Adaptive time-stepping Monte Carlo integration of Coulomb collisions
K.Särkimäkia,∗,E.Hirvijokib,J.Teräväa
aAaltoUniversity,Espoo,Finland
bPrincetonPlasmaPhysicsLaboratory,Princeton,NJ,USA
7
1 Abstract
0
WereportanaccessibleandrobusttoolforevaluatingtheeffectsofCoulombcollisionsonatestparticleinaplasma
2
thatobeysMaxwell-Jüttnerstatistics. TheimplementationisbasedontheBeliaev-Budkercollisionintegralwhichal-
n
lowsboththetestparticleandthebackgroundplasmatoberelativistic.Theintegrationmethodsupportsadaptivetime
a
J stepping,whichisshowntogreatlyimprovethecomputationalefficiency. TheMonteCarlomethodisimplemented

8 forboththethree-dimensionalparticlemomentumspaceandthefive-dimensionalguidingcenterphasespace.
1 Detaileddescriptionisprovidedforboththephysicsandimplementationoftheoperator. Thefocusisinadaptive
integration of stochastic differential equations, which is an overlooked aspect among existing Monte Carlo imple-
]
h mentations of Coulomb collision operators. We verify that our operator converges to known analytical results and
p demonstratethatcarelessimplementationoftheadaptivetimestepcanleadtoseverelyerroneousresults.
- Theoperator isprovidedas aself-containedFortran 95moduleandcan beincludedinto existingorbit-following
m
toolsthattraceeitherthefullLarmormotionortheguidingcenterdynamics. Theadaptivetime-steppingalgorithm
s is expected to be useful in situations where the collision frequencies vary greatly over the course of a simulation.
a
l Examples include the slowing-down of fusion products or other fast ions, and the Dreicer generation of runaway
p
electronsaswellasthegenerationoffastionsorelectronswithionorelectroncyclotronresonanceheating.
.
s
c Keywords: Coulombcollision,MonteCarlo,Fokker-Planckequation,Milsteinmethod
i
s
y
h PROGRAMSUMMARY Fokker-Planck coefficients, a feature useful for constructing
p
Program Title: AMCC – (A)daptive (M)onte-(C)arlo accurate orbit averaged collision operators. The package
[
(C)oulombcollisions also provides explicit one-step sympletic integrator for

1 Licensingprovisions:LGPL-2.1 the relativistic Lorentz force that can be used for tracing
v Programminglanguage:Fortran95 test-particlesingivenelectromagneticbackgrounds.
3 Nature of problem: Test-particle tracing is a common feat
4 within existing fusion applications. While efficient adaptive
0
methodsexistforintegratingtheincompressibleHamiltonian
5 flow, the effects of Coulomb collisions are commonly 1. Introduction
0
implementedwithfarlesssophisticatedalgorithms.
. Standard Runge-Kutta methods developed for ordi-
1 Solution method: The relativistic Fokker-Planck equation nary differential equations (ODEs), are based on first
0 for test-particles in Maxwell-Jüttner background plasmas
7 is converted into a stochastic differential equation. The convertingtheODEintoanintegralequation,andthen
1 stochastic differential equation is solved using adaptive discretizing the integral. Stochastic differential equa-
v: Monte Carlo techniques. Methods to evaluate the effect of tions (SDEs), often encountered when stochastic pro-
i Coulomb collisions for both the three-dimensional particle cessessuchascollisionsareofinterest,canequivalently
X
momentum space and the five-dimensional reduced guiding betransformedintointegralequations, withtheexcep-
r centerphasespaceareincluded. tion that the resulting integral is non-Riemannian [1].
a
Additional comments including Restrictions and Unusual The discretization rules used for ODEs then no longer
features: apply and, instead, either Itô or Stratonovich calculus
The package includes optionality to evaluate the relativistic mustbeadoptedtoobtainnumericalmethodsforsolv-
ingSDEs.
∗Correspondingauthor. Perhaps the most straight-forward discretization
E-mailaddress:konsta.sarkimaki@aalto.fi method for SDEs is the Euler-Maruyama method.
PreprintsubmittedtoComputerPhysicsCommunications January19,2017
While simple to implement, the Euler-Maruyama in section 5 for both particle and guiding center pic-
method suffers from constant time-step requirement in ture. The collision operators are verified in section 6,
the sense that one is not allowed to simply recompute where we also compare the adaptive method with the
the step if the estimated integration error turns out too fixed time step scheme, and confirm that the adaptive
large,asthismightleadtoaconvergencetoawrongso- method isboth faster andmore accurate thanthe com-
lution[2].ThereasonisthatwhileEuler-Maruyamahas monlyusedEuler-Maruyamamethod.
weakorderofconvergence1.0,itsstrongorderofcon-
vergenceisonly0.5. ForadaptiveintegrationofSDEs,
2. Test-particleFokker-Planckequation
a method must have strong convergence of at least 1.0
to guarantee that it converges to a correct solution. A In the limit of binary collisions, the collisional evo-
common application, where time steps are tested and lution of the particle distribution function of species a,
rejectedwhenneeded,istest-particletracingwherethe interactingwithspeciesb,isdeterminedby
Hamiltonianmotionoftheparticleissolvedwithadap-
∂f (cid:88)
tiveRunge-Kuttamethods.Therefore,totreatcollisions a = C [f , f ]. (1)
∂t ab a b
withamannerconsistentwiththeHamiltonianmotion, b
a higher order scheme than the Euler-Maruyama must
In a plasma, the collisions are dominated by the small
beused.
anglescatteringevents,sothatthecollisionoperatorcan
Further, the test-particle collision frequencies in a
bewrittenintheLandau-Fokker-Planckform[6]
plasma can vary greatly during the course of follow-
ing the particle trajectory: energetic particles start as Γ ∂ (cid:90) (cid:32) f¯ ∂f f ∂f¯ (cid:33)
C [f , f ]= ab · du¯ U · b a − a b ,
nearly collisionless and slowly drift toward the more ab a b 2ma∂u R3 BB ma ∂u mb ∂u¯
collisional bulk population, while some particles from (2)
the bulk may accelerate to high energies. Therefore, whereU istherelativisticBeliaev-Budkertensor.The
BB
adaptivetime-steppingmethodsshouldbeadoptedalso coordinateu = p/m cisparticleamomentumnormal-
a
forthestochasticcontributiontotheparticlemotion,to ized to the rest mass m and speed of light c, while
a
reduce both the convergence errors and computational u¯ = p¯/m c denotes the same for particle of species
b
costs. b. The quantities with an overbar are evaluated at u¯.
We thus see fit to introduce an adaptive time- The operator also depends on the species charges, q
a
stepping algorithm that would allow incorporating the andq ,throughthecoefficientΓ =q2q2lnΛ /(4πε2)
b ab a b ab 0
Coulombscatteringintotracingofeithertheparticleor where ε is the vacuum permittivity, and lnΛ is the
0
guiding-centerdynamics. OurMonteCarloimplemen- Coulomblogarithm. TheCoulomblogarithmdescribes
tationevaluatesthemomentumchangeofatest-particle theratioofminimumandmaximumimpactparameters
when it collides with a Maxwell-Jüttner background lnΛ = lnr /r and as such it indicates by which
max min
plasma, consisting of electrons and possibly multiple factor small angle scattering dominates the large angle
ion species. The operator is based on the Beliaev- scattering. Naturalchoiceforthemaximumimpactpa-
Budker collision integral [3] and, therefore, is applica- rameteristheDebyelength, λ = (cid:113)ε (cid:80) (cid:16)T /n q2(cid:17),
bleevenifeitherorboththetestparticleandtheback- D 0 b b b b
where n is density and T is temperature. Depending
ground plasma populations are relativistic. Advanced
which one is smaller, the minimum impact parame-
MonteCarloCoulombcollisionoperatorshavebeende-
ter is determined either via classical electron radius,
velopedquiterecently[4,5],buttheoperatordeveloped
r = q q /(4πε m v˜2), where m = m m /(m + m )
herediffersfromthoseinthatitisrelativistic,adaptive, cl a b 0 r r a b a b
is the reduced mass and v˜ = (cid:104)|v −v |(cid:105) is the mean
andalsoapplicableinguidingcenterdynamics. a b
relative velocity, or by a quantum mechanical limit,
The paper is organized as follows. In section 2, we
r = (cid:126)/(2m v˜), where (cid:126) is the reduced Planck con-
begin by discussing how the Monte Carlo operator is qm r
stant. Infusionplasmas,lnΛhasvaluesintherangeof
obtainedfrom theBeliaev-Budkercollision integralby
10 – 20, but the uncertainties in estimating r means
finding the related SDE. The corresponding guiding- min
that the Coulomb logarithm is only accurate to within
centeroperatorispresentedinsection3. Discretization
1/lnΛ.
oftherelativisticparticlecollisionoperatorisfirstdone
ReturningtotheBeliaev-Budkertensor,wefindthat
with the Euler-Maruyama method, in section 4, where
ithasarathercomplicatedexpression[3]
wealsoshowhowtocomputethecollisioncoefficients
efficientlyasthesearealsoneededfortheadaptivetime- r2 (cid:16) (cid:17)
U = w2I−uu−u¯u¯ +r(uu¯ +u¯u) , (3)
stepping. The adaptive scheme is described in detail BB γ¯γw3
2
√ √
withr = γγ¯ −u·u¯, w = r2−1, andγ = 1+u2 is ThecoefficientsareillustratedinFig.1asafunctionof
theLorentzfactor. Inthenon-relativisticlimit,c → ∞, u.
the Beliaev-Budker tensor reduces to the better known Thespecialfunctionsµ (u;Θ),µ (u;Θ),andµ (u;Θ)
0 1 2
Landautensor[7] aregivenby
(cid:32) (cid:33)
UL = |u−1u¯| I− (u−|uu¯−)(uu¯|2−u¯) . (4) µ0 = γ2L0−ΘeL11/Θ+K(Θ(1−/Θγ))ue(1−γ)/Θ, (14)
2
From now on, we will focus solely on the relativis- γ2L −ΘL +(Θγ−1)ue(1−γ)/Θ
tic expressions which, of course, are valid also in the µ1 = 1 e10/ΘK (1/Θ) , (15)
non-relativisticregime. Inordertoobtainthetestparti- 2
2ΘγL +(1+2Θ2)ue(1−γ)/Θ
clecollisionoperator,wefirstwritetheBeliaev-Budker µ = 1 , (16)
collision integral, Eq. (2), in anexplicit Fokker-Planck 2 Θe1/ΘK2(1/Θ)
form wherethefunctionsL (u;Θ)andL (u;Θ)are
0 1
∂f ∂ ∂ ∂
a =− ·(K f )+ :(D f ), (5) √
∂t ∂u a a ∂u∂u a a (cid:90) u e(1− 1+s2)/Θ
cw(cid:80)iehbneDtrseabKt[hafebb][vfaebrc]etaonsrudKmDama=ba[tfi(cid:80)bo]nb,sKanoadbf[ftfuhbne]castnpiodenctaiheless-otwefnitssheoercdDoiseatffir=i-- LL01 ==(cid:90)0uddsse(1√−√11++s2s)2/Θ., ((1178))
bution f . The expressions for the species-wise coeffi- 0
b
cientsare Our notation here differs from the Ref. [8] as we
Γ (cid:90) havemultipliedboththenumeratoranddenominatorin
Dab[fb]= 2mmab2aΓ R3(cid:90)du¯ UBBf¯b, ∂f¯ ∂ (6) Ewqhse.n(Θ14i)s-sm(1a6ll).by e1/Θ to avoid floating point errors
K [f ]= a ab du¯ U · b + ·D [f ]. The Monte Carlo simulation of a given test parti-
ab b mb2m2a R3 BB ∂u¯ ∂u ab b cle distribution is based on the connection between
(7) the Fokker-Planck equation and stochastic differential
From now on, we will assume that the distributions equations. If the distribution function f (u,t) satisfies
a
f areMaxwell-Jüttnerdistributions equation (5), then an individual sample particle from
b
√ f (u,t)obeysthefollowingstochasticdifferentialequa-
f (u)= nbe− 1+u2/Θb . (8) tiaon, known as the Langevin equation, of Itô kind [9]
b 4πΘ K (1/Θ )
b 2 b
HereKν(x)(nottobeconfusedwiththecoefficientKab) du(t)= Ka(u(t),t)dt+σa(u(t),t)·dW, (19)
is the νth order modified Bessel function of the second
wheretherank-2tensorσ satisfiesthecondition
kind,andΘ =T /m c2 isthenormalizedtemperature. a
b b b
In this case, the species-wise diffusion tensor Dab and σ σ(cid:124) =2D . (20)
a a a
theforceK becomeisotropic[7]
ab
D = D (u)uˆuˆ +D (u)(I−uˆuˆ), (9) Here dW is a differential of a vector-valued, uncorre-
ab ab,(cid:107) ab,⊥
lated standard Wiener processes W ∼ N(0,tI), with
K = K (u)uˆ, (10)
ab ab N being the standard multivariate normal distribution.
where uˆ ≡ u/u is the unit vector parallel to u. The Since the diffusion tensor is diagonal, the decomposi-
coefficients K , D , and D are defined in terms tionEq.(20)iseasytoaccomplish,andwefind
ab ab,(cid:107) ab,⊥
of three special functions µ (u;Θ ), µ (u;Θ ), and
0 b 1 b (cid:112) (cid:112)
µ2(u;Θb), and they are given by the following expres- σa = 2Da,(cid:107)(u(t))uˆuˆ + 2Da,⊥(u(t))(I−uˆuˆ). (21)
sions[8]
Γ n 1 (cid:32)µ m (cid:33) 3. Guiding-centertest-particleoperator
K =− ab b 0 + aµ , (11)
ab m2c3 u2 γ m 1
a b
Γ n Θ γ In many applications the rapid oscillation of a
Dab,(cid:107) = ma2bc3b ub3 µ1, (12) charged particle in plane perpendicular to a magnetic
a field is of little interest. If the magnetic field is slowly
Γ n 1 (cid:16) (cid:17)
D = ab b u2(µ +γΘ µ )−Θ µ . (13) varying, one can resort to guiding center formalism
ab,⊥ m2ac3 2γu3 0 b 2 b 1 whichomitsthisgyromotionanddependsonlyonthe
3
a) d)
K K
0 2 u
/s Θ=10−1 /s 1
u] -1 Θ=10−1.5 u]
[ [
30 Θ=10−2 30 0
1 -2 1
-1
b) e)
D ν
150 || 4
s 3
2[u]/s 10500 og1/1012
l
0 0
c) f)
150 D⊥ 0.5 DX
2[u]/s 10500 320m/s 0.25
1
0 0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
u u
Figure1: Collisioncoefficientsasafunctionofnormalizedmomentumu = p/mcfordifferentvaluesofnormalizedtemperatureΘ = T/mc2.
(a)–(c)Thecoefficientsrelatedtotheparticlecollisionoperator,Eq.(19).(d)–(f)Thecoefficientrelatedtotheguidingcentercollisionoperator
Eqs.(43)–(44). Thevaluesareforanelectrontestparticleinteractingwithanelectronplasmawithdensityn=1020m−3. Notethelogarithmic
scalein(e).
gyro-averagedquantities,thusreducingthe6Dparticle wherethemodifiedfrictioncoefficientisdefinedby
phasespaceintoa5Dguidingcenterphasespace.
(cid:88)
Rigorous transformation of the particle phase-space Qa(u)= Qab(u)uˆ, (25)
Fokker-Planck equation (5) into guiding-center phase- b
spacewasfirstcarriedoutinRef.[10]. Thetransforma- Q (u)= K − ∂Dab,(cid:107). (26)
tion is based on the observation that the noncanonical ab ab ∂u
particle phase-space Poisson bracket could be used to
Now,usingthePoissonbrackets,oneobtains
expressthemomentum-spacederivativesaccordingto
∂f
∂· ={x, ·}, (22) ∂ta =−{xi,Qa,ifa}+{xi,Da,ij{xj, fa}}, (27)
∂u
wherethesubscriptsi, jdenotetheCartesianindicesfor
thevectorandtensorcomponents,andsummationover
where x is the spatial coordinate. This fact can be put
repeatedindicesisassumed.
to use by first observing that the isotropic test-particle
diffusiontensorsatisfies The guiding-center transformation then follows via
Lie-transform (see Ref. [10] for details), and the fi-
nal form of the collision operator for the gyroangle-
∂ ∂D
∂u ·Dab = ∂uab,(cid:107)uˆ, (23) independent guiding-center distribution function F be-
comes
(cid:32) (cid:33)
sothattheparticlephase-spacecollisionoperatorcanbe ∂F 1 ∂ ∂F
a =− · JQ F +JD · a , (28)
writteninaform ∂t J ∂Z a a a ∂Z
(cid:32) (cid:33)
∂f ∂ ∂ ∂f whereZareguidingcenterphase-spacecoordinatesand
a =− ·(Q f )+ · D · a , (24)
∂t ∂u a a ∂u a ∂u J is the transformation Jacobian. The guiding center
4
frictionanddiffusioncoefficientsare where Ω = q B/m is the gyrofrequency. In uniform
a a
magnetic field this corresponds to classical diffusion.
Q =(cid:104)∆ ·T−1Q (cid:105) (29)
a,α α a Second, the momentum magnitude and direction now
D =(cid:104)∆ ·T−1D ·∆ (cid:105), (30) haveseparatecoefficientsforthediffusion,
a,αβ α a β
where T−1 is the guiding-center push-forward, (cid:104)·(cid:105) de- D =2D , (37)
ab,u ab,(cid:107)
notesagyroaverage,and∆ aretheso-calledprojection
α D =(1−ξ2)ν , (38)
vectors,definedaccordingto ab,ξ ab
∆ ={T−1X,Z } , (31) andforthedrift
α α gc
2D
with{F,G}gc beingtheguidingcenterPoissonbracket. Kab,u =Qab,u+D(cid:48)ab,(cid:107)+ uab,(cid:107), (39)
Hereindicesα,βdenotetheguidingcentercoordinates.
Writingtheguiding-centerFokker-Planckequationina Kab,ξ =−ξνab, (40)
formsimilartoEq.(5),wefind
where
∂∂Fta =−J1 ∂∂Z ·(JKaFa)+ J1 ∂∂Z∂∂Z :(JDaFa), νab = 2Dua2b,⊥ (41)
(32) is the pitch collision frequency. Third, both ν and
ab
wherethedriftcoefficientis K diverge at u = 0, which makes the particle col-
ab,u
1 ∂ lisionoperatormoreattractivewhensimulatingthermal
Ka =Qa+ J ∂Z ·(JDa). (33) particles. Also the guiding-center specific coefficients
D ,K ,andν areillustratedinFig.1. Notethat
ab,X ab,u ab
Theguidingcenterformalismrestsontheassumption K = 0, and that we have implicitly assumed uni-
ab,X
thatthemagneticmomentisinvariant. Thisinvariance formmagneticfieldwhenwritingdownthecoefficients
can be used to reduce the number of equations of mo-
above,seeRef.[11]fordetails.
tion to four by choosing the magnetic moment as one
The guiding center collision operator we have is in
coordinate. However, thediffusiontensorisnotdiago-
curvilinear coordinates so obtaining the corresponding
nalinthiscase[11]but,aswesawinthelastsection,a
Langevinequationisnotastrivialasitwasintheparti-
diagonalbasisisdesiredwhenconsideringthenumeri-
clepicture.DetailsonthederivationoftheMonte-Carlo
calimplementation. Fortunately, thereexistsasuitable
operator are found in Appendix A, and here we only
setofcoordinateswherethediffusiontensoris(almost)
showtheresult
diagonal. This basis is Z = (X,u,ξ) where X is the
guiding center location, u is the magnitude of the nor- dZ=K dt+Σ ·dW, (42)
a a
malizedmomentum,andpitchisξ =uˆ·bˆ,wherebˆ isthe
unitvectorparalleltothemagneticfieldB(X). Inthese whereΣ isagaineasilyobtainedfromthedecomposi-
a
curvilinearcoordinatestheJacobianisJ = m Bu,and tion (1/2)Σ Σ(cid:124) = D . Written explicitly, the guiding
a a a a
thediffusiontensorbecomes centercollisionoperatorisasetofequations
(cid:88) (cid:88) (cid:88)
Da = Dab,X(IX−bˆbˆ)+ Dab,uuˆuˆ+ Dab,ξξˆξˆ, dX= (cid:112)2Da,X(I−bˆbˆ)·WX, (43)
b b b (34) du=Ka,udt+ (cid:112)2Da,(cid:107)dWu, (44)
withdiagonalmatrixI havingnon-zeroelementsonly (cid:113)
in coordinates X. TheXfriction coefficient has only one dξ =−ξνadt+ (1−ξ2)νadWξ, (45)
component
with W , W , and W being independent Wiener pro-
X u ξ
Γ n 1 (cid:32)m (cid:33) cesses.
Q =− ab b aµ . (35)
ab,u m2c3 u2 m 1
a b
4. MonteCarloalgorithm
Theguidingcentercollisionoperatorhasfewnotable
differences to the particle collision operator. First, the
The equation (19) does not have a known analytical
collisionsnowcausealsospatialdiffusion,withadiffu-
solutionand,therefore,numericalmethodsarerequired.
sioncoefficient[11]
Giveninitialconditionu(t ),thenumericalapproxima-
0
(cid:34) (cid:35)
1−ξ2 c2 tionforu(t)isobtainedbydiscretizingthetimecoordi-
Dab,X = (Dab,(cid:107)−Dab,⊥) 2 +Dab,⊥ Ω2, (36) nateastk+1 =tk+∆t,k=0,...,n,andevaluatingu(tk+1)
5
at each step using the solution of the previous step as a)
an initial condition. The evaluation is commonly done
withtheEuler-Maruyamamethodwhich,whenleaving
outthespeciessubscriptaandswitchingtoEinsteinno-
tationforbrevity,reads
(cid:88)3
ui(tk+1)=ui(tk)+Ki(u(tk))∆t+ Gij(u(tk))∆Wj (46)
j=1
wheretheindexi=1,2,3representsthecoordinatesin b)
Cartesianbasis. Inthisbasis,thematrixGreads
(cid:112) (cid:112) (cid:104) (cid:105)
G = 2D uˆ δ + 2D δ −uˆ uˆ , (47)
ij (cid:107) i ij ⊥ ij i j
whereδ istheDiracdelta.Thediscretizeddifferentials
ij
ofWienerprocessesaredrawnfromanormaldistribu-
tion∆Wi =Wi(tk+1)−Wi(tk)∼N(0,∆t).
Implementing Euler-Maruyama method is straight-
forward, and all the struggle is in evaluating the coef-
ficientsKandDthat,throughthespecialfunctionsµ0, Figure2: Integrals17and18asafunctionofuandΘ. Blacklines
µ , and µ , depend on the integrals L and L . These arecontoursoftherespectiveintegralswhereasredlineisthecurve
1 2 0 1 √
integrals(Eqs.(17)and(18))cannotbesolvedanalyti- u = Θ2+2Θ. Thevaluestabulatedforthisfigureareenoughto
covermostofthetokamakplasmasas, forexample, aslow10eV
cally and, therefore, we evaluate them via the adaptive
tungstenimpurityhasu ≈ 1×10−5,afast1GeVrunawayelectron
Simpson’smethod. TheadaptiveSimpson’smethoddi- hasu≈2×103whilecold10eVheliumhaveΘ≈3×10−9,andhot
videstheintegrationintervalintosubintervalsuntilthe 50keVelectronshaveΘ≈1×10−1.Notethattheintegralshavevery
differenceintheresultingnumericalapproximationbe- similarvalues.
tween the successive divisions is less than a given tol-
erance. However,mostofthecontributiontoL andL
0 1
5. Integrationwithanadaptivetimestep
comes from small values of u, and the adaptive Simp-
son’smethodconvergestoanincorrectvalueiftheup-
per limit for u is large while Θ is small. This can be
avoidedbynotingthatbothintegrandsequaltounityat The collision operator introduced in the last section,
u=0and,astheintegrandsdecayrapidly,theintegrals Eq. (46), was based on the Euler-Maruyama method
havepracticallyconstantvaluesbeyondacertainpoint. and, therefore, is not applicable for an adaptive time
Therefore,wecanmaketheintegrationrobustbydefin- stepscheme. Anadaptiveintegrationisdesirablesince
ingacut-offlimitforuas especiallythepitchcollisionfrequencyhasastrongde-
√ pendencyontestparticlemomentum(recallFig.1).For
(cid:112)
e(1− 1+u2c)/Θ =(cid:15) ⇒ uc = (1−Θln(cid:15))2−1, (48) example,inaproton-electronplasmawithT = 10keV,
n = 1020 m−3, a thermal electron has a pitch collision
wheretheaccuracyiscontrolledwithparameter(cid:15), and frequencyintheorderofν≈103 s−1,while,forahigh-
thenapplyingtheadaptiveSimpson’smethodseparately energyelectronwithE =1MeV,thefrequencyistwo
kin
onintervals[0,uc]and[uc,u](whenuc <u). ordersofmagnitudelessν≈10s−1.
However, evaluating the integrals this way at every
The simplest method which is suitable for the adap-
integration time step would be very inefficient. There-
tiveschemeistheMilsteinmethodwhichhasbothweak
fore,wecalculateandtabulatetheL andL valuesfor
0 1 and strong convergence of 1.0. When the evolution of
a wide range of u and Θ values before proceeding to
an N-dimensional stochasticvariable Y isgiven bythe
solve Eq. (46). From the cut-off limit, Eq. (48), we
followingItôform
can deduce that the integrals are sensitive to u and Θ
√
only near the curve u = Θ2+2Θ, but not elsewhere
dY =p(Y)dt+g(Y)dW, (49)
as Fig. 2 illustrates. This means that L and L do not
0 1
havetobetabulatedwithhighresolutionand,therefore,
the tabulated values can be used to efficiently compute theEuler-Maruyamamethodisacquiredbysimpledis-
thecoefficientsKandDateverytimestep. cretization: dY → ∆Y,dt → ∆t,anddW → ∆W. The
6
Milsteinmethod 5.1. Particleoperator
The Milstein method for the particle operator is ob-
(cid:88)N
Yi(tk+1)=Yi(tk)+pi(Y(tk))∆t+ gij(Y(tk))∆Wj tainedfrombyfirstsubstitutingthetensorG,Eq.(47),
intheextraterminEq.(50). Thenexttaskwouldbeto
j=1
discretize the double integral I, but this turns out to be
1(cid:88)N (cid:88)N ∂
+ g (Y(t )) g (Y(t ))I (t ), (50) fataltoourapproach.ThediagonalcomponentsofIcan
2 j=1 n=1 nj k ∂Yn ij k nj k be shown to be exactly Iii = (∆Wi)2 −∆t, but the non-
diagonaloneshavetobeapproximated[14,15]. There
hasaformsimilartotheEuler-Maruyama,withthedif- arewaystodothisbuttheycomewithaheftyprice:the
ference being the additional last term, which contains strongorderofconvergewouldbereducedto0.5.Since
thedoubleItôintegral wechosetheMilsteinmethodfortheexactreasonthat
its strong order is 1.0, this is obviously unacceptable.
I (t )=(cid:90) tk+1dW(t)dW (t). (51) Ontheotherhand,approximationsthatwouldmaintain
ij k i j
tk the strong order 1.0 are complicated and expensive to
compute, reducing the efficiency gains of the adaptive
In the adaptive scheme, the integration time step is
integration, so we seek a better way to implement the
varied so that local truncation error stays within user-
operator.
defined boundaries. For our purposes, we follow the
ThefactthatGisdiagonalin(u ,u )basisprovides
(cid:107) ⊥
scheme proposed in Ref. [12], which was later refined
uswithasuitablepath. Withonlydiagonalcomponents
in Ref. [13], in which two different error controls are
included,i.e. thenoisebeingcommutative,theMilstein
used. The first error control, defined in terms of the
methodreads
deterministiccomponentp,is
εdrift =i=m1,a...x,N(cid:40)2ε1abs,i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)pi∂∂Ypii(∆t)2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:41), (52) Yi(tk+1+)1=gYi((Ytk()t+))p∂i(Yg(tk()Y)∆(tt)+)(cid:104)g(i∆i(WY(t)k2)−)∆∆Wt(cid:105)i , (56)
2 ii k ∂Y ii k i
i
whereε isthetoleratederrorforthevariableY. The
abs,i i
seconderrorcontrolisdefinedintermsofthestochastic sothatthecollisionoperatorbecomes
componentg,
(cid:112)
εdiff,1 =i=m1,a...x,N6ε1abs,i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)gii(cid:32)∂∂gYiii(cid:33)2(∆Wi)3(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12). (53) u⊥u,(cid:107)j((ttkk++11))==1K(2uD(t⊥k)()u∆(ttk+))∆(cid:112)W2Dj,(cid:107)(ju=(tk1),)∆2,W3 (57)
+ D(cid:48)(u(t ))[(∆W )2−∆t], (58)
2 (cid:107) k 3
Alternatively,theseconderrorcontrolcanbedefinedin
termsofbothpandg (cid:88)2
ui(tk+1)=ui(tk)+u(cid:107)(tk+1)uˆi+ u⊥,j(tk+1)⊥ˆj (59)
(cid:40) 1 (cid:12)(cid:12) ∂p (cid:12)(cid:12)(cid:41) j=1
εdiff,2 =i=m1,a...x,N 6εabs,i (cid:12)(cid:12)(cid:12)gii∂YiiJi(cid:12)(cid:12)(cid:12) , (54) where (uˆ,⊥ˆ ,⊥ˆ ) form an orthogonal basis, and prime
1 2
whichisusefulwheng(cid:48)iszeroorexpensivetocompute. denotes partial derivative with respect to u. Equa-
HereJistheStratonovichintegral,J =(cid:82)tk+1(cid:82)t◦dWdt, tion (57) has formally reduced to the Euler-Maruyama
whichhasthevalue[13] i tk tk i method as the noise is additive, i.e., ∂D⊥/∂u⊥ = 0, so
this is a special case where even the Euler-Maruyama
1 (cid:32) ω (cid:33) method has strong convergence of 1.0. The error
J = ∆t ∆W + √i , (55)
i i terms, Eqs. (52) and (53), are now solely determined
2 3
byEq.(58):
√
where ω ∼ N(0,∆t). Approximating |ω| ≈ ∆t, we (cid:112)
get|J|=i(1/2)∆t(|∆W|+ √∆t/3). Onecainseethatthe εabs, =εtol(|K∆t|+ 2D(cid:107)∆t), (60)
i i
errorcontrolsgivenbyEqs.(53)and(54)areoftheor- |KK(cid:48)|(∆t)2
ε = , (61)
derO((∆t)3/2),whileEq.(52)hashigherorderO((∆t)2). drift 2ε
abs,
Thelatterthenseemsunnecessary. However,thisisnot (D(cid:48))2|∆W |3
thecaseinthelimitofweakdiffusion,whereg≈0,and εdiff = (cid:107) (cid:112) 3 . (62)
6ε D
thedeterministiccomponentpdominates. abs, (cid:107)
7
wherewehavechosenEq.(53)astheseconderrores- Equation (69) was reduced to the Euler-Maruyama
timate instead of Eq. (54). This is because K(cid:48) can be form as we have assumed uniform magnetic field, i.e.,
zero(recallFig.1),inwhichcasebotherrorestimates, (∂/∂X)(·) = 0. Inadditiontotheseequations,thereare
Eqs.(52)and(54),wouldyieldzerovalues. also boundary conditions. The particle pitch is limited
TheMilsteinmethodrequiresadditionalcomputation totheinterval[−1,1]whichcanbeenforcedwithare-
of D(cid:48),while K(cid:48) isneededfortheerrorestimates. Both flectingboundarycondition,
(cid:107)
can be evaluated analytically (D(cid:48) is included for com-
⊥
pletenesssake): ξ(tk+1)=sgn(ξ(tk+1))(2−|ξ(tk+1)|), (72)
∂K = Γabnb 1 (cid:34)2(cid:32)µ0 + maµ (cid:33) applied if |ξ(tk+1)| > 1. Similarly, u cannot have nega-
∂u m2c3 u3 γ m 1 tivevalues,soareflectingboundaryconditionshouldbe
a b
−u(cid:32)µγ(cid:48)0 + mmaµ(cid:48)1(cid:33)+u2µγ30(cid:35), (63) s(reetcaatlluF=ig0..1H),oswetetvinegr,tahsebroetflheνctainngdbKouunddivaerrygecoantduit=ion0
b
∂D Γ n Θ (cid:104) (cid:105) toasmallpositivev√aluebutstillbelowthe√thermalmo-
(cid:107) = ab b b uγ2µ(cid:48) −(1+2γ2)µ , (64) mentumvalueu ≈ 2Θ,e.g. atu = 0.05 2Θ,ensures
∂u m2c3 γu4 1 1
a that the time step in the adaptive scheme does not be-
∂D⊥ = Γabnb 1 (cid:104)(4γ2−1)Θ µ −uΘ γ2µ(cid:48) comeextremelysmall,andthatthemethodisstableifa
∂u m2c3 2γ3u4 b 1 b 1 fixedtimestepisused.
a
(cid:105)
−u2((2γ2−1)µ +Θ γ3µ )+u3γ2(µ(cid:48) +Θ γµ(cid:48)) , TheerrorestimatesfromEqs.(52)-(54)are
0 b 2 0 b 2
(65) (cid:16) (cid:112) (cid:17)
ε =ε |K |∆t+ 2D ∆t , (73)
abs,u tol u (cid:107)
wherethederivativesofthespecialfunctionsare ε =max(cid:12)(cid:12)(cid:12)QuQ(cid:48)u(cid:12)(cid:12)(cid:12), (cid:12)(cid:12)(cid:12)ξν2(cid:12)(cid:12)(cid:12) (∆t)2, (74)
µµ(cid:48)1(cid:48)0 == 22ΘΘγγuuLL10ΘΘ++γγe(cid:16)(eγ111//Θ+−ΘKK222ΘΘ2(cid:16)2(cid:16))Θ1u(cid:17)Θ1(cid:17)2u(cid:17)e2(e1(−1γ−)γ/Θ)/Θ, , ((6667)) εddriiffft =(cid:112)m1a−xξ2ν(cid:12)(cid:12)(cid:12)(cid:12)62D3εε/2(cid:107)aa2b(bs(cid:12)(cid:12)(cid:12)∆s∆,,uuWW(cid:112)u2ξD)ε3+(cid:107)a(cid:12)(cid:12)(cid:12)(cid:12)b,s√,ξ∆t/3)(cid:12)(cid:12)(cid:12)∆t
(cid:16) (cid:17) . (75)
2Θ2uL1+ 2Θ3γ+2Θ2+Θγ−u2 e(1−γ)/Θ 2εabs,ξ 
µ(cid:48)2 = Θ2γe1/ΘK (cid:16)1(cid:17) .
2 Θ Notethatwehavetwotolerances: ε forthemomen-
abs,u
(68) tum and ε for the pitch which we set ε = ε
abs,ξ abs,ξ tol
as the pitch values are bounded in an interval. Strictly
Onecannotethatµ(cid:48) =(u/γ)µ .
1 2 speaking,theerrortoleranceEq.(74),shouldbedefined
in terms of K , but using Q we only have to evaluate
5.2. Guidingcenteroperator u u
thederivative
WithMilsteindiscretization,theguidingcentercolli-
sionoperator,Eqs.(43)-(45),becomes ∂Qab,u =−Γabnb 1 ma(uµ(cid:48) −µ ). (76)
∂u m2c3 u3m 1 1
a b
(cid:88)3 (cid:113) (cid:16) (cid:17)
Xi(tk+1)= Xi(tk)+ 2DX,ij(Z(tk)) δij−bˆibˆj ∆WX,j, Qu is the dominant term in Ku outside the diffusion
j=1 dominated regime so this alteration is justifiable. For
(69) theξerrorterminEq.(75),wehadtochoosediffusion
(cid:112)
u(tk+1)=u(tk)+Ku(Z(tk))∆t+ 2D(cid:107)(Z(tk))∆Wu errorestimateEq.(54)insteadbecausetheestimatede-
rivedfromEq.(53)
1 (cid:104) (cid:105)
+ D(cid:48)(Z(t )) (∆W )2−∆t , (70)
ξ(tk+1)=ξ2(cid:113)(tk(cid:107))−ξk(tk)ν(Z(tuk))∆t (cid:113)2Dξ∂(cid:112)∂2ξDξ2 = (cid:112)1ξ−2 ξ2ν3/2, (77)
+ (1−ξ2(t ))ν(Z(t ))∆W
k k ξ
1 (cid:104) (cid:105) diverges when |ξ| = 1. No separate error limit was set
− ξ(t )ν(Z(t )) (∆W )2−∆t . (71)
k k ξ forthespatialcoordinateX.
2
8
5.3. OptimaltimestepandBrownianbridge The bias is avoided by introducing the conditioned
Nowthatthecriteriafortimesteprejectionareestab- probabilitydistributionknownastheBrownianbridge.
lished,thenexttaskistochooseanoptimaltimestepto First, letWt− andWt+ beadjacentrealizedWienerpro-
minimizethenumberofrejections.Agoodguessforthe cesseswitht− <t+. Whent− <t<t+,theprocessWt is
itnhietiearlrostreepsitsim∆attineist =co(cid:15)mt3o/pl2l/icνa,tbeustfitnhdeipnrgestheencoepotifm|∆alWti|mine nfoollloowngseardnioffremreanlltyndoirsmtriabludteisdtraibsuNtio(Wn wt−,h(etr−e tth−e)Im),ebaunt
step.Thesimplestscheme,knownastheBrowniantree, is
t−t
isbasedonhalvinganddoublingthecurrenttimestep. µ=W +(W −W ) − , (83)
However,thisschemeisfarfromoptimaland,therefore, t− t+ t− t+−t−
wechoosetoimplementthealgorithm,describedinde- andthevariance
tvaailluienoRfetfh[e12er]r,owrheesrteimthaetense.xItntitmheersetgeipmdeepoefnwdesaokndthife- Σ= (t−t−)(t+−t). (84)
t+−t−
fusion,εdrift >εdiff,wecouldtreatthecollisionoperator
asanODE,andchoosethenexttimestepas OnecanobservethatthevarianceofWt hasitsmaxima
at the center of the interval [t−,t+] while the expected
∆t(cid:48) =min(1.5,βε−1/2)∆t, (78) valuesfollowastraightlinefromW toW .Onlywhen
drift t− t+
where β < 1 is a safety factor for which we set value Wt+ fort+ >tdoesnotexist,Wt ∼N(Wt−,(t−t−)I).
β = 0.9. However, εdrift > εdiff could also be due to
extraordinary small |∆W|, not because drift dominates, 6. Verificationandbenchmark
sothenexttimestepischosenas
To summarize, we have now developed following
∆t(cid:48)
∆t =n , (79) collision operators: the fixed step Euler-Maruyama
next
3 method in the particle phase space, Eq. (46), the
where adaptive Milstein method in the particle phase space,
Eqs. (57) – (59), and the adaptive Milstein method in
n=max{l: |∆Wn|<∆Wopt, ∀l=1,...,3}. (80) theguidingcenterphasespace,Eqs.(69)–(71). From
nowon, werefertotheseasFEP,AMP,andAMG,re-
Here ∆W = βε−1/3|∆W| is the estimate for the “op-
opt diff spectively. These operators should yield equivalent re-
timal” value of Wiener process. In the diffusion domi-
sults which preferably are the same as those obtained
natedregime,εdiff > εdrift,weagaindeterminethenext analytically–atopicweinvestigatehere.
stepiterativelyfrom
Ourfirsttaskistoverifythatagiventestparticlepop-
∆t ulationrelaxestoMaxwell-Jüttnerdistribution,Eq.(8).
∆t =n , (81)
next 3 Inequilibrium,themagnitudeofthemomentumisdis-
tributedas
wherenisdeterminedbythecondition u∼u2e−√1+u2/Θ, (85)
n=max{l: |∆Wn|<∆Wopt, ∀l=1,...,lmax}, (82) and the pitch is distributed as ξ ∼ U(−1,1), where U
wherel√max =2ifthecurrentstepwasrejected,lmax =4 is the uniform distribution. According to Eq. (85), the
if ∆W/ ∆t < 2, i.e., the current Wiener process value test particle mean momentum should converge to 0.56
wasnotanoutlier,andotherwisel =6. andvarianceto3.62whenconsideringatestcasewhere
max
bothbackgroundandthetestparticlepopulationconsist
WhenevernewWienerprocessesaregenerated,these
of(relativistic)electronswithΘ = 1×10−1. Likewise
must always be stored – even when the time step they
mean pitch should converge to 0 and variance to 1/3.
areassociatedwithisrejectedortheyaregeneratedfor
From Fig. 3 we see that this is the case: All operators
thesolepurposeofdeterminingthenexttimestep.Only
whenintegrationhasreachedtimetcanprocessesW(t(cid:48)), convergetotheequilibriumvaluesatthesamerateand,
t(cid:48) < t,bediscarded. Thereasonforthisisthatthereal- thus,areverifiedinthisregard.
Theabovetestdidnotverifythespatialcollisionop-
izedvaluesconditionthedistributionoftheWienerpro-
erator, Eq. (69), which, in a uniform magnetic field,
cesses,sodiscardingthemcanleadtoabiasifthedis-
shouldcorrespondtotheclassicaldiffusiongivenbythe
cardmechanismissystematic.Asystematicmechanism
coefficient
canarise, e.g., fromerrorestimateEq.(53)asdiscard-
ingWienerprocesseswouldleadtoover-representation (cid:90) ∞ 1
ofsmallvaluesof∆W whengandg(cid:48)arelarge. Dc = du2ρ2(u)ν(u)f(u), (86)
0
9
a) b)
where ρ is the Larmor radius, ν is the pitch collision 0.9 0.3
frequency,and f(u)isthemomentumdistributionfunc-
0.25
tion. The spatial diffusion coefficient D can be esti- 0.8
B
wmhaeterdeDwijtharaetMesotnptaertCicalrelodimffeutshioondcaoseDffiBci=en(cid:80)tsn,j=a1nDdjn/nis, uE[]0.7 uSD[]0.01.52 FEP
AMP
either the number of test particles or, for a single par-
0.6 0.1 AMG
ticle, the number of time steps. It can be shown that Analytical
theratio Dj/DX obeysχ2 distributionwith1degreeof 0.5 0.05
0 0.01 0.02 0 0.01 0.02
freedom [16]. With the AMG operator, D is easy to
j Time(s) Time(s)
calculateasD =(∆X )2/2∆t ,where∆X isthechange c) d)
j j j
inguidingcenterposition,alongsomepredefineddirec- 0.2 0.6
tion,duringtime∆t.WiththeAMPorFEPoperator,the
0
collisionsonlyaffectparticlemomentumand,hence,do
-0.2 0.4
noopterlaetaodrtiosscpoautpialleddiwffuitshiotnh,euLnoleressnttzhefopracret.icFleorconlulimsieorn- ξE[]-0.4 ξD[]
S
ical evaluation of the Lorentz force, we use the energy -0.6 0.2
conservingscheme[17]. -0.8
Usingatestcasewherethetestparticlesarethermal-
-1 0
ized, and evaluating diffusion coefficient with different 0 0.01 0.02 0 0.01 0.02
operators, we obtain the test particle diffusion coeffi- Time(s) Time(s)
cient distributions shown in Fig. 4. We thus confirm
Figure3:Timeevolutionofthenumericalsolutionoftheparticledis-
thatalloperatorsyieldresultsmatchingwelltotheana- tributioncomputedwithdifferentmethods(FEP,AMP,andAMG)
lyticalresultcalculateddirectlyfromEq.(86). andthethermalequilibrium(Analytical). (a)Meanand(b)standard
Wehavenotyetshownthattherateatwhichtestpop- deviationofthemomentumdistribution. (c)Meanand(d)standard
deviation of the pitch distribution. Initially each test particle had
ulationrelaxestoequilibriumiscorrect,orthattheop- (cid:112)
erators are valid also in the regime of low diffusivity, u= (1+3Θ)2−1andξ=−1.
where the deterministic parts in the collision operators
dominate. Totacklebothissues,weperformaslowing-
(recallFig.1). Figure5showstheslowingdowndistri-
downsimulationusingthesamebackgroundplasmabut
thistimewithtestparticleshavinginitiallyu=5. Inthe butionestimatedthiswayoverlappingwiththedistribu-
tions obtained numerically with the different collision
slowing-down simulation, fast test particles are simu-
operators. Wecanthereforeconcludethattheoperators
lateduntiltheyreachacertainmomentumbelowwhich
yield a correct relaxation rate and are valid also in the
they can be considered part of the thermal population.
regimeoflowdiffusionality.
Here the common figure of merit is the slowing-down
time,i.e.,thetimeittakesforaparticletothermalize. Having verified the collision operators, it is time to
There are ways to analytically estimate the slowing- benchmarkthemtoseewhatcanbegainedbyusingthe
down time [18], but here we derive a simple estimate adaptivetimestep. Tothisend,weintroduceoneaddi-
from the theory of stochastic processes. The slowing- tionaloperator,FMG,whichistheguidingcenterAMG
downprocessisanalogoustotheproblemoffindingthe operator but with a fixed time step. The test case con-
so-called first passage time. Consider a stochastic 1D sists of simulating u = 5 electrons in Θ = 1 × 10−2
process defined by the Langevin equation (49) where plasma until they slow down below the energy corre-
coefficients pandgareconstant. Ifthisprocesshasan sponding to the background temperature. We measure
initialvalueY(t )=Y ,thenthetime,τ,ittakestoreach how the mean slowing-down time converges when de-
0 0
value Y +α for the first time, obeys an inverse Gaus- creasingthetimestepinthefixedschemesortheerror
0
siandistributionτ ∼ IG(α/p,α2/g2). Herethe1Dpro- tolerance in the adaptive ones. Here we noticed that
cessinquestionistheguidingcentermomentumequa- theoperatorsFEPandAMPconvergedto0.670swhile
tion (70), u is the particle initial momentum and u FMG and AMG converged to 0.673 s. This difference
0 min
is the momentum value below which a particle can be couldoriginatefromtheguidingcentertransformation,
considered thermalized. Setting u = 5 and u = 1, butitisinsignificantastheCoulomblogarithmisonly
0 min
wehaveα = −4. Furthermore, settingthebackground accuratetowithin1/lnΛ. Therateofconvergencefor
temperature to Θ = 1 × 10−2, the coefficients Q and alloperatorsisshowninFig.6,wheretheerrorisplot-
D areapproximatelyconstantontheinterval[u ,u ] ted with respect to the elapsed cpu time. In all cases,
(cid:107) min 0
10