Diffusion and guiding center approximation for particle transport in strong magnetic fields

The diffusion limit of the linear Boltzmann equation with a strong magnetic field is performed. The giration period of particles around the magnetic field is assumed to be much smaller than the collision relaxation time which is supposed to be much smaller than the macroscopic time. The limiting equation is shown to be a diffusion equation in the parallel direction while in the orthogonal direction, the guiding center motion is obtained. The diffusion constant in the parallel direction is obtained through the study of a new collision operator obtained by averaging the original one. Moreover, a correction to the guiding center motion is derived.


Introduction
The motion of charged particles under the action of a strong magnetic field is an important phenomena encountered in artificial and natural plasmas (Tokamakas, ionospheric plasmas, etc). In the presence of constant magnetic field B, a charged particle (electron in our case) of mass m and charge q with velocity v follow a helical trajectory around magnetic field line. The gyration radius of the particles, called Larmor radius r L , is inversely proportional to the magnitude of the magnetic field (r L = mv qB ). This means that when the magnetic field becomes large, the particles get trapped along the direction of B. In addition, when an electric field E is applied, the particles experience a drift of the instantaneous centers of Larmor circles, called in the literature the guiding center motion, in the direction perpendicular to both the electric and magnetic fields. The velocity of this drift is given by the following formula Direct simulations of Vlasov or Boltzmann equations in the presence of such large magnetic fields requires the numerical resolutions of the small position and time scales induced by the gyration along the magnetic field. Hence, the question of deriving approximate models, numerically less expensive, is of great importance. Various models have been developed to this aim like the "guiding-center" and the gyrokinetic models. The guiding center approximation consists in averaging the motion over the gyroperiod when supposing B goes to +∞ (r L → 0) [26] and the gyrokinetic approximation is a generalization of the guiding center one in the case of kinetic equations. For a complete physical review, we refer for instance to [27,15,28,31]. Mathematically speaking, homogenization techniques were used in [19,20] by E. Frénod and E. Sonnendrücker to justify the high magnetic field limit of Vlasov and Vlasov-Poisson systems. It is proven in [19] that the tridimensional guiding center approximation leads to a one-dimensional kinetic model in the direction of the magnetic field. The drift phenomenon perpendicular to the magnetic field is rigourously obtained in [20] for the two dimensional Vlasov equation on a sufficiently long time scale.
Other asymptotic regimes of the Vlasov equation with strong magnetic fields are studied [21,18,22,16,17]. The gyrokinetic approximation of the Vlasov-Poisson system has also been intensively studied by F. Golse and L. Saint-Raymond for the two and tridimensional motions in various asymptotic regimes [24,33,25,34].
In all the above references, the transport is assumed to be ballistic, which means that particles do not suffer any collision during their motion. We are here interested in regimes were collisions are important and we propose to study the interplay between the fluid limit induced by collisions and the high frequency gyrations around the magnetic field. To simplify the study, we shall only consider collisions with a thermal bath at a given temperature. In this case, the diffusion limit of the obtained Boltzmann equation leads to the so-called drift diffusion equation. The question of deriving diffusions equations from the Boltzmann equation has been widely studied during the last three decades in the context of radiative transfer [5,6,35,1,7], in semiconductors and plasmas [32,23,30,2,3,4,8,9,10] and many other contexts that we do not mention here. In the presence of a magnetic field, an important quantity to be considered is the gyroperiod measuring the time that a particle of mass m with charge q and submitted to the constant magnetic field B makes a 2π rotation around this field (T c = 2πm/qB). When the gyroperiod is much larger than the relaxation time, then the magnetic field effect disappear in the diffusion limit as will be seen later on. When T c is of the same order of magnitude as the relaxation time, the diffusion matrix has an antisymmetric component, generated by the magnetic field (this has been proven for collision with walls by P. Degond & al [8,11,14]). Formal results for binary gas mixtures can also be found in [29,12,13,9]. We shall consider in this paper the situation where the gyroperiod is much smaller than the relaxation time. In that case, one expects that the diffusion matrix is obtained through the analysis of an averaged collision operator. Also, the high frequency oscillations only occur in the direction perpendicular to the magnetic field while parallel velocity stay unaffected. Therefore, the behaviours of the solutions in these directions are expected to be different. We shall show in this paper that the transport is diffusive along the magnetic field, while it is dominated by the guiding center drift in the orthogonal direction. The diffusion coefficient in the parallel direction will be obtained by averaging the collision cross section around the magnetic field.
The outline of the paper is as follows. In the next section, we introduce the scal-ing, notations and state the main results. The proofs rely on the study of the operator accounting for collisions and gyration around the magnetic field. The analysis of this operator is done in Section 3. Section 4 is devoted to the asymptotics of this operator when the scaled gyration period tends to zero. It is shown that this limit is well described by a collision operator with a cross section averaged around the magnetic field. This result is then used in Section 5 in order to prove the main results of the paper. Some concluding remarks are listed in the last section of the paper.

