New fast accurately conservative scheme for solving numerically the time-dependent isotropic Fokker–Planck equation

Accepted Manuscript New fast accurately conservative scheme for solving numerically the time-dependent isotropic Fokker-Planck equation Charlotte Boukandou-Mombo, Hassan Bakrim, Jonathan Claustre, Joëlle Margot, Jean-Pierre Matte, François Vidal

PII: DOI: Reference:

S0010-4655(17)30211-4 http://dx.doi.org/10.1016/j.cpc.2017.07.004 COMPHY 6263

To appear in:

Computer Physics Communications

Received date : 21 March 2016 Revised date : 8 June 2017 Accepted date : 6 July 2017 Please cite this article as: C. Boukandou-Mombo, H. Bakrim, J. Claustre, J. Margot, J. Matte, F. Vidal, New fast accurately conservative scheme for solving numerically the time-dependent isotropic Fokker-Planck equation, Computer Physics Communications (2017), http://dx.doi.org/10.1016/j.cpc.2017.07.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

New fast accurately conservative scheme for solving numerically the time-dependent isotropic Fokker-Planck equation Charlotte Boukandou-Momboa , Hassan Bakrima , Jonathan Claustrea , Jo¨elle Margotb, Jean-Pierre Mattea , Fran¸cois Vidala,c a Institut

national de la recherche scientifique-Centre ´ energie, mat´ eriaux et t´ el´ ecommunications, 1650 boul. Lionel Boulet, Varennes, Qu´ ebec, Canada J3X 1S2 b Groupe de physique des plasmas, Universit´ e de Montr´ eal, C.P. 6128, succ. centre-ville, Montr´ eal, Qu´ ebec, Canada H3C 3J7 c Author to whom correspondence should be addressed: [email protected]

Abstract We present a new numerical method for solving the time-dependent isotropic Fokker-Planck equation. We show analytically and numerically that the numerical scheme provides accurate particle and energy density conservation in practical conditions, an equilibrium solution close to the Maxwellian distribution, and the decrease of entropy with time. The slight nonconservation of particle and energy density is only due to the finite value of the upper bound of the energy grid. Additionally, the totally implicit scheme proves to provide positive solutions and to be unconditionally stable. The implicit forms of the scheme can be set as a nonlinear tridiagonal system of equations and solved iteratively. For a uniform grid in energy with N points, the number of operations required to compute the solution at a given time is only O(N ), in contrast to the totally explicit variant, which requires O(N 3 ) operations due to the restriction on the time step. The time-centered variant is more accurate than the totally implicit one, and uses an equivalent CPU time, but does not provide positive solutions for very large timesteps. The results of the method are analyzed for the classical problem of an initially Gaussian distribution as well as for an initially quasi-truncated Maxwellian distribution. Keywords: Fokker-Planck equation, finite differences, conservative scheme, entropy property, positivity, iterative solution, numerical stability

1. Introduction In physics the Fokker-Planck equation (FPE) describes the time evolution of the velocity or energy distribution function of particles subjected to many small changes of velocity or energy [1, 2], Preprint submitted to Computer Physics Communications

June 8, 2017

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

or of other properties [3] . In plasma physics, it is most often used to describe the cumulative effect of small deviation Coulomb collisions on the velocity distribution function of a given species of charged particles [4, 5, 6, 7]. When Coulomb collisions are strongly dominant, the energy distribution function tends towards a Maxwell-Boltzmann form fM (u) ∝ exp(−u/T ), where u is the energy and T =

2 3

< u > is the temperature, as this is the equilibrium solution of the FPE. The form

of the energy distribution in a plasma is the result of the competition between this collisional relaxation and other physical mechanisms. For example, in the Boltzmann modeling of microwave discharges, inelastic collisions tend to deplete the distribution function at energies above the first excitation threshold [8, 9, 10], and the heating (by microwaves or by DC fields) can also tend to give non-Maxwellian distribution functions [11, 12]. In this paper we discuss the FPE in the context of a gas of electrons without interactions with other species in the plasma (these interactions can be taken into account in separate steps in the simulation). Our treatment is limited to the isotropic part of the FPE. Even in simulations where electron transport is significant, the isotropic part of the FPE is an important part of the computation [13, 14, 15]. The isotropic part of the FPE reads: ∂G(u, t) ∂f (u, t) = −u−1/2 ∂t ∂u

(1)

where f (u, t) is the distribution function depending on energy u and time t, and   ∂f (u, t) ν(u) 3/2 I(u, t)f (u, t) + J(u, t) G(u, t) = −2 u n ∂u defines the collision term with I(u, t) =

Z

u

f (u′ , t)u′1/2 du′

(2)

(3)

0

2 J(u, t) = 3

Z

u

′

f (u , t)u

′3/2

′

du + u

3/2

0

Z

∞

′

f (u , t)du

u

′



(4)

The electron-electron Coulomb collision frequency is given by ν(u) =

4πne4 ln(Λ) 2

(4πε0 m) v 3

(5)

where n is the (conserved) electron density, m is the electron mass, ln(Λ) is the Coulomb logarithm, and v = (2u/m)1/2 is the velocity. 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In plasma simulations, electron-electron collisions have to be calculated frequently because electrons determine the plasma dynamics and their collisional thermalization, as described by the FPE, is among the fastest processes taking place in the plasma. This is true even for plasma discharges with a relatively low degree of ionization [10]. Solving the isotropic FPE in the context of complex time-dependent plasma simulations, such as cold plasmas which involve many atomic excited states [8, 9, 10], therefore requires an efficient numerical scheme with solutions respecting the key properties of the continuous FPE, namely positivity of the solutions, Maxwellian equilibrium, conservation of particle density n and of energy density ε, and the decrease of entropy S(t) with time (Boltzmann H theorem). The three latter properties read dn/dt = 0, dε/dt = 0, and dS/dt ≤ 0, respectively, where n= ε= S(t) =

Z

Z

∞

0 Z ∞

f (u, t)u1/2 du

(6)

f (u, t)u3/2 du

(7)

0

∞

f (u, t) ln (f (u, t)) u1/2 du

(8)

0

Explicit numerical schemes satisfying these key properties totally or partly have been developed for several years. [16, 17, 18] However the stability of explicit schemes imposes a severe constraint of the form ∆t < C(∆u)2 on the time step ∆t, where ∆u is the step in energy and C is a coefficient depending on f (u, t). [17] For N energy grid points, the latter inequality implies O(N 2 ) time steps. Since at least O(N ) operations are needed at every time step, a number of operations O(N 3 ) must therefore be performed to reach a given simulation time. Strategies have however been proposed to relax this constraint in the context of explicit methods. [19] Implicit schemes are generally more robust than explicit schemes. Several authors have developed finite difference implicit schemes for solving the FPE [20, 21, 22, 23, 24, 17, 25, 26, 27], but in these works either the above-mentioned key properties of the FPE are not all verified or efficient numerical schemes have not been demonstrated in practice. Based on ideas from Chang and Cooper [20], Bobylev and Chuyanov [21], and Langdon [22], first tested by Kho [23], Epperlein developed, and aplied to two-dimensional plasma simulations, numerical algorithms involving a nonuniform grid in velocity space [28, 24]. A similar approach has been applied in several studies of electron transport in laser heated plasmas, in one or two spatial dimensions, with or without self-consistent magnetic fields. [29, 30, 26, 31, 32] Besides conserving particles and energy (see however the critical

3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

discussion in [27]), this kind of approach also provides positivity and Maxwellian equilibirum due to the use of the Chang and Cooper method [20]. Somewhat later, Buet and Le Thanh [27, 33] investigated several new conservative implicit schemes and discussed how they respect the key properties of the FPE. In particular their schemes S1 and S3 , which have been put to work in [33], were shown to be conservative and to provide positivity and Maxwellian equilibrium. Most of the schemes mentioned above amount to solving a nonlinear tridiagonal system of equations, which can likely be solved in O(N ) operations using the fixed-point iterative method. A shortcoming of most iterative methods is that acceptable conservation is obtained only after a sufficient number of iterations, depending on the nature of the problem and the time step. Nevertheless, accurate energy conservation has been reported in test problems [23, 34] and in practical computations [26] using the fixed-point iterative method. As an alternative to iterative methods, conservative implicit schemes involving the single inversion of a dense matrix at each time step, have been developed [24, 35]. While such schemes require at least O(N 3 ) operations per time step when using standard resolution methods, the use of a fast Krylov solver, originally proposed in [36], allowed to reduce the number of operations to ≤ O(N 2 ) without preconditionning [35]. Epperlein has reported that when the timestep ∆t is considerably longer than the thermal collision time τ , it can be more advantageous to perform such a full matrix inversion than to iterate [24]. We have found that, with our scheme, even when τ ≪ ∆t, there is a considerable increase of the number of iterations required, but the scaling is still more favorable than O(N 2 ).

Besides finite difference schemes, other approaches based for instance on eigenfunction expansion, path integrals, Green’s functions, the method of moments, the Monte Carlo method, or wavelets [37] have been used for solving the FPE (see [38] and references therein). In this article, we present a simple finite difference method, using a uniform grid in energy, whose implicit forms can be set as a nonlinear tridiagonal system of equations and solved iteratively using the fixed-point iterative method. We show analytically and numerically that the above-mentioned key properties of the FPE (conservation of particles and energy, positivity, Maxwellian equilibirum and decrease of entropy) are all verified provided the upper bound of the energy grid N ∆u is large enough. This is made possible simply by discretizing appropriately the integrals I(u, t) and J(u, t). As calculating the matrix elements and solving the tridiagonal system of equations both require O(N ) operations, the total number of operations for each time step as well as for the whole simulation are also O(N ). This was actually checked by measuring the CPU time as a

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

function of N . The test problems anayzed in section 4 show that the totally implicit scheme (TIS) is unconditionally stable, and converges toward the equilibrium solution for arbitrarily large time steps. Moreover, the accuracy criterion for iterations, as given by Eq. (18), is improved by several orders of magnitude at each iteration. Although many of the most interesting features of our scheme are likely shared by some of those discussed above ([23, 24, 26, 33]), ours is different, likely simpler, its properties are shown to satisfy all the key properties of the FPE, and we make a more extensive numerical validation of it. In addition, we investigate the stability and accuracy of the scheme for different parameters θ, which characterizes the relative amount of explicit and implicit features of the scheme. This article is organized as follows. In section 2 we present the numerical scheme. In section 3 we analyze the properties of the discretized FPE. In section 4, we apply the numerical scheme to two tests problems and we evaluate the performance of the method in terms of CPU time as well as conservation and other properties. Section 5 is a brief summary.

2. The numerical scheme We discretize energy as ui = (i − 1/2)∆u, with i = 1, 2, ..., N and time as tk = k∆t, where k = 0, 1, 2, ..., M . The finite-difference equation corresponding to Eq. (1) is i 1 h fik+1 − fik −1/2 k+θ Gk+θ − G = −ui i+1/2 i−1/2 ∆u ∆t

(9)

where fik ≡ f (ui , tk ) and Gk+θ i+1/2 ≡ G(ui +∆u/2, tk +θ∆t) with 0 ≤ θ ≤ 1. The distribution function

fik is assumed to be uniform in each energy cell ui−1/2 ≤ u < ui+1/2 , where ui±1/2 = ui ± ∆u/2.