Setting of the problem and main results
The problem we are interested in is a singular perturbation of the three dimensional Boltzmann equation with a constant large magnetic field. The position variable is denoted by r = (x, y, z), the velocity variable is denoted v = (v x , v y , v z ). The magnetic field is assumed to be constant and parallel to the z axis. We shall denote by r ⊥ = (x, y) and v ⊥ = (v x , v y ), the orthogonal variables. The distribution function is a solution of the following singularly perturbed Boltzmann equation where T ⊥ and T z are the respectively perpendicular and parallel parts of the transport operator The collision operator is given by The function M is the normalized Maxwellian: and the cross section σ is symmetric and bounded from below and above by two positive constants. The parameter ε is intended to go to zero and represents the scaled mean free path. The scaled magnetic field is given by ez The gyroperiod is εη 2 and we shall consider the situation where η and ε go to zero simultaneously and the situation where ε tends to zero for η given then η tends to zero. Since the magnetic field strongly confines the motion of the particles in the perpendicular direction (x, y), we have rescaled this direction with the prefactor η. This is why the orthogonal operator T ⊥ comes with a different scaling from the parallel one T z . By doing so, we shall obtain a non trivial motion in the orthogonal direction (guiding center).

Scaling
In this subsection, we shall explain the hypotheses that lead to the scaled Boltzmann equation (2.1). In particular, we highlight the differences of the position scales in the parallel and perpendicular directions (with respect to the magnetic field). The starting point is the diffusion scaling in the absence of a magnetic field. The unscaled Boltzmann equation in this case reads where (t, X, V ) are the physical time position and velocity variables, while Φ is the electrostatic potential, e is the elementary charge and m the mass of the electron (to fix the ideas, we consider electrons). The collision operator Q has the form Σ is the transition rate and M T is the Maxwellian at the temperature T K B being the Boltzmann constant. The thermal energy and velocity are respectively defined by The electrostatic energy is assumed to be of the order of the thermal energy and varies over a macroscopic length scale L : The transition rate Σ is the inverse of a time. The typical value of this time, the scattering time, is denoted by τ s . So we write where σ is its scaled version assumed to be of order 1. The mean free path L s will be defined as the distance that a thermal electron crosses during the scattering time τ s This distance is assumed to be much smaller than the macroscopic distance L and we shall denote it by The only thing left to know is the macroscopic time scale. Since the kernel of the collision operator is generated by the Maxwellian which carries no flux, one has to go to the diffusion scaling and define the macroscopic time τ m by Using (τ m , L, V th , U th ) as units for time, position, velocity and energy, the Boltzman equation can be written Let us now look at the effect of a large magnetic field B = Be z . The scaled Boltzmann equation contains the additional term − q m (V × B) · ∇ V F. An electron submitted to this magnetic field spins around it following the equation The period of this precession, called the cyclotron period is given by We shall assume that this period is much smaller than the scattering time and we write Of course the motion along the magnetic field stays unchanged and we shall not change the length scale in this direction. In the orthogonal direction however, the combined action of the electric field and the magnetic field results in the drift of the guiding center (the electron motion is a rotation around a center which is now drifting under the action of the electric field). The velocity of this drift is given by Its order of magnitude is then given by wher E ⊥ is the perpendicular component of the electric field.
The last hypothesis that we make is the potential Φ varies in the orthogonal direction on a lengthscale L ⊥ which is not equal to the parallel lengthscale L. Moreover, L ⊥ is such that the guiding center moves during the macroscopic time by a distance of the order of L ⊥ itself. In other words, we want Replacing the thermal energy by its expression given above in this section, we find after straightforward computations that Now rescaling the orthogonal position variable by L ⊥ , we obtain the Boltzmann equation (2.1).
Remark 2.1. From the scaling hypotheses, we see that between two successive collisions, the electron precesses a large number of times around the magnetic field. Therefore, one might expect that the collision cross section will be averaged along these precessions. On the other hand, since the orthogonal lengthscale has been designed in order to see the guiding center drift, the limiting equation should exhibit this term. Finally, since the parallel motion to the magnetic field does not feel the magnetic field, one is entitled to expect a non vanishing diffusion in this direction. This is exactly what we shall prove in this paper : drift in the orthogonal direction, finite diffusion in the parallel direction governed by a collision operator averaged over precessions around the magnetic field.