This avoids having a point at u = 0 which would give a sigularity in Eq. (9). We define k+1 k Gk+θ i+1/2 ≡ (1 − θ)Gi+1/2 + θGi+1/2

where Gki+1/2

= −2g

"

k Ii+1/2

k fi+1 + fik 2

!

+

k Ji+1/2

(10)

k fi+1 − fik ∆u

!#

(11)

with g = ν(u)n−1 u3/2 , a constant. The discretized FPE (9) becomes k+1 k+1 Ak+1 fi−1 + Bik+1 fik+1 + Cik+1 fi+1 = Dik i

5

(12)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where

  ∆t −1/2 1 k+1 1 k+1 − Ii−1/2 ui + Ji−1/2 ∆u 2 ∆u     1  k+1 ∆t −1/2 1 k+1 tk+1 k+1 − Ii+1/2 − Ii−1/2 Ji+1/2 + Ji−1/2 ui Bik+1 = 1 − 2gθ ∆u 2 ∆u   ∆t −1/2 1 k+1 1 k+1 Cik+1 = −2gθ u I + J ∆u i 2 i+1/2 ∆u i+1/2  
1−θ k k 1−θ 1−θ k k 1 − Bik fik − Ai fi−1 + 1 + Ci fi+1 Dik = − θ θ θ Ak+1 = −2gθ i

(13)

k+1 k+1 Note that Ak+1 = Ak1 = 0 (since, from Eqs. (3) and (4), I1/2 = 0 and J1/2 = 0), so that 1 k+1 f0k+1 does not need to be defined. At the opposite boundary condition, we normally set fN +1 = 0, k+1 k+1 k+1 but it will be seen below that setting fN fN with the appropriate choice of αk+1 such +1 = α

that 0 < αk+1 < 1 enforces the conservation of particles or energy. One thus needs to make the k+1 k+1 k+1 ← BN + αk+1 CN . substitution BN

The way of discretizing the integrals I(u, t) and J(u, t) is critical to meet the key properties of the continuous FPE, especially the conservation of energy density (see section 3). The main idea of the numerical method described here is to define these integrals as  2 X k  3/2 3/2 fj uj+1/2 − uj−1/2 3 j=1 i

k = Ii+1/2

k = Ji+1/2



i X



1

(14)

N X



2 3/2 3/2 3/2 fjk  ∆u + uj−1/2 + ui+1/2 fk u 3 j=1 j 2 j+1/2 j=i+1

(15)

k k with I1/2 = 0 and J1/2 = 0 and ui+1/2 = i∆u. Eqs. (14) and (15) have been obtained through the

following approximations, which become exact when ∆u → 0,  2 1  3/2 1/2 3/2 ui ≃ ui+1/2 − ui−1/2 3 ∆u  1  3/2 3/2 3/2 ui+1/2 + ui−1/2 ui ≃ 2

1/2

in the trivial discretized forms of I(u, t) and J(u, t) containing ui

(16) 3/2

and ui .

k , fik and In the case θ = 0, Eqs. (12) and (13) indicate that fik+1 is obtained directly from fi−1 k fi+1 . This corresponds to the totally explicit scheme (TES). However, for θ > 0 Eq. (12) appears

to be a nonlinear system of equations which can be set in a tridiagonal matrix form and solved iteratively. Defining the upper indices in Eq. (17) as the iteration number, Eq. (12) is, for a given k, formally equivalent to m m Am−1 fi−1 + Bim−1 fim + Cim−1 fi+1 = Di0 i

6

(17)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

with m taking successively the values 1, 2, 3, ..., the upper index 0 corresponding to the time tk . Equation (17) is solved for each value of m. This is the fixed-point method. The iteration process is stopped when the relative change of the distribution function satisfies the criterion m fi f m−1 − 1 ≤ ǫ i

(18)

for all i, where ǫ is the stopping parameter.

The stability analysis of the numerical scheme is difficult due to the nonlinearity of Eq. (12). In section 4 we check numerically that the TES (θ = 0) requires O(N 3 ) operations to achieve a given simulation time using the largest time step ∆t for which the scheme remains stable. This indicates that the TES is subject to a stability condition of the form ∆t < C(∆u)2 , where C is some prefactor, which is typical of diffusion-type equations. Furthermore, we check that the TIS (θ = 1) is unconditionally stable and requires only O(N ) operations for reaching a given simulation time. In addition to stability, it is however essential that the numerical scheme possesses the same key properties as the continuous FPE. This is the subject of the next section. 3. Properties of the discretized Fokker-Planck equation In this section we discuss the preservation of the key properties of the continuous FPE (1) by the discretized FPE (12), namely the conservation of particle and energy density, the Maxwellian equilibrium, the positivity of the solutions and the decrease of entropy. 3.1. Particle density conservation From Eq. (9), one sees immediately that nk+1 − nk = −∆tGk+θ N +1/2 where nk =

PN

j=1

(19)

1/2

uj fjk ∆u. Therefore, the particle density is well conserved when uN is large

enough so that the right-hand side of Eq. (19) is much smaller than n0 /M . One notes that k k k k exact particle density conservation can be enforced by setting fN +1 = αn fN , with αn such that

Gk+θ N +1/2 = 0. From Eqs. (11), (14) and (15), one finds αkn so that 0 < αkn < 1, as expected.

=

PN

3/2

j=1 PN j=1

7

fjk uj−1/2 3/2

fjk uj+1/2

(20)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.2. Energy density conservation From Eq. (9), one sees that 

where εk =

PN

j=1

3/2

εk+1 − εk = ∆t 

N X j=1



k+θ  Gk+θ j+1/2 ∆u − GN +1/2 uN +1

(21)