Notations
The analysis of the limit ε → 0 naturally leads to the study of the operator while the limit η → 0, involves the cylindrical averaging around the vector e z . To this aim, we denote by R(τ ) the rotation around e z with angle τ which is the multiplication operator by the following matrix The cylindrical average of a function h = h(v) is defined by It is clear that A is the orthogonal projector (in the L 2 sense) on the set of cylindrically invariant functions. We also define the partial average We shall moreover define the averaged function σ and the symmetric average σ by and the corresponding collision operators We define the functional spaces and its cylindrically symmetric subspace We also define the space We shall make the following hypotheses Assumption 1. The cross-section σ belongs to W 2,∞ (R 6 ) and is supposed to be symmetric and bounded from above and below: We finally define the space (2.14)

Main results
The main results of this paper are summarized in the following two theorems. The first one deals with the limit ε → 0 while η > 0 is kept fixed.
Theorem 2.1. Let η > 0 be fixed. Let T ∈ R * + and assume that the initial data of (2.1), f 0 , belongs to L 2 M V . Then, with Assumptions 1 and 2, the problem (2.1) has a unique weak solution in and the diffusion matrix D η is given by the formula being the only solution of (3.7).
Remark 2.2. In the relaxation time case σ(v, v ′ ) = 1 τ , the matrix D η can be computed explicitly: In the limit η → 0, the upper 2 x 2 block of the matrix reduces to the antisymmetric matrix , whereas its symmetric part is of order η 2 . The following proposition shows that these features are still valid in the case of a non constant scattering section σ.
we have the following expansions in H 2 M as η tends to zero.
. The diffusion matrix D η is a definite positive matrix. Its symmetric and antisymmetric parts, D η s and D η as , have respectively the following expansions Remark 2.3. The above theorem shows that the diffusion is of order O(1) in the parallel direction while it scales like η 2 in the orthogonal one, exactly like in the relaxation time coefficient. The antisymmetric part of the diffusion matrix acts essentially on the orthogonal direction and leads to the guiding center motion. However, the above formula shows a correction of order η of the guiding center motion, which is due to collisions. Indeed, it is readily seen that for any three dimensional vector Z Theorem 2.2. Under the same hypotheses as for Theorem 2.1, we now let ε and η simultaneously and independently tend to zero (there is no assumption on the relative scale between ε and η). In this case, the sequence (f εη ) εη of weak solutions of (2.1) The parallel current density, J z is given by where D z is the diffusion constant given by and X z is the zero-th order term of X η z (defined by (4.7)).
3 Analysis of the operator Q η .
Let us consider the space L 2 M introduced in (2.10) and define the scalar product on this space by Let D(Q η ) be defined by We have the following result  Q η (f )dv = 0 (mass conservation).
3. The kernel of Q η is the real line spanned by the Maxwellian : 4. Let P be the orthogonal projection on KerQ η . The following coercivity inequality holds for any function f ∈ D where α 1 > 0 is the lower bound of σ (2.13).
5. The range of Q η , denoted by Im(Q η ), is the set of functions g ∈ L 2 M satisfying the solvability condition: In addition, for all g ∈ Im(Q η ) there exists a unique function f ∈ D satisfying both Q η (f ) = g and Proof. First of all, we recall that the operator Q satisfies all items of Proposition 3.1 (by substituting D by L 2 M ). Now Items 1., 3. and 4. of the Proposition follow immediately from these properties since we have In order to prove Item 2., we need to prove that I − Q η is one-to-one from D to L 2 M . In order to do so, we rewrite the equation f − Q η (f ) = g under the following form The characteristics of the equation are defined by The solution of the equation can then be defined by integration over the characteristics and finally leads to A fixed point argument leads to the existence and uniqueness of f ∈ L 2 M solving the equation f − Q η (f ) = g in the distributional sense. Since Q is a bounded operator on L 2 M , we immediately obtain that (v × e z ) · ∇ v f ∈ L 2 M so that the constructed solution is in the domain D of the operator Q η . Let us now prove item 5. Consider g ∈ L 2 M such that R 3 g(v)dv = 0. From the maximal monotonicity of Q η , for each λ > 0, λId − Q η is surjective. Therefore, for all λ > 0, there exists f λ ∈ D(Q η ) satisfying Integrating this equation with respect to v gives λ Taking the scalar product of (3.6) with f λ in L 2 M , one obtains This implies that (f λ ) λ>0 is bounded in L 2 M . Up to an extraction of a subsequence, there exists f ∈ L 2 M such that (f λ ) converges weakly to f . By passing to the limit λ → 0 in the weak form of equation (3.6), we get In addition, Q(f ) and g belong to L 2 M , then (v × e z ) · ∇ v f ∈ L 2 M and then f ∈ D(Q η ). Of course R 3 f dv = 0 since this is true for f λ . Now if there is another solution of −Q(·) = g with zero velocity average, the difference of this solution with f is in the kernel of Q η and has a zero average. It is necessarily equal to zero.
The following corollary is a direct consequence of Proposition 3.1.
and from the coercivity inequality (3.2), (X η z ) η and (η 2 X η ⊥ ) η are bounded sequences in L 2 M . Namely, we have the following a priori estimates The following proposition allows to reformulate the problem Q η (f ) = −g by using the characteristics.
Proposition 3.2. For all g ∈ Im(Q η ), we have the following equivalence where Q + is given by (3.4) and L η is the operator on L 2 M given by the expression: Proof. Let (S η , D(S η )) be the unbounded operator on L 2 M defined by With the decomposition (3.3), we have Q η = Q + − S η . Proceeding as in the proof of Proposition 3.1 (replacing ν + 1 by ν), we show that S η is invertible and its inverse is which concludes the proof.
Proof. This estimate follows immediately from the definition of L η (3.10) using Assumption 1.

Proposition 3.3.
Under Assumption 1, the solution X η of (3.7) belongs to (H 2 M (R 3 )) 3 and X η and all its derivatives with respect to v of order less than or equal 2 (up to dividing by the Maxwellian) are polynomially increasing when |v| goes to +∞. Namely, we have for all α = (α 1 , α 2 ) ∈ N 2 with α 1 + α 2 ≤ 2 and where ∂ α : Here, Q η 1 and Q η 3 are two polynomials of degrees 1 and 3 respectively (depending on η).
Proof. From (3.9), we deduce the identity and with (3.13) one obtains (3.14). Estimate (3.15) follows by derivation of the above equation with respect to v and by using the fact that the cross-section σ(v, v ′ ) belongs to W 2,∞ (R 6 ).

Expansion of X η with respect to η
A rigourous expansion of X η around η = 0 will be carried out in this section. Corollary 3.1 provides us with bounds of order O(1) for X η z and of order O( 1 η 2 ) for X η ⊥ . We shall prove in this section that X η ⊥ is actually bounded in H 2 M . This is due to the fast precession around the axis e z generated by the differential term 1 η 2 (v × e z ) · ∇ v . Filtering out these oscillations is done through the reformulation (3.9) of the equation satisfied by X η . Moreover, we shall expand X η in powers of η 2 up to the first order. We now list the properties of Q which are completely inherited from those of Q.  3. The null set of Q is given by be the orthogonal projection on N(Q), the following coercivity inequality holds, 5. The range of Q in L 2 M is given by In the sequel, we shall also need the properties of Q considered on the whole space L 2 M (whereas Q is restricted to cylindrically symmetric functions). Some of the properties of Q are inherited from those of Q but some of these properties like the symmetry are lost. The following proposition summarizes these properties.

The range of Q in L 2
M is the following set 3. For each g ∈ Im(Q), there exists a unique function f ∈ L 2 M satisfying such that Q(f ) = g. The solution can be written

2)
and where we have denoted Q = Q + − ν. In addition, there exists a constant C > 0 independent of g (and f ) such that Proof. Let us begin with item 1. It is clear that if f is Maxwellian, it is in the kernel of Q. Let now f ∈ N(Q). We have Q + (f ) − νf = 0. Since the Maxwellian is cylindrically symmetric, it is readily seen that Q + (f ) is always cylindrically symmetric. Therefore,