uj fjk ∆u. Using Eqs. (14) and (15) as well as the following general identity for

arbitrary ai and bi (i = 1, 2, ..., N ) N X i=2

one obtains

ai

i−1 X j=1

bj =

N −1 X i=1

bi

N X

aj

(22)

j=i+1

N X

N 4 k X k 3/2 Gkj+1/2 ∆u = −g fN ∆u f u 3 +1 j=1 i j+1/2 j=1

(23)

Thus, as for particle density conservation, the energy density is well conserved provided uN is large k k k enough so that the right-hand side of Eq. (21) is much smaller than ε0 /M . Setting fN +1 = αε fN ,

one finds that εk+1 − εk = 0 provided that αkε = (uN +1 /uN )αkn , where αkn is given by Eq. (20). One

notes that 0 < αkn < αkε < 1. The fact that αkε 6= αkn prevents the simultaneous exact conservation

of particles and energy density. It is worth noting that in all practical kinetic simulations, uN /T must be taken large enough to capture the physical effects taking place in the high-energy range (such as multiple superelastic collisions between electrons and highly excited atomic states, for instance). The slight nonconservation of energy and particles is not due to the finite difference scheme, but rather to the finite value of the maximum energy uN . Even with a continuum representation of the distribution up to a maximum energy, there would be a similar nonconservation, because diffusion in energy space would diffuse particles beyond this maximum energy. Therefore enforcing particles or energy consevation for inadequate values of uN /T can provide physically wrong results. For values of uN /T typically used in plasma kinetic simulations, particle and energy density are well conserved using k the boundary condition fN +1 = 0, as will be illustrated in section 4.

3.3. Maxwellian equilibrium One can show that G(fM ) = 0, where fM is the Maxwell distribution u 2 n − T √ e fM (u) = π T 3/2 8

(24)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

with T = (2/3)ε/n, the temperature. Hence fM (u) is an equilibrium solution of the continuous FPE (1). Substituting fM (ui ) in Gi+1/2 given by Eq. (11), one obtains for i ≪ N : Gi+1/2 ∼ nν(T )

∆u  ui 3/2 −2ui /T e T T

(25)

For i ∼ N the corresponding expression depends on the boudary condition used for fN +1 . From

Eq. (25) we conclude that the equilibrium solution fi∗ of the discretized FPE coincides with fM (ui )

when ∆u/T → 0. Numerical results show that fi∗ and fM (ui ) actually look always close even for values of ∆u/T not much smaller than 1, particularly if we use fN +1 = αn fN or fN +1 = αε fN instead of fN +1 = 0. 3.4. Positivity Chang and Cooper [20] have shown that the solution fik+1 of the tridiagonal system of equations ≤ 0, Bik+1 ≥ 0, Cik+1 ≤ 0, and Dik ≥ 0. One can show that the (12) is positive provided that Ak+1 i first three inequalities are always verified here using the discretizations (14) and (15). However, the condition Dik ≥ 0 is guaranteed only when θ = 1 (TIS), in which case Dik = fik . In that case, an

initial fi0 > 0 will thus provide fik+1 > 0 for k ≥ 0. When θ < 1 the sign of Dik is undetermined and

therefore positivity is not assured. As discussed in section 4, in the case θ = 1/2, which corresponds to the time-centred implicit scheme (CIS) (accurate to O(∆t2 ) in contrast to the TIS which is only accurate to O(∆t)) we found that fik+1 shows oscillations when ∆t is greater than some critical value and then negative values for larger ∆t. This critical ∆t is however much larger than for the TES (θ = 0) for the same parameters, and, as will be seen below, it does not depend on ∆u, contrary to the TES. 3.5. Decrease of entropy The solution of the continuous FPE obeys the Boltzmann H theorem dS/dt ≤ 0, where entropy S(t) is given by Eq. (8), with the equality reached at equilibrium, that is when f (u, t) is the Maxwell distribution (24). The entropy property can be set in the more convenient form Z ∞ dS(t) ∂f (u, t) = ln (f (u, t)) u1/2 du ≤ 0 dt ∂t 0

(26)

Using Eq. (9) the discretized form of Eq. (26) reads (we omit the time index k + θ) N −1 X i=1

Gi+1/2 ln



fi+1 fi

 9

− GN +1/2 ln(fN ) ≤ 0

(27)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

We prove that this inequality holds when replacing ln(fi+1 /fi ) by the approximation fi+1 /fi − 1, which is accurate when fi+1 /fi ≈ 1, that is

  fi+1 + GN +1/2 (ln(fN ) − 1) ≥ 0 Gi+1/2 1 − fi i=1

N X

(28)

Indeed, using the general equality (22) and Eqs. (11), (14) and (15), one shows that the first term of Eq. (28) is a positive constant multiplied by N −1 X

N X

i=2 j=i+1

(i − 1)3/2

N X fj2 (fj fi−1 − fi fj−1 )2 (j − 1)3/2 + fN fj−1 fi−1 fj−1 j=2

which is clearly positive when fi > 0. Consequently the inequality (28) is verified when uN is large enough so that GN +1/2 (ln(fN ) − 1) is negligible compared to the first term. Since Gi+1/2 is

in some sense small at the equilibrium solution fi∗ , one sees that the rate of decrease of entropy (proportional to the left-hand side of Eq. (27)) vanishes as fi approaches fi∗ . We note that the

assumption fi+1 /fi ≈ 1 used above amounts to assuming that the energy grid is fine enough to describe accurately the distribution function. However, a more general demonstration of Eq. (27) could possibly be found.

4. Test problems In the following we use the normalizations u/u0 ← u, where u0 is an arbitrary reference en3/2