is cylindrically symmetric and thus
In view of Proposition 4.1, f is proportional to the maxwellian.
In order to prove items 2. and 3., we consider g ∈ Im(Q), and f ∈ L 2 M such that Q(f ) = g. Decomposing f and g as The left hand side of the above equality is cylindrically symmetric, whereas cylindrical averages of the right hand side vanish. Therefore, both are equal to zero. We deduce that which can be rewritten Q(f ) = g − Q + ( g−g ν ) and also In addition, since h f = g − g ν and with (2.13), there exists C 2 > 0 such that h f M ≤ C 2 g M and estimate (4.3) holds.
Definition 4.2. We define the gyration operator G by The kernel of this operator is nothing but the set of cylindrically symmetric functions. It is clear that its range is contained in the set of L 2 M functions which have zero cylindrical averages. The following lemma shows that both sets coincide. in nonempty. Any solution f can be written where f 1 si an arbitrary cylindrically symmetric function (which means that A(f 1 ) = f 1 ) and the average operator A 1 is defined by (4.5) Moreover, we have A(νA 1 (g)) = 0.
Proof. Let us consider the problem This problem has a unique solution f η defined by where C η (v) = e η 2 2π 0 ν(R(s)v)ds − 1 = e 2πη 2 ν − 1. Passing to the limit η → 0, we find that f η converges almost everywhere towards A 1 (g). Moreover A(η 2 νf η ) = A(g) = 0, which proves after passing to the limit η → 0 that A(νA 1 (g)) = 0. Of course, any other solution of the equation G(f ) = g is obtained by adding an arbitrary cylindrically symmetric function.
Let us now consider the problem with f 1 is an arbitrary cylindrically symmetric function. Moreover, the mapping (g, h) → f is linear and continuous both on First of all, Lemma 4.1 shows that f = f 1 + A 1 (g). Besides, we have Since f 1 is cylindrically symmetric, then A(νf 1 ) − νf 1 = 0, which leads, thanks to the identity A(νA 1 (g)) = 0, to Q(f ) = −νA 1 (g) + h. It is now enough, in view of the solvability condition (4.1) to check that ν ν (−νA 1 (g) + h)dv = 0, which is readily seen.

Expansion of X η z
The idea is to make a Hilbert expansion of X η z in powers of η 2 : We have the following equations z )) = 0, i ≥ 1. Therefore, the terms of the expansion can be successively computed by solving the problems As for X η z , the above system can be reformulated as follows Proposition 4.5. The expansion The proof of this proposition is identical to that of Proposition 4.4 and is skipped. The only thing which has to be checked is the expression of the anisotropic part of X

Proof of the main theorems
We begin this section by making precise the definition of weak solutions of the Boltzmann equation and give quite standard a priori estimates on this solution. For the moment, ε and η are arbitrary and all the constants are, unless specified, independent of these parameters.
Theorem 5.1. Assume that f 0 ∈ L 2 M V (drdv). Then, ∀ε, η > 0 and ∀T ∈ R * + , there exists a unique weak solution f εη ∈ C 0 ([0, T ], L 2 M V (drdv)) of (2.1). In addition, the following mass conservation equation holds and there exists a constant C > 0 independent of ε and η such that An immediate consequence of the integrability properties of the function f εη , with respect to v is that the weak formulation (5.1) is satisfied by test functions which are not necessarily compactly supported in velocity. More precisely, we have We shall not prove this corollary which is straightforward. The proof of Theorem 5.1 is also classical (see for example [32]). We show here how the a priori estimates (5.4) and (5.5) can be obtained. We simply multiply the Boltzmann equation (2.1) by f εη M V and integrate with respect to r and v. Straightforward computations show that this leads to the entropy inequality d dt Thanks to the boundedness of the potential V and its derivatives, estimate (5.4) follows by applying Gronwall Lemma and then the bound (5.5) follows immediately.

Proof of Theorem 2.1
In this section, we are interested in the limit ε → 0, while η is fixed. To this aim, we follow the usual approach based on the moment method. Namely, thanks to estimates (5.4) and (5.5), the sequences f εη , (ρ εη ), (J εη z ) and (J εη ⊥ ) are bounded with respect to ε. One can find ρ η (t, x) ∈ L 2 ((0, T ) × R 3 ) such that ρ εη converges weakly to ρ η in L 2 ((0, T ) × R 3 ) and f εη ⇀ ρ η M in L ∞ (0, T ; L 2 M (R 6 )). Similarly, J εη converges weakly to J η in L 2 ((0, T )×R 3 ). We have of course, the mass conservation equation ∂ t ρ η + div r J η = 0. Now the only thing left to show is to identify J η . The way to proceed for identifying J η is standard (see for instance [1,5,35,7,32,10,2,3,4]). The Boltzmann equation can be rewritten where D (1) It is clear that D Let us now expand D η z . We have The first term of the right hand side converges as ε and η tend to zero towards