ergy, t/τ ← t, where τ = ν −1 (u0 ) is the collision time, and u0 f (u/u0)/n ← f (u). Using these normalizations, the discretized FPE still looks as Eqs. (9-15) except that g = 1 in Eq. (11). In the following test problems, we set the system of nonlinear equations (12) obtained for θ > 0 in a tridiagonal matrix form and solve it iteratively using the fixed-point method, described by Eq. (17). The properties of the discretized FPE discussed in the previous section hold for an exact solution of Eq. (12), so that a very small stopping parameter ǫ must be used in the stopping criterion (18). We first apply the numerical method to the relaxation of an initially Gaussian energy distribution "  2 # u − uG f (u, t = 0) = CG exp − (29) δG for u ≥ 0, analogous to the classical test problem treated initially by MacDonald, Rosenbluth and Chuck [16] and subsequently by Kho [23] and by Epperlein [34]. In our example, we use 10

Initial distribution 0.05 0.2

1

0.6 1 5

Maxwellian distribution

M

f(u) / f (0)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

Normalized energy, u

Figure 1: Energy distribution at different normalised times for an initially Gaussian function. The calculation parameters are given in the text.

uG = 1 and δG = 0.25 as the energy spread. The energy grid is chosen such that uN = 25 and N = 1000, implying that ∆u ≈ 0.025. The discrete integrals of Eq. (29) yield CG ≈ 2.2659 for n = 1 and ε ≈ 1.0159, or T ≈ 0.6773. In the following we used the boundary condition fN +1 = 0, unless otherwise specified, since for the ratio uN /T ≈ 37 used here, there is no advantage in using the boundary condition fN +1 = αfN discussed in sections 3.1 and 3.2. This ratio is typical of electron kinetic simulations. The value ∆u ≈ 0.025 was selected to provide a good sampling of the distribution functions and smooth graphs of fi vs ui . However, we have checked that the following numerical results remain nearly unchanged when using a tenfold larger ∆u, with, for instance, uN = 25 and N = 100 or uN = 250 and N = 1000. Figure 1 shows the distribution function for different times t on a linear scale while figure 2 shows the same data on a logarithmic scale to emphasize the high energy behaviour. Here we used the TIS with ∆t = 1/300 and ǫ = 10−10 . On the plots, the distributions are normalized with respect to fM (u = 0). One sees that at t = 5 the distribution function can hardly be distinguished from the Maxwellian distribution fM (u) with T ≈ 0.6773 in the energy range u . 5. The distribution function would need, however, more time to reach equilibrium at higher energies due to the smaller values of the Coulomb collision frequency ν(u), Eq. (5). We checked that the calculated distribution function gets closer and closer to the Maxwellian distribution at high energy as time t increases. However, with the fN +1 = 0 boundary condition, this is not the case very close to uN , whereas

11

it does hold all the way up to uN if a boundary condition fN +1 = αn fN or fN +1 = αε fN is used. Specifically, we found that the numerical solution matched the exact equilibrium solution (Maxwellian) extremely well for all energies up to uN = 25 at t ≈ 100, for all the values of the timestep which were tried, from ∆t =

1 30

up to ∆t = 100, although in the latter case, 2000 iterations

were required to satisfy the very strict tolerance which we demanded, ǫ = 10−10 .

1

Initial distribution 0.05 0.1

0.2 0.6

0.01

1 5

Maxwellian distribution

1E-3

M

f(u) / f (0)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1E-4

1E-5

1E-6

1E-7 1

2

3

4

5

6

7

8

9

10

11

12

Normalized energy, u

Figure 2: Energy distribution at different normalised times for an initially Gaussian function. As in figure 1 but on a logarithmic scale and for u ≥ 1.

The relative conservation error on the energy density as a function of time is shown in figure 3 for different time steps. The relative error grows steadily, reaching O(10−11 ) at t = 5 depending somewhat on the time step. The relative error on the particle density (not shown) is rather of O(10−14 ). We checked that these relative conservation errors increase when uN is reduced, in accordance with the discussions of sections 3.1 and 3.2. For example, for the boundary condition fN +1 = 0 and the parameters ∆u = 0.025, ∆t = 1/5, t = 5 and uN = 12.5, the relative errors on particle and energy density are O(10−10 ) and O(10−9 ), respectively, which are still acceptable. However, for uN = 6.25, all the other parameters being the same, the corresponding values are both O(10−3 ). In the latter case, when using the boundary condition fN +1 = αn fN (see section 3.1), the relative errors become O(10−14 ) and O(10−4 ), respectively. When using the boundary condition fN +1 = αε fN (see section 3.2) the relative errors become O(10−4 ) and O(10−10 ), respectively. In the latter case, the accuracy of the energy density conservation can be improved by decreasing the stopping parameter ǫ. For example, for ǫ = 10−14 the relative errors become O(10−4 ) and O(10−14 ), 12

respectively.

)

3

1/150 1/20 1/5

2

1

1

R

elative error on energy density (10

-11

1/300

0 0

1

2

3

4

5

Normalized time, t

Figure 3: Relative conservation error on energy density for different time steps ∆t using ǫ = 10−10 .

10

9

8

Number of iterations

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

7

6

5

4

3

2

1

0 0

1

2

3

4

5

Normalized time, t

Figure 4: Number of iterations for different stopping parameters ǫ using ∆t = 1/300. Lines are drawn to guide the eye.

The number of iterations required to satisfy the convergence criterion (18) as a function of time is shown in figure 4 for different values of the stopping parameter ǫ. As expected, the number of iterations decreases as ǫ increases, and as the distribution function approaches equilibrium. Not more than 9 iterations are required per time step for the smallest stopping parameter considered. In addition, there is a difference of only two iterations between ǫ = 10−6 and 10−10 . Since the 13

relative error on the energy density at t = 5 appears to grow almost proportionaly to ǫ, increasing ǫ excessively is not an option. The TIS proves to allow arbitrarily large time steps while preserving all of the key properties of the FPE discussed in section 3. For instance, a solution close to the expected Maxwellian distribution is obtained using a single time step with ∆t ≫ 1. The number of iterations corresponding

to the stopping parameter ǫ = 10−10 is represented as a function of time in figure 5 for different time steps ∆t. As expected, the number of iterations increases with the time step. One notes that, between ∆t =1/300 and ∆t = 1, the number of iterations increases only by a factor of 7-8, which indicates that increasing the time step is clearly advantageous in terms of calculation time.

70

1/300 60

1/20 1/2

Number of iterations

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

50

1

40

30

20

10

0 0

1

2

3

4

5

Normalized time, t

Figure 5: Number of iterations for different time steps ∆t using ǫ = 10−10 . Lines are drawn to guide the eye.

The latter advantage is however offset by a loss in precision, which is all the more important that the TIS is accurate only to O(∆t). Figure 6 shows the relative differences between the distribution functions calculated with different time steps and a reference distribution corresponding to ∆t = 1/300 at the same time t = 0.2 (therefore still far from equilibrium). As required by particle and energy density conservation, the relative difference oscillates around the reference distribution. A maximum relative difference of about 10% is obtained in the case of a single time step of ∆t = 0.2. As a result, a compromise must be found between speed and accuracy. The previous calculations, done with the TIS (θ =1), were redone by using the CIS (θ = 1/2), which has an accuracy of O(∆t2 ). When ∆t is small, no significant differences are observed between the solutions of the two schemes. The differences increase, however, when ∆t increases, as can be 14

CIS 1/150

elative difference of solutions

0.10

R

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

TIS 1/150 CIS 1/20 TIS 1/20

0.05

TIS 1/5

0.00

-0.05

-0.10

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Normalized energy, u

Figure 6: Relative differences between the solutions obtained at t = 0.2 using different ∆t for the TIS and the CIS with respect to the solution obtained with ∆t = 1/300. The curve for the CIS with ∆t = 1/5 is not shown since the scheme is unstable.

seen in figure 6. For a given time step the solutions obtained from the CIS are closer to the reference solution as compared to the TIS. However, for ∆t & 1/6, the CIS becomes unstable and starts producing large oscillations around the expected solution. The solution eventually becomes negative as ∆t further increases. For this reason the case ∆t = 1/5 for the CIS is not shown in figure 6. Besides producing more accurate solutions within its stability regime, the CIS proves to require fewer iterations to satisfy the convergence criterion as shown in figure 7. However, we found little difference in the overall computation times between the two schemes. To quantify the performances of the numerical method in terms of calculation time, we have measured the computation time (CPU time) as a function of the number of points N in energy, using uN = 25 in all cases, for the TES, TIS and CIS. As discussed above, the stability of the TES requires that ∆t < C(∆u)2 , with C depending on fik [17] However, because of the nonlinearity of the FPE, the exact condition is difficult to determine analytically. The search for the maximum allowed ∆t was therefore conducted empirically, i.e. ∆t was gradually increased until the appearance of large oscillations around the expected solution. For example, for the standard parameters used in the Gaussian test problem, the maximum time step for the TES is estimated to be ∆t ≈ 1/4000. Figure 8 represents the CPU time as a function of N for the TIS and the TES at t = 5. The CPU times of the CIS are practically the same as those of the TIS and are therefore not represented. All 15

20

TIS 1/300 CIS 1/300 TIS 1/10

Number of iterations

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

15

CIS 1/10

10

5

0 0

1

2

3

4

5

Normalized time, t

Figure 7: Number of iterations as a function of time for the TIS and the CIS for ∆t = 1/300 and 1/10 using ǫ = 10−10 . Lines are drawn to guide the eye.

the curves follow approximately straight lines on a log-log scale. As can be seen directly, the slope of the curve corresponding to the TES is about 3 while the slopes for the TIS (and for the CIS) are parallel and close to 1, indicating that the number of operations is O(N 3 ) and O(N ), respectively. The latter result is essentially due to the facts that the limits on ∆t do not depend on ∆u and that only O(N ) operations are necessary to evaluate the matrix elements and to solve the tridiagonal system of equations. As a second test problem we consider the case of an initially quasi-truncated Maxwellian, defined as

   u   u ≤ uQ exp −     T1  f (u, 0) = CQ u − uQ uQ   exp − exp − u > uQ T1 T2

(30)

with uQ = 1, T1 = 0.25 and T2 = 0.005. This example is relevant to simulations of cold plasmas where the distribution function is often almost cut off near the energy corresponding to the first excited state of the gas atoms [8, 9, 10]. Here, we are interested in the relaxation of this initial distribution towards equilibrium using the TIS. As in the Gaussian test problem, we use N = 1000, uN = 25 and ǫ = 10−10 . CQ ≈ 9.4380 provides n = 1 while the energy density is ε ≈ 0.3312, or T ≈ 0.2208. Figure 9 represents the distribution function for different times obtained with ∆t = 1/300. One observes the gradual convergence towards the Maxwellian distribution corresponding to the 16

TIS 1/300

CPU Time (seconds)

1000

TIS 1/20 TIS 1 TES (variable)

100

10

1

0.1

0.01

1000

2000

4000

6000

8000 10000

Number of points in energy, N

Figure 8: Measured CPU time at t = 5 for the TES (maximum time step ∆t preserving stability) and the TIS for different ∆t using ǫ = 10−10 . The CPU times of the CIS (not shown) are practically the same as those of the TIS.

100

Initial distribution

1

5

0.0 0.01

0.2 0.6

1E-4

M

f(u) / f (u=0)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1

5

1E-6

M

axwellian distribution

1E-8 1E-10 1E-12 1E-14 1E-16 1E-18 1E-20 0

1

2

3

4

5

6

7

8

9

10

Normalized energy, u

Figure 9: Energy distribution at different normalised times for an initially quasi-truncated Maxwellian distribution. The calculation parameters are given in the text.

temperature T ≈ 0.2208. The sharp variation of the initial distribution does not affect the stability of the calculation. Indeed, whatever the time step, the TIS was perfectly stable, providing positive solutions and always leading ultimately to a solution close to the expected Maxwellian distribution after a few collision times (longer time is need at high energies, as mentioned previously). In contrast, the TES shows instabilities for ∆t & 1/6000. All the observations made for the Gaussian test problem still hold here. In particular the relative error on the energy density is O(10−10 ). 17

-100 -10

Gaussian Quasi-truncated

-1

M

axwellian

-0.1

dS/dt (arb. units)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

-0.01 -1E-3 -1E-4

5

-1E-

-1E-6

7

-1E-

-1E-8

9

-1E-

-1E-10 -1E-11 0

1

2

3

4

5

Normalized time, t

Figure 10: Time evolution of the rate of entropy decrease for the initially Gaussian distribution and the quasitruncated Maxwellian distribution using the TIS. The calculation parameters are the same as for Figs. 1,2 and 9.

Finally we show in figure 10 the time evolution of the rate of change of entropy as given by the left-hand side of Eq. (27) for the two test problems considered. As expected, there is a steady decrease of entropy, approaching a constant value as time increases. We note that in the case of 0 the initially quasi-truncated Maxwellian the ratio fi+1 /fi0 is O(10−4 ) in the fast transition region,

near u = 1, so that the distribution function is not well described in this area. Nevertheless, the overall time evolution of entropy is not affected by this initial lack of accuracy.

5. Conclusion In summary, we have described a new fast and efficient numerical method for solving the isotropic FPE on a uniform energy grid. The fully implicit version of the method retains all the key properties of the continuous FPE, namely the accurate conservation of particle and energy density, positivity of solutions, decrease of the entropy in time and a Mawellian equilibrium solution. This was made possible by discretizing appropriately the integrals I(u, t) and J(u, t) defining the collision term in the FPE. These key properties are satisfied provided some conditions on the energy grid are respected, in particular the upper bound of the energy grid must be large enough. A parameter θ allows adjusting centering in time. The implicit forms of the method amount to solving a nonlinear tridiagonal system of equations, which can be done iteratively with relatively few iterations. As

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

expected, the totally implicit scheme (θ = 1) proves to be unconditionally stable. Although the time-centered implicit scheme (θ = 1/2) proves to be more accurate (O(∆t2 )) than the totally implicit scheme (O(∆t)), it produces negative solutions for large time steps. However the threshold value for instability is much larger than for the totally explicit scheme (θ = 0) and it is independent of the energy grid resolution. The numerical method was tested on two problems : the initially Gaussian distribution and the initially quasi-truncated Maxwellian. The CPU time for the totally explicit scheme grows as O(N 3 ) within the stability limit, and as O(N ) for the totally implicit and time-centered schemes. The CPU time measured for the time-centered implicit scheme appears to be practically the same as that of the totally implicit scheme, for the same time step and energy grid. The gain in calculation time provided by the implicit schemes as compared to the explicit scheme is however offset by a loss of precision in the solutions when excessively large time steps are used. Therefore, a compromise must be adopted between speed and accuracy. Such a fast practically conservative algorithm constitutes a significant asset in time-dependent simulations of cold plasmas involving a great number of atomic excited states and reactions. The next useful step for this purpose would be to extend the scheme to a non-uniform energy grid in order to decrease the density of grid points at high energy, and hopefully speed-up the caclulations in one- or two-dimensional geometries.

ACKNOWLEDGEMENTS

This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Strategic Project Grants Program (Grant number 447624-13).

19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

References [1] A. D. Fokker, Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld, Ann. Phys. 43 (1914) 810–820. ¨ [2] M. Planck, Uber einen Satz der statistischen Dynamik and seine Erweiterung in der Quantentheorie, Sitzungsberichte der Preussischen Akademie der Wissenschaften 24 (1917) 324–341. [3] J. Y. Rempel, M. G. Bawendi, K. F. Jensen, Insights into the kinetics of semiconductor nanocrystal nucleation and growth, J. Am. Chem. Soc. 131 (12) (2009) 4479–4489. [4] I. P. Shkarofsky, T. W. Johnston, M. P. Bachynski, The Particle Kinetics of Plasmas, AddisonWesley, Reading, Massachusetts, 1966. [5] J. Oxenius, Kinetic Theory of Particles and Photons, Springer-Verlag, New York, 1986. [6] A. G. Peeters, D. Strintzi, The Fokker-Planck equation, and its application in plasma physics, Ann. Phys. 17 (2008) 142–157. [7] A. G. R. Thomas, M. Tzoufras, A. P. L. Robinson, R. J. Kingham, C. P. Ridgers, M. Sherlock, A. R. Bell, A review of Vlasov-Fokker-Planck numerical modeling of inertial confinement fusion plasma, J. Comput. Phys. 231 (2012) 1051–1079. [8] L. L. Alves, C. M. Ferreira, Electron kinetics in weakly ionized helium under DC and HF applied electric fields, J. Phys. D: Appl. Phys. 24 (4) (1991) 581. [9] L. L. Alves, G. Gousset, C. M. Ferreira, A collisional-radiative model for microwave discharges in helium at low and intermediate pressures, J. Phys. D: Appl. Phys. 25 (12) (1992) 1713. [10] M. Santos, C. Noel, T. Belmonte, L. L. Alves, Microwave capillary plasmas in helium at atmospheric pressure, J. Phys. D: Appl. Phys. 47 (26) (2014) 265201. [11] M. J. Druyvesteyn, F. M. Penning, The mechanism of electrical discharges in gases of low pressure, Rev. Mod. Phys. 12 (1940) 87–174. [12] K. Behringer, U. Fantz, Spectroscopic diagnostics of glow discharge plasmas with nonMaxwellian electron-energy distributions, J. Phys. D: Appl. Phys. 27 (10) (1994) 2128–2135.

20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[13] A. R. Bell, R. G. Evans, D. J. Nicholas, Elecron energy transport in steep temperature gradients in laser-produced plasmas, Phys. Rev. Lett. 46 (1981) 243–246. [14] J. P. Matte, J. Virmont, Electron heat transport down steep temperature gradients, Phys. Rev. Lett. 49 (1982) 1936–1939. [15] F. Alouani-Bibi, M. Shoucri, J. P. Matte, Different Fokker-Planck approaches to simulate electron transport in plasmas, Comput. Phys. Commun. 164 (1-3) (2004) 60–66. [16] W. M. MacDonald, M. N. Rosenbluth, W. Chuck, Relaxation of a system of particles with Coulomb interactions, Phys. Rev. 107 (2) (1957) 350–353. [17] C. Buet, S. Cordier, Numerical analysis of the isotropic Fokker-Planck-Landau equation, J. Comput. Phys. 179 (1) (2002) 43 – 67. [18] N. Crouseilles, F. Filbet, Numerical approximation of collisional plasmas by high order methods, J. Comput. Phys. 201 (2) (2004) 546 – 572. [19] F. Filbet, L. Pareschi, A numerical method for the accurate solution of the Fokker-PlanckLandau equation in the nonhomogeneous case, J. Comput. Phys. 179 (1) (2002) 1 – 26. [20] J. S. Chang, G. Cooper, A pratical difference scheme for Fokker-Planck equations, J. Comput. Phys. 6 (1969) 1–16. [21] A. V. Bobylev, V. A. Chuyanov, On the numerical solution of Landau’s kinetic equation, USSR Comput. Math. Math. Phys. 16 (1976) 121–130. [22] A. B. Langdon, Conservative differencing of the electron Fokker-Planck transport equation, CECAM Report of Workshop on The Flux Limiter and Heat Flow Instabilities in Laser-Fusion Plasmas (1981) 69. [23] T. H. Kho, Relaxation of a system of charged particles, Phys. Rev. A 32 (1) (1984) 666–669. [24] E. M. Epperlein, Fokker-Planck modeling of electron transport in laser-produced plasmas, Laser Part. Beams 12 (1994) 257–272. [25] A. Arnold, A. Unterreiter, Entropy decay of discretized Fokker-Planck equations I — temporal semidiscretization, Comput. Math. Appl. 46 (10-11) (2003) 1683 – 1690. 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[26] R. J. Kingham, A. R. Bell, An implicit Vlasov-Fokker-Planck code to model non-local electron transport in 2-D with magnetic fields, J. Comput. Phys. 194 (2004) 1–34. [27] C. Buet, K.-C. Le Thanh, About positive, energy conservative and equilibrium state preserving schemes for the isotropic Fokker-Planck-Landau equation, CEA technical report hal00092543v2 (2006) 1–42. [28] E. M. Epperlein, R. J. Rickard, A. R. Bell, Two-dimensional nonlocal electron transport in laser-produced plasmas, Phys. Rev. Lett. 21 (1988) 2453–2456. [29] T. H. Kho, M. G. Haines, Nonlinear electron transport in magnetized laser plasmas, Phys. Fluids 29 (8) (1986) 2665–2671. [30] E. M. Epperlein, M. G. Haines, Plasma transport coefficients in a magnetic field by direct numerical solution of the Fokker-Planck equation, Phys. Fluids 29 (4) (1986) 1029–1041. [31] A. P. L. Robinson, A. R. Bell, R. J. Kingham, Fast electron transport and ionization in a target irradiated by a high power laser, Plasma Phys. Control. Fusion 48 (8) (2006) 1063–1076. [32] A. R. Bell, A. P. L. Robinson, M. Sherlock, R. J. Kingham, W. Rozmus, Fast electron transport in laser-produced plasmas and the KALOS code for solution of the Vlasov-Fokker-Planck equation, Plasma Phys. Control. Fusion 48 (3) (2006) R37–R57. [33] C. Buet, K.-C. Le Thanh, Positive,conservative, equilibrium state preserving and implicit difference schemes for the isotropic Fokker-Planck-Landau equation, CEA technical report hal-00142408v3 (2007) 1–26. [34] E. M. Epperlein, Implicit and conservative difference scheme for the Fokker-Planck equation, J. Comput. Phys. 112 (1994) 291–297. [35] M. Lemou, L. Mieussens, Fast implicit schemes for the Fokker-Planck-Landau equation, C. R. Acad. Sci. Paris, Ser. I 338 (10) (2004) 809 – 814. [36] L. Chac´on, D. C. Barnes, D. A. Knoll, G. H. Miley, An Implicit Energy-Conservative 2D Fokker-Planck Algorithm: II. Jacobian-Free Newton-Krylov Solver, J. Comput. Phys. 157 (1999) 654–682.

22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[37] X. Antoine, M. Lemou, Wavelet approximations of a collision operator in kinetic theory, C. R. Math. 337 (2003) 353–358. [38] H. Hang, N. M. Ghonien, Formulation of a moment method for multidimentional Fokker-Planck equations, Phys. Rev. E 51 (6) (1995) 5251–5260.

23