1 Introduction
The progress in computational technology made in recent decades has stimulated the use of the Monte Carlo approach in many areas of research. Plasma physics is not an exception. In plasma theory, the Monte Carlo method is used to simulate the random motion of test charged particles and to extrapolate these results to sample ensembles and bulk plasma behaviour. The basis of this approach is the interpretation of the kinetic equations for electron and ion distribution functions, including the operators of pair Coulomb collisions in the Fokker–Planck (FP) form, as a representation of the forward Kolmogorov equations for transition probabilities of the Markov single-particle trajectories in the space of coordinates and velocities,
$X(t)=(\boldsymbol{r}(t),\boldsymbol{v}(t))$
. Such a description mathematically justifies the use of stochastic differential equations (SDEs) in the form initiated by Ito (Reference Ito1944), which is now generally accepted. Today, the literature on this topic is quite rich. We will refer to the first monograph reprinted in Gikhman & Skorohod (Reference Gikhman and Skorohod1972) where SDE theory was presented in complete form, and a comprehensive textbook (Kloeden & Platen Reference Kloeden and Platen1995). All the coefficients of the Ito stochastic differential are strictly defined by the coefficients of the FP operator; there is no problem writing down the Ito SDEs for single-particle motion in the most general form in which the diffusion in the velocity space is determined by the three-dimensional (3-D) Ito process. However, these equations are too complicated for a numerical solution, but with reference to magnetic confinement they can be reduced to a 5-D set of drift SDEs to provide a description of the guiding centre diffusion in the configuration space. In this paper we propose such a transformation.
In a magnetized plasma, the transition to the reduced (guiding centre) phase space is well studied in the drift theory of single-particle motion in an inhomogeneous magnetic field in the absence of collisions. Drift equations are obtained by the averaging method in theory of oscillations (Bogoliubov Reference Bogoliubov and Mitropolski1961). The drift analysis uses a small parameter
$\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D70C}/L$
, where
$\unicode[STIX]{x1D70C}$
is the (Larmor) gyro-radius,
$L$
is the inhomogeneity scale and the dynamics of
$\unicode[STIX]{x1D707}(t)$
,
$\unicode[STIX]{x1D700}(t)$
(the magnetic moment and kinetic energy of particle) is determined only to zero order with respect to parameter
$\unicode[STIX]{x1D6FF}$
together with the one-dimensional motion along the magnetic force lines, while the drift across magnetic field is defined in the first order with respect to
$\unicode[STIX]{x1D6FF}$
. The mixing of orders in the description of the spatial motion creates difficulties in the kinetic theory, especially with taking into account the spatial diffusion, because the coefficients of collisional diffusion bring into the FP operator terms of the second order in
$\unicode[STIX]{x1D6FF}$
according to the drift theory.
This led to the emergence of two approaches to the description of spatial diffusion. In the Longmire and Rosenbluth paper (Longmire & Rosenbluth Reference Longmire and Rosenbluth1956) the kinetic theory of diffusion of a fully ionized plasma in a strong uniform magnetic field is presented. This result was obtained by transforming the FP collision operator in to 6-D space
$X=(\boldsymbol{R},\boldsymbol{v})$
with gyro-averaging being performed when calculating the classical diffusion coefficients,
$D_{c}=\unicode[STIX]{x1D708}\unicode[STIX]{x1D70C}^{2}$
, under condition
$\unicode[STIX]{x1D708}\ll \unicode[STIX]{x1D6FA}$
, where
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D6FA}$
are the collisional and cyclotron frequencies, defined as
$\unicode[STIX]{x1D6FA}=eB/(mc)$
. Later, in the neoclassical approach for axisymmetric magnetic configurations, simplified gyro-averaged kinetic equations are formulated in the drift 5-D space
$\bar{X}=(\boldsymbol{R},\unicode[STIX]{x1D707},\unicode[STIX]{x1D700})$
where the spatial cross-field drift is taken into account to the first order in
$\unicode[STIX]{x1D6FF}$
, while the kinetics of
$\unicode[STIX]{x1D707}$
,
$\unicode[STIX]{x1D700}$
and the FP collision operator defined in the space
$\bar{X}=(\boldsymbol{R},\unicode[STIX]{x1D707},\unicode[STIX]{x1D700})$
for the gyro-tropic background plasma considered in zero approximation with respect to
$\unicode[STIX]{x1D6FF}$
. Thus, classical diffusion was excluded from consideration, but the collisional spatial diffusion of trapped and passing particles is formulated in the neoclassical drift time scale that allows averaging over the periods of drift particle motion. In many papers this approach is used to describe the relatively large time scale diffusion of the closed orbits of guiding centres circulating around the magnetic axis or trapped in the region of a weak magnetic field in the toroidal tokamak plasma due to rare pair collisions. In this way, the banana-plateau diagram is determined for the neoclassical diffusion coefficient,
$D_{neo}=\unicode[STIX]{x1D708}^{\ast }\unicode[STIX]{x1D70C}^{2}$
where
$\unicode[STIX]{x1D708}^{\ast }$
the collisional frequency is renormalized for various neoclassical diffusion regimes in Hinton & Hazeltine (Reference Hinton and Hazeltine1976) and Galeev & Sagdeev (Reference Galeev, Sagdeev and Leontovich1979). In the papers of Xu & Rosenbluth (Reference Xu and Rosenbluth1991), Eriksson & Helander (Reference Eriksson and Helander1994), Tessarotto, White & Zheng (Reference Tessarotto, White and Zheng1994), Boozer (Reference Boozer2002) and others, the coefficients of the spatial FP operator were determined, which were called the Monte Carlo operators, and model SDEs were formulated in 3-D space
$(\bar{r},\unicode[STIX]{x1D707},\unicode[STIX]{x1D700})$
where
$\bar{r}$
marks the guiding orbit-averaged banana radial parameter without references to the theory of diffusion processes.
Another approach is realized by Brizard (Reference Brizard2004), where a complete representation of the gyro-kinetic equation in the form of a five-dimensional Fokker–Planck operator is obtained directly by means of the Lie transform of the original Fokker–Planck operator in the phase
$X=(\boldsymbol{R},\boldsymbol{v})$
into the drift space
$\bar{X}=(\boldsymbol{R},\unicode[STIX]{x1D707},\unicode[STIX]{x1D700})$
, that makes it possible to exclude a fast variable similar to how it is realized within the Bogolyubov–Mitropol’skii method in the theory of single-particle motion (Bogoliubov Reference Bogoliubov and Mitropolski1961). Thus, the work (Brizard Reference Brizard2004) can be considered as an expansion of the approach early proposed by Longmire & Rosenbluth (Reference Longmire and Rosenbluth1956). The results (Brizard Reference Brizard2004) were implemented in the Monte Carlo method in the form of a complete set of SDEs for the variables
$\bar{X}$
(Hirvijoki et al.
Reference Hirvijoki, Brizard, Snicker and Kurki-Suoni2013). In this context, the paper of Chang, Hong & Chenhao (Reference Chang, Hong and Chenhao2011) should be also mentioned, where the 5-D FP operator was derived with coefficients directly calculated accounting for the effect of Coulomb collisions on the drift variables in a Lorenz plasma.
In this paper, we propose an alternative derivation of Brizard’s SDEs, which does not use any transformations of the kinetic equations in the 6-D space,
$X$
, and does not use the reduced 5-D FP operator in general. We obtain the drift SDEs from the most general non-gyro-averaged SDEs formulated in the phase space
$X$
by applying the Ito formula to calculate the stochastic differentials of functions of random variables when they determine a nonlinear transformation of the diffusive components. Then, gyro-averaged SDEs are finally obtained by calculating stochastic integrals in the basic drift time scale. This procedure reduces to a simple gyro-averaging. Proposed approach provides simpler and more effective algorithm to calculate the drift stochastic particle motion naturally combining the results of the classical drift theory and classical diffusion description. It should be pointed out, that in proposed approach both the drift and random parts of the stochastic differentials are determined in the same order of the drift theory (preliminary results were published in Gurin & Yavorskij (Reference Gurin and Yavorskij2017)). In the present paper, some coefficients of Hirvijoki et al. (Reference Hirvijoki, Brizard, Snicker and Kurki-Suoni2013) are reduced to the relevant orders of the drift theory. We give a full set of SDEs, which can serve as a tool for simulating infinite drift diffusion trajectories of the guiding centres that have escaped from the toroidal plasma without using analytical models of banana and transit orbits. Finally, we present SDEs adopted for the practical solutions in toroidal quasi-cylindrical coordinates in axisymmetric magnetic configurations. These equations are applicable for the description of both bulk plasma and ion spices such as injected ions or charged fusion products (excluding their self-collisions).
The paper is organized as follows. In § 2 we present the basic equations unifying the SDEs and plasma kinetic theories. In § 3 we formulate the SDEs in drift approximation in velocity space for gyro-tropic magnetized plasmas in accordance with the kinetic theory of fully ionized plasma. In § 4 the general formulations for the commonly used approximation of isotropic plasmas are provided. In § 5 the SDEs describing diffusion of guiding centres in physical space are presented. For the case of axial symmetry of the magnetic configuration, we propose a model for spatial diffusion exposing the random Brownian trapped/passing walk of the guiding centres in poloidal plane. This model, being numerically implemented, can serve as an explicit demonstration of the neoclassical ‘banana’ diffusion in tokomaks. In § 6 we summarize the results and discuss some features of the derived equations.
2 Main formulas of stochastic and kinetics
The Ito SDE for Markov multidimensional diffusion process
$X(t)=(x_{1}(t),\ldots ,x_{n}(t))$
, is
$$\begin{eqnarray}\text{d}x_{i}=a_{i}(t,X(t))\,\text{d}t+\mathop{\sum }_{j=1}^{m}b_{ij}(t,X(t))\,\text{d}W_{j}(t),\end{eqnarray}$$
where
$W_{j}(t)$
is the Wiener processes, and the forward Kolmogorov equation for the transition probability of Markov continuous process is
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}p_{s,Y}(t,X)+\mathop{\sum }_{i=1}^{n}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{i}}[a_{i}(t,X)p_{s,Y}(t,X)]=\frac{1}{2}\mathop{\sum }_{i,j=1}^{n}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x_{i}\unicode[STIX]{x2202}x_{j}}[d_{ij}(t,X)p_{s,Y}(t,X)],\end{eqnarray}$$
for
$t>s$
and the initial condition
$p_{s,Y}(s,X)=\unicode[STIX]{x1D6FF}(X-Y)$
, represent the alternative approaches to the complete description of the Markov process (Gikhman & Skorohod Reference Gikhman and Skorohod1972; Kloeden & Platen Reference Kloeden and Platen1995). The coefficients
$a,b$
are determined at the left edge of the time interval (
$t,t+\text{d}t$
) where the independence of the Wiener increments,
$\text{d}W_{i}(t)=W_{i}(t+\text{d}t)-W_{i}(t)$
, provides the Markovian solution of SDE. The particularity of Ito’s analysis is the possibility of considering the quadratic infinitesimal form
$\text{d}W_{i}\,\text{d}W_{j}$
as a non-random quantity:
$$\begin{eqnarray}\text{d}W_{i}\,\text{d}W_{j}=\unicode[STIX]{x1D6FF}_{ij}\,\text{d}t.\end{eqnarray}$$
In the same sense, with probability 1, the following equations take place:
$$\begin{eqnarray}\text{d}x_{i}\,\text{d}x_{j}=\mathop{\sum }_{k=1}^{m}b_{ik}b_{jk}\,\text{d}t=\text{d}_{ij}\,\text{d}t.\end{eqnarray}$$
Equation (2.4) determines the correspondence between the diffusion coefficients
$b_{ij}$
and
$d_{ij}$
in the Ito and the Kolmogorov equations. The coefficients
$a_{i}$
are strictly equal in (2.1) and (2.2).
Equations (2.3) and (2.4) provide the rule for calculating the stochastic differential of an arbitrary function
$g(X)$
, known as the Ito formula:
$$\begin{eqnarray}\text{d}g(X)=\mathop{\sum }_{i=1}^{n}\,\text{d}x_{i}\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}x_{i}}+\frac{\text{d}t}{2}\mathop{\sum }_{i=1}^{n}\mathop{\sum }_{j=1}^{n}d_{ij}\frac{\unicode[STIX]{x2202}^{2}g}{\unicode[STIX]{x2202}x_{i}x_{j}},\end{eqnarray}$$
where the stochastic differential (2.1) should be substituted in (2.5). The Ito formula plays an important role in analysis allowing us to enter new variables into description.
The parameters
$m$
and
$n$
are arbitrary in (2.1), usually
$m<n$
. Also, not all of components
$\text{d}x_{i}$
can contain a fluctuation Wiener part. Such a situation is realized in the theory of fully ionized plasma in which the spatial trajectories of particles are taken to be differentiable, while dynamics in velocity space is considered as a random diffusion process due to the pair Coulomb collisions of charged particles. The standard kinetic description of a fully ionized plasma is formulated in terms of the electrons and ions distribution functions
$f_{\unicode[STIX]{x1D6FC}}(t,\boldsymbol{r},\boldsymbol{v})$
given by the system of equations (we use the non-index notation of the three-dimensional vector algebra):
Here
$l^{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}}$
is a Coulomb logarithm,
$\unicode[STIX]{x1D711}$
and
$\unicode[STIX]{x1D713}$
been the Rosenbluth potentials. With some redesignations (Rosenbluth, MacDonald & Judd Reference Rosenbluth, MacDonald and Judd1957) we have:
$$\begin{eqnarray}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FD}}(\boldsymbol{v})=\int \,\text{d}\boldsymbol{v}^{\prime }|\boldsymbol{v}-\boldsymbol{v}^{\prime }|f_{\unicode[STIX]{x1D6FD}}(\boldsymbol{v}^{\prime }),\quad \unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}}(\boldsymbol{v})=\int \text{d}\boldsymbol{v}^{\prime }f_{\unicode[STIX]{x1D6FD}}(\boldsymbol{v}^{\prime })/|\boldsymbol{v}-\boldsymbol{v}^{\prime }|.\end{eqnarray}$$
In accordance with definition of the distribution function as a density of a particles in the phase space,
$X=(\boldsymbol{r},\boldsymbol{v})\,f_{\unicode[STIX]{x1D6FC}}$
can be expressed through the single-particle transition probability
$p_{s,Y}(t,X)$
if trajectories
$X(t)$
are Markov processes:
$$\begin{eqnarray}f_{\unicode[STIX]{x1D6FC}}(t,X)=\int \,\text{d}Y\,f_{\unicode[STIX]{x1D6FC}}(s,Y)p_{s,Y}(t,X).\end{eqnarray}$$
So, the distributions
$f_{\unicode[STIX]{x1D6FC}}$
and
$p$
turn out to be solutions of the forward Kolmogorov equation, corresponding to their normalizations:
$\int \text{d}Xp_{s,Y}(t,X)=1$
,
$\int \text{d}X\,f_{\unicode[STIX]{x1D6FC}}(t,X)=N_{\unicode[STIX]{x1D6FC}}$
, where
$N_{\unicode[STIX]{x1D6FC}}$
is the full number of particles of the
$\unicode[STIX]{x1D6FC}$
species in the plasma volume. The Markov interpretation gives a broader meaning to the value
$p$
than a fundamental solution of the kinetic equation (2.6), which is suitable only for linear problems whereas collision operator (2.6) introduces a nonlinearity into the kinetic theory, if even the magnetic and electric fields are assumed to be known. Equation (2.9) can serve as the basis for applying an iterative procedure for determining the quasi-equilibrium distribution functions, calculating step-by-step the transition probability
$p$
for an ensemble of particles with test initial distributions and potentials by the Monte Carlo method. However, if single-particle trajectories against the background of plasma with fixed distributions are of interest, the solution may be obtained by one-step Monte Carlo procedure. All the information required for the Monte Carlo method is provided by the full SDE system that can be written according to the definitions (2.1), (2.2), (2.4), (2.6), (2.7) as:
Here the symbol
$T$
denotes transposition, since the matrix 
may not be symmetric.
3 Diffusion of velocity in gyro-tropic plasma
The SDEs (2.10) look very laconic, but are not easy to solve numerically. In the general case, the factorization of the
$3\times 3$
diffusion matrix 
according to (2.11) turns out to be a laborious task. In the plasma confined in a strong magnetic field, the distribution functions and, consequently, the Rosenbluth potentials depend on two velocity components. The analytic transformations of SDEs are possible by means of Ito’s formula when new variables are introduced. In this section the factorization of diffusion matrix will be carried out analytically. As well the stochastic description of test particle motion will complete the standard drift theory of the collisionless motion of charged particles.
As in the drift theory (Bogoliubov Reference Bogoliubov and Mitropolski1961), the parallel and transverse velocity components in a magnetic field
$\boldsymbol{B}(\boldsymbol{r})$
may be expressed as
$u=\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{h}$
,
$w=|\boldsymbol{v}-u\boldsymbol{h}|$
, where
$\boldsymbol{h}=\boldsymbol{B}/B$
(the possible non-stationary of the field is neglected below). In the velocity space with a positive local triplet
$(\boldsymbol{h}(\boldsymbol{r}),\boldsymbol{n}_{1}(\boldsymbol{r}),\boldsymbol{n}_{2}(\boldsymbol{r}))$
, the vector
$\boldsymbol{v}$
and operator
$\unicode[STIX]{x1D735}_{v}$
take the form, which explicitly reflects the dependence on the gyro-angle
$\unicode[STIX]{x1D6FC}$
:
$\boldsymbol{v}=u\boldsymbol{h}+w\boldsymbol{e}_{w}$
,
$\boldsymbol{e}_{w}=\boldsymbol{n}_{1}\cos \unicode[STIX]{x1D6FC}+\boldsymbol{n}_{2}\sin \unicode[STIX]{x1D6FC}$
,
$\unicode[STIX]{x1D735}_{v}=\boldsymbol{h}(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}u)+\boldsymbol{e}_{w}(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}w)+\boldsymbol{e}_{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x2202}/w\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})$
,
$\boldsymbol{e}_{\unicode[STIX]{x1D70C}}=\boldsymbol{h}\times \boldsymbol{e}_{w}=(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})\boldsymbol{e}_{w}$
. The SDEs for
$u,w$
can be obtained by using the Ito rules (2.5). According to the full set of (2.10), the differential
$\text{d}\boldsymbol{r}$
does not contain the Wiener part, and variable
$u$
depends linearly on
$\boldsymbol{v}$
, so
$u$
is differentiated according to the ordinary rules, whereas the variable
$w$
introduces nonlinearity and is differentiated by Ito’s formula. So, the differentials
$\text{d}u$
,
$\text{d}w$
take the form:
In explicit form, the second-order derivatives
$\unicode[STIX]{x1D735}_{\boldsymbol{v}}\unicode[STIX]{x1D735}_{\boldsymbol{v}}w=\boldsymbol{e}_{\unicode[STIX]{x1D70C}}\boldsymbol{e}_{\unicode[STIX]{x1D70C}}/w$
give the additional
$\text{d}t$
-contribution into second equation (3.1):
To obtain the drift version of these equations, it suffices to consider the stochastic integrals defining increments
$\unicode[STIX]{x0394}u=\int _{t}^{t+\unicode[STIX]{x0394}t}\,\text{d}u$
,
$\unicode[STIX]{x0394}w=\int _{t}^{t+\unicode[STIX]{x0394}t}\,\text{d}w$
within the drift scale time interval which includes many gyro-periods,
$\unicode[STIX]{x1D6FA}\unicode[STIX]{x0394}t\gg 1$
but it is relatively small in a time scale of variations of differential coefficients in (3.2),
$\unicode[STIX]{x1D708}\unicode[STIX]{x0394}t\ll 1$
, where
$\unicode[STIX]{x1D708}$
is estimated through
$a$
-terms in (2.11) as
$\unicode[STIX]{x1D708}=|a/v|$
. In the collisionless case, this procedure leads to substituting the value
$\langle \boldsymbol{e}_{w}\boldsymbol{\cdot }(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{h}\rangle =w\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{h}/2$
(brackets
$\langle \ldots \rangle$
will denote averaging over
$\unicode[STIX]{x1D6FC}$
). The diffusion terms in (3.2), being represented on the triplet
$(\boldsymbol{e}_{u},\boldsymbol{e}_{w},\boldsymbol{e}_{\unicode[STIX]{x1D70C}})$
,
$\boldsymbol{e}_{u}=\boldsymbol{h}$
and do not depend on
$\unicode[STIX]{x1D6FC}$
thus manifest the result of averaging with only some redefinition of the Wiener increments. The three-component Wiener differential
$\text{d}\boldsymbol{W}$
, defined with indices
$i=(u,w,\unicode[STIX]{x1D70C})$
in the basis
$\boldsymbol{e}_{i}$
instead the original triplet
$(\boldsymbol{h},\boldsymbol{n}_{1},\boldsymbol{n}_{2})$
turns out the same standard three-component Wiener differential. It is easy to see that the components
$(\text{d}W^{u},\text{d}W^{w},\text{d}W^{\unicode[STIX]{x1D70C}})$
obey the same algebra as in (2.3) regardless of the value
$\unicode[STIX]{x1D6FC}$
:
$$\begin{eqnarray}\text{d}W^{i}\,\text{d}W^{j}=(\boldsymbol{e}_{i}\boldsymbol{\cdot }\,\text{d}\boldsymbol{W})(\boldsymbol{e}_{j}\boldsymbol{\cdot }\,\text{d}\boldsymbol{W})=(\boldsymbol{e}_{i}\boldsymbol{\cdot }\,\boldsymbol{e}_{j})\,\text{d}t=\unicode[STIX]{x1D6FF}_{ij}\,\text{d}t.\end{eqnarray}$$
In gyro-tropic plasma, the distribution functions and Rosenbluth potentials depend only on
$u,w$
; therefore, the vector notation of the diffusion tensor (2.7) takes the form:
It can be represented by a block-diagonal matrix:
The diagonal block
$2\times 2$
can easily be decomposed into two identical symmetric matrix multipliers by various methods. Hear we use the algebra of Pauli matrices and reduce the problem of determining the elements of
$b$
-block to the solution of a biquadratic equation. The result can be formulated using auxiliary quantities
$d_{0}=(d_{uu}d_{ww}-d_{uw}^{2})^{1/2}$
,
$b_{0}=(2d_{0}+d_{uu}+d_{ww})^{1/2}$
as
$$\begin{eqnarray}b_{uu}=(d_{0}+d_{uu})/b_{0},\quad b_{ww}=(b_{0}+d_{ww})/b_{0},\quad b_{uw}=b_{wu}=d_{uw}/b_{0},\quad b_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}=\sqrt{d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}}.\end{eqnarray}$$
One can verify that these solutions satisfy (2.10), i.e.
$$\begin{eqnarray}d_{uu}=b_{uu}^{2}+b_{uw}^{2},\quad d_{ww}=b_{ww}^{2}+b_{uw}^{2},\quad d_{uw}=b_{uw}(b_{uu}+b_{ww}).\end{eqnarray}$$
Thus, we obtain approximate equations for the stochastic increments
$\unicode[STIX]{x0394}u$
,
$\unicode[STIX]{x0394}w$
that preserve the Ito’s differential structure with
$\unicode[STIX]{x0394}t$
-and
$\unicode[STIX]{x0394}W$
-terms as differentials. They do not contain the parameter
$\unicode[STIX]{x1D6FC}$
and can be regarded as infinitesimal in the drift time scale with the symbol substitution
$\unicode[STIX]{x0394}\rightarrow d$
:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{d}u=\left({\displaystyle \frac{e}{m}}\boldsymbol{h}\boldsymbol{\cdot }\boldsymbol{E}+{\displaystyle \frac{w^{2}}{2}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{h}+a_{u}^{(c)}\right)\,\text{d}t+b_{uu}\,\text{d}W_{u}+b_{uw}\,\text{d}W_{w},\\ \text{d}w=\left(-{\displaystyle \frac{uw}{2}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{h}+a_{w}^{(c)}+{\displaystyle \frac{1}{2w}}d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}\right)\,\text{d}t+b_{wu}\,\text{d}W_{u}+b_{ww}\,\text{d}W_{w}.\end{array}\right\}\end{eqnarray}$$
Note that
$\text{d}\boldsymbol{r}=\langle \boldsymbol{v}\rangle \,\text{d}t=u\boldsymbol{h}\,\text{d}t$
in this approximation.
The drift approach is usually formulated in terms of constants of drift motion, i.e. the adiabatic invariant
$\unicode[STIX]{x1D707}=w^{2}B_{0}/2B$
, where
$B_{0}$
is an arbitrarily chosen magnetic field magnitude, and
$\unicode[STIX]{x1D700}=v^{2}/2$
. In Ito’s formulation, it follows:
$\text{d}B=u\boldsymbol{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\,\text{d}t$
,
$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D707}=\frac{B_{0}}{B}\left[\left(wa_{w}^{(c)}+\frac{1}{2}(d_{ww}+d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}})\right)\,\text{d}t+w(b_{wu}\,\text{d}W_{u}+b_{ww}\,\text{d}W_{w})\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{d}\left(\unicode[STIX]{x1D700}+\frac{e}{m}U\right) & = & \displaystyle \left[\left(ua_{u}^{(c)}+wa_{w}^{(c)}+\frac{1}{2}(d_{uu}+d_{ww}+d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}})\right)\,\text{d}t\right.\nonumber\\ \displaystyle & & \displaystyle +\left.(ub_{uu}+wb_{wu})\,\text{d}W_{u}+(ub_{uw}+wb_{ww})\,\text{d}W_{w}\vphantom{\left(\frac{1}{2}\right)}\right].\end{eqnarray}$$
Here the potentiality of the electric field
$\boldsymbol{E}=-\unicode[STIX]{x1D735}U,\text{d}U=u\boldsymbol{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}U\,\text{d}t$
, is assumed. All the coefficients in these equations are explicitly expressed in terms of the Rosenbluth potentials according to the definitions (2.6), (2.7), (3.5), and (3.6). Final formulas for the
$\text{d}t$
-terms in (3.9), (3.10) are simplified with considering the equality
$\unicode[STIX]{x0394}_{v}\unicode[STIX]{x1D711}=2\unicode[STIX]{x1D713}$
, with
$\unicode[STIX]{x0394}_{v}(\ldots )=\unicode[STIX]{x1D735}_{\boldsymbol{v}}\unicode[STIX]{x1D735}_{\boldsymbol{v}}(\ldots )$
:
$$\begin{eqnarray}\displaystyle & \displaystyle wa_{w}^{(c)}+\frac{1}{2}(d_{ww}+d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}})=\mathop{\sum }_{\unicode[STIX]{x1D6FD}}L^{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}}\left[\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}}-\frac{1}{2}(\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FD}})_{uu}^{\prime \prime }+\left(1+\frac{m_{\unicode[STIX]{x1D6FC}}}{m_{\unicode[STIX]{x1D6FD}}}\right)w(\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}})_{w}^{\prime }\right] & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle ua_{u}^{(c)}+wa_{w}^{(c)}+\frac{1}{2}\text{tr}(\{)=\mathop{\sum }_{\unicode[STIX]{x1D6FD}}L^{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}}\left[\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FD}}+\left(1+\frac{m_{\unicode[STIX]{x1D6FC}}}{m_{\unicode[STIX]{x1D6FD}}}\right)(u(\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}})_{u}^{\prime }+w(\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}})_{w}^{\prime })\right]. & \displaystyle\end{eqnarray}$$
The
$\text{d}W$
-terms in (3.9), (3.10) are written in terms of the Rosenbluth potentials with substitution of the results (3.5), (3.6).
4 Velocity diffusion in isotropic plasma
All formulas are greatly simplified if the plasma is isotropic, that is, if the distribution functions
$f_{\unicode[STIX]{x1D6FD}}$
and the Rosenbluth potentials defining the collision operator (2.7) depend only on the velocity
$v=(u^{2}+w^{2})^{1/2}$
. In this case, the diffusion tensor has the bi-vector form:
The factoring of the diffusion bi-vector has the obvious solution:
Equation (4.1) defines in the
$\boldsymbol{e}$
-basis the block-diagonal matrix containing the non-zero blocks
$$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}d_{uu} & d_{uw}\\ d_{wu} & d_{ww}\end{array}\right)=\frac{1}{v^{2}}\left(\begin{array}{@{}cc@{}}w^{2}D_{\bot }+u^{2}D_{\Vert } & uw(D_{\Vert }-D_{\bot })\\ uw(D_{\Vert }-D_{\bot }) & u^{2}D_{\bot }+w^{2}D_{\Vert }\end{array}\right),\end{eqnarray}$$
with
$d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}=D_{\bot }$
. The block of non-zero elements of the
$b$
-matrix in (4.2) has the same form in (4.3) with a simple substitution
$\sqrt{D_{\Vert }}$
,
$\sqrt{D_{\bot }}$
instead of
$D_{\Vert }$
,
$D_{\bot }$
. Such a representation would have solved the problem of the factoring of the matrix in (4.3). But the isotropy of plasma gives the possibility for subsequent simplifications. It is easy to see that the matrices (4.2) and (4.3) can be subsequently factorized: 
where the matrices 
and 
have diagonal asymmetric blocks:
$$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}c_{uu} & c_{uw}\\ c_{wu} & c_{ww}\end{array}\right)=\frac{1}{v}\left(\begin{array}{@{}cc@{}}w\sqrt{D_{\bot }} & u\sqrt{D_{\Vert }}\\ -u\sqrt{D_{\bot }} & w\sqrt{D_{\Vert }}\end{array}\right),\quad \left(\begin{array}{@{}cc@{}}e_{uu} & e_{uw}\\ e_{wu} & e_{ww}\end{array}\right)=\frac{1}{v}\left(\begin{array}{@{}cc@{}}w & -u\\ u & w\end{array}\right)\end{eqnarray}$$
$e_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}=1$
,
$c_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}=\sqrt{D_{\bot }}$
. These relations define the differential 
by the linear forms
$\text{d}W_{\unicode[STIX]{x1D706}}\equiv (w\,\text{d}W_{u}-u\,\text{d}W_{w})/v$
,
$\text{d}W_{\unicode[STIX]{x1D700}}\equiv (u\,\text{d}W_{u}+w\,\text{d}W_{w})/v$
which turns out standard Wiener components satisfying conditions (2.3), in particular
$\text{d}W_{\unicode[STIX]{x1D706}}\,\text{d}W_{\unicode[STIX]{x1D700}}=0$
. This transformation turns out to be especially effective for the variables
$\unicode[STIX]{x1D700}=v^{2}/2$
,
$\unicode[STIX]{x1D706}=w^{2}/v^{2}$
which are chosen in many papers for studying orbital effects of collisionless motion. We can use (3.10) with substituting
$uc_{uu}+wc_{wu}=0$
,
$uc_{uw}+wc_{ww}=v\sqrt{d_{\Vert }}$
that follows from (4.4) and gives an equation containing only the Wiener differential
$\text{d}W_{\unicode[STIX]{x1D700}}$
:
Finally, the right-hand side of (4.5) can be expressed explicitly in terms of the Rosenbluth potentials:
$$\begin{eqnarray}\text{d}\left(\unicode[STIX]{x1D700}+\frac{e_{\unicode[STIX]{x1D6FC}}}{m_{\unicode[STIX]{x1D6FC}}}U\right)=\mathop{\sum }_{\unicode[STIX]{x1D6FD}}L^{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}}\left[\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}}(v)+\left(1+\frac{m_{\unicode[STIX]{x1D6FC}}}{m_{\unicode[STIX]{x1D6FD}}}\right)v\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}}^{\prime }(v)\right]\,\text{d}t+v\left[\mathop{\sum }_{\unicode[STIX]{x1D6FD}}L^{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FD}}^{\prime \prime }(v)\right]^{1/2}\,\text{d}W_{\unicode[STIX]{x1D700}}.\end{eqnarray}$$
To have the SDE for the pitch parameter
$\unicode[STIX]{x1D706}$
, we write it down as follows:
$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D706}=\frac{1}{v^{4}}(v^{2}\,\text{d}w^{2}-w^{2}\,\text{d}v^{2})-\frac{\text{d}v^{2}}{v^{6}}(v^{2}\,\text{d}w^{2}-w^{2}\,\text{d}v^{2}).\end{eqnarray}$$
This ordinary expansion is equivalent to the Ito formula if, using rules (2.4), we write the bi-linear term in the form
$D_{\unicode[STIX]{x1D700}\unicode[STIX]{x1D706}}\,\text{d}t/\unicode[STIX]{x1D700}=\tilde{\text{d}}\unicode[STIX]{x1D706}\,\tilde{\text{d}}\unicode[STIX]{x1D700}/\unicode[STIX]{x1D700}$
, where
$\tilde{\text{d}}$
denotes the Wiener part of the differential. Herein
$\tilde{\text{d}}\unicode[STIX]{x1D706}$
is defined by the linear part of (4.7):
$$\begin{eqnarray}\tilde{\text{d}}\unicode[STIX]{x1D706}=2\frac{uw}{v^{3}}[(uc_{wut}-wc_{uut})\,\text{d}W_{\unicode[STIX]{x1D706}}+(wc_{uw}-uc_{ww})\,\text{d}W_{\unicode[STIX]{x1D700}}].\end{eqnarray}$$
Substitution of coefficients (4.4) gives as a result
$uc_{wu}-wc_{uu}=-v\sqrt{D_{\bot }}$
and
$wc_{uw}-uc_{ww}=0$
. We avoided direct differentiation of
$\unicode[STIX]{x1D706}$
by Ito as it is obvious that the last term in (4.6) is zero due to the condition
$\tilde{\text{d}}\unicode[STIX]{x1D706}\,\tilde{\text{d}}\unicode[STIX]{x1D700}\sim \,\text{d}W_{\unicode[STIX]{x1D706}}\,\text{d}W_{\unicode[STIX]{x1D700}}=0$
. Therefore, the mean part of the differential,
$\bar{\text{d}}\unicode[STIX]{x1D706}$
is determined by the linear ordinary form:
$$\begin{eqnarray}\bar{\text{d}}\unicode[STIX]{x1D706}=\frac{1}{v^{2}}\left[-\unicode[STIX]{x1D706}u\left(\frac{2e}{m}\boldsymbol{h}\boldsymbol{\cdot }\boldsymbol{E}+v^{2}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{h}\right)+\frac{u^{2}}{v^{2}}(d_{ww}+d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}})-\unicode[STIX]{x1D706}d_{uu}\right]\,\text{d}t.\end{eqnarray}$$
After all substitutions, the total result takes the form:
$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}\left(\frac{\text{d}B}{B}+\frac{2e}{mv^{2}}\,\text{d}U\right)+\frac{2}{v^{4}}\left[\left(u^{2}-\frac{w^{2}}{2}\right)D_{\bot }\,\text{d}t-uvw\sqrt{D_{\bot }}\,\text{d}W_{\unicode[STIX]{x1D706}}\right].\end{eqnarray}$$
If the variations of the electrostatic potential on the magnetic surface are neglected (usually this field is weak (Galeev & Sagdeev Reference Galeev, Sagdeev and Leontovich1979)), the value
$B_{0}/B$
becomes an integrating factor for the stochastic differential
$\text{d}\unicode[STIX]{x1D706}$
. Then (4.10) can be rewritten in the most convenient form for the orbital integral
$\unicode[STIX]{x1D6EC}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D700}=\unicode[STIX]{x1D706}B_{0}/B$
:
$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D6EC}=\frac{B_{0}}{v^{2}B}\{(2-3\unicode[STIX]{x1D706})D_{\bot }\,\text{d}t+2v[\unicode[STIX]{x1D706}(1-\unicode[STIX]{x1D706})D_{\bot }]^{1/2}\,\text{d}W_{\unicode[STIX]{x1D706}}\}.\end{eqnarray}$$
Sometimes the ratio
$\unicode[STIX]{x1D709}=u/v=\pm \sqrt{1-\unicode[STIX]{x1D706}}$
is used as a pitch parameter. The equation for
$\unicode[STIX]{x1D709}$
may be found by use of the Ito formula
$\text{d}\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}^{\prime }\,\text{d}\unicode[STIX]{x1D706}+D_{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D706}}^{\prime \prime }/2$
, with
$D_{\unicode[STIX]{x1D706}}=\{2v^{-1}[\unicode[STIX]{x1D706}(1-\unicode[STIX]{x1D706})D_{\bot }]^{1/2}\}$
, and reads as follows:
$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D709}=-\frac{\unicode[STIX]{x1D706}}{2\unicode[STIX]{x1D709}}\unicode[STIX]{x1D706}\frac{\text{d}B}{B}-\frac{D_{\bot }}{v^{2}}\unicode[STIX]{x1D709}\,\text{d}t-\sqrt{\frac{D_{\bot }}{v^{2}}\unicode[STIX]{x1D706}}\,\text{d}W_{t}.\end{eqnarray}$$
It should be noted that the choice of the sign in random parts of stochastic differentials containing independent Wiener sources is arbitrary. In the last equation, the minus is adopted, since the expression for the differential
$\text{d}\unicode[STIX]{x1D706}$
contains the same Wiener noise
$\text{d}W_{\unicode[STIX]{x1D706}}$
with a sign plus. In this case the differentiation of the equality
$\unicode[STIX]{x1D706}+\unicode[STIX]{x1D709}^{2}=1$
by Ito’s rule gives zero as verification:
$$\begin{eqnarray}\text{d}(\unicode[STIX]{x1D706}+\unicode[STIX]{x1D709}^{2})=\text{d}\unicode[STIX]{x1D706}+2\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D709}+\text{d}t\unicode[STIX]{x1D706}D_{\bot }/v^{2}=0.\end{eqnarray}$$
The use of quantity
$\unicode[STIX]{x1D6EC}$
is preferable because it is the orbital invariant.
5 Spatial diffusion
The above analysis matches to the limit
$B\rightarrow \infty$
because formulation of the collisional diffusion in velocity space needs only the vector
$\boldsymbol{h}(\boldsymbol{r})$
, which determines the local direction of the magnetic field. These SDEs correspond to the Monte Carlo operators obtained in Xu & Rosenbluth (Reference Xu and Rosenbluth1991), Eriksson & Helander (Reference Eriksson and Helander1994), Tessarotto et al. (Reference Tessarotto, White and Zheng1994), Boozer (Reference Boozer2002) and others in various ways. Below the supplementary equations of spatial diffusion, which approve and refine the results (Hirvijoki et al.
Reference Hirvijoki, Brizard, Snicker and Kurki-Suoni2013), are derived by stochastic transformations of (2.10).
In this zero-order drift approximation, the spatial motion is a one-dimensional displacement along the magnetic field line
$\text{d}\boldsymbol{r}=u\,\text{d}t\boldsymbol{h}$
. The first-order drift theory of the spatial diffusion is based on the transition to the coordinates of the guiding centre of the Larmor circle of a charged particle:
$\boldsymbol{R}=\boldsymbol{r}-\unicode[STIX]{x1D746}$
, where
$\boldsymbol{g}(\boldsymbol{r})=\boldsymbol{h}/\unicode[STIX]{x1D6FA}$
. Since
$\unicode[STIX]{x1D746}$
depends linearly on the velocity, and (2.10) for
$\text{d}\boldsymbol{r}$
do not contain the Wiener terms, the differentiation of the function
$\boldsymbol{R}(\boldsymbol{r},\boldsymbol{v})$
according to Ito’s rules looks like in the ordinary analysis:
The matrix representation of the tensor 
in the basis
$(\boldsymbol{h},\boldsymbol{e}_{w},\boldsymbol{e}_{\unicode[STIX]{x1D70C}})$
in (4.11) gives:
We highlight once more that the Wiener components
$\text{d}W_{u}$
,
$\text{d}W_{w}$
,
$\text{d}V_{\unicode[STIX]{x1D70C}}$
arise as a result of a transition to a rotating triplet
$(\boldsymbol{h},\boldsymbol{e}_{w},\boldsymbol{e}_{\unicode[STIX]{x1D70C}})$
, but turns out to be independent of
$\unicode[STIX]{x1D6FC}$
together with the block-diagonal elements
$b_{wu},b_{ww},b_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}$
. Moreover the random walks of guiding centres have to be considered only in a fixed basis associated with the local value of the magnetic field. Thus, the random infinitesimal displacement (5.2) is determined by the Wiener differentials, the magnitude of which is modulated by the instantaneous values of the oscillating components of the vectors. The determination of the effective Wiener shift in the drift time scale makes sense only in a fixed basis associated with the local value of the magnetic field. To derive the drift approximation of the rotating spatial displacement, it is sufficient to consider the Ito stochastic integral over time within the defined above time interval
$(t,t+\unicode[STIX]{x0394}t)$
. The linear form (5.2) has being used as an integrand of the stochastic integral
$\unicode[STIX]{x0394}R=\int _{t}^{t+\unicode[STIX]{x0394}t}\tilde{\text{d}}\,\boldsymbol{R}(\unicode[STIX]{x1D70F})$
, determines a Gaussian random vector with zero mean components. To completely describe the Gaussian process in the drift time scale, it suffices to consider the mean expectation of a bi-linear form
$\mathit{E}\{\unicode[STIX]{x0394}\boldsymbol{R}\unicode[STIX]{x0394}\boldsymbol{R}\}$
. In view of the independence of the elementary increments of all Wiener components, the mean value of the double stochastic integral can be easily reduced to a single definite integral according to the following rule (Gikhman & Skorohod Reference Gikhman and Skorohod1972). So, it follows:
$$\begin{eqnarray}\mathit{E}\{\unicode[STIX]{x0394}\boldsymbol{R}\unicode[STIX]{x0394}\boldsymbol{R}\}=\mathit{E}\left\{\int _{t}^{t+\unicode[STIX]{x0394}t}\int _{t}^{t+\unicode[STIX]{x0394}t}\,\text{d}\boldsymbol{R}(\unicode[STIX]{x1D70F}_{1})\,\text{d}\boldsymbol{R}(\unicode[STIX]{x1D70F}_{2})\right\}=\frac{1}{\unicode[STIX]{x1D6FA}^{2}}\int _{t}^{\text{d}t}\text{d}\unicode[STIX]{x1D70F}[\boldsymbol{e}_{\unicode[STIX]{x1D70C}}\boldsymbol{e}_{\unicode[STIX]{x1D70C}}(b_{wu}^{2}+b_{ww}^{2})+\boldsymbol{e}_{w}\boldsymbol{e}_{w}b_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}^{2}].\end{eqnarray}$$
The integrand (5.3) contains the time dependence of the bi-linear forms of the oscillating vectors
$\boldsymbol{e}_{w}$
,
$\boldsymbol{e}_{\unicode[STIX]{x1D70C}}$
. If
$\unicode[STIX]{x1D6FA}\unicode[STIX]{x0394}t\gg 1$
, the integral (5.3) is asymptotically determined by the value of the bi-linear forms averaged over the angle
$\unicode[STIX]{x1D6FC}$
:
Finally, the integral (5.3) turns out to be proportional to
$\unicode[STIX]{x0394}t$
and does not depend on the gyro-oscillations:
$$\begin{eqnarray}\mathit{E}\{\unicode[STIX]{x0394}\boldsymbol{R}\unicode[STIX]{x0394}\boldsymbol{R}\}=\frac{\unicode[STIX]{x0394}t}{2\unicode[STIX]{x1D6FA}^{2}}(b_{wu}^{2}+b_{ww}^{2}+b_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}^{2})(\boldsymbol{n}_{1}\boldsymbol{n}_{1}+\boldsymbol{n}_{2}\boldsymbol{n}_{2}).\end{eqnarray}$$
Bearing in mind the relation (2.11) between the coefficients
$b$
and
$d$
, we can express the result (5.5) in terms of the diffusion coefficients in a gyro-tropic plasma:
$b_{wu}^{2}+b_{ww}^{2}+b_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}^{2}=d_{ww}+d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}}$
. In the drift time scale, in which the dynamics of the coefficients
$d$
and
$b$
is realized, the quantity (5.5) can be assumed to be infinitesimal, and (5.5) is satisfied with probability 1 without involving the averaging procedure
$\mathit{E}\{\ldots \}$
. In the drift time scale the expression (5.5) can be rewritten as follows:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\text{d}}\boldsymbol{R}\,\tilde{\text{d}}\boldsymbol{R}=(\boldsymbol{n}_{1}\boldsymbol{n}_{1}+\boldsymbol{n}_{2}\boldsymbol{n}_{2})D_{c}\,\text{d}t,\\ D_{c}={\displaystyle \frac{1}{2\unicode[STIX]{x1D6FA}^{2}}}(d_{ww}+d_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}})={\displaystyle \frac{1}{\unicode[STIX]{x1D6FA}^{2}}}\mathop{\sum }_{\unicode[STIX]{x1D6FD}}L^{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}}\left(\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}}-{\displaystyle \frac{1}{2}}(\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FD}})_{uu}^{\prime \prime }\right).\end{array}\right\}\end{eqnarray}$$
Equation (5.6) is given for the gyro-tropic case and can be also put as
In the isotropic plasma it takes the form
$D_{c}=[(1-\unicode[STIX]{x1D706}/2)D_{\bot }+(\unicode[STIX]{x1D706}/2)D_{\Vert }]/\unicode[STIX]{x1D6FA}^{2}$
. This result corresponds to the representation of the stochastic process
$\boldsymbol{R}(t)$
as diffusive two-dimensional Brownian walks across the magnetic field with the classical coefficient
$D_{c}$
:
$$\begin{eqnarray}\tilde{\text{d}}\boldsymbol{R}=\sqrt{D_{c}}(\boldsymbol{n}_{1}\,\text{d}\bar{W}_{1}+\boldsymbol{n}_{2}\,\text{d}\bar{W}_{2}).\end{eqnarray}$$
New independent Wiener components
$\text{d}\bar{W}_{1}$
,
$\text{d}\bar{W}_{2}$
are asymptotically introduced on a drift time scale, in contrast to the initial
$\text{d}W_{1}$
,
$\text{d}W_{2}$
and
$\text{d}W_{u}$
,
$\text{d}W_{w}$
,
$\text{d}W_{\unicode[STIX]{x1D70C}}$
, defined locally in time. Components
$\text{d}\bar{W}_{i}$
and
$\text{d}W_{i}$
do not correlate on the drift time scale and obey the conditions (2.3). This follows from the calculation of the mean value of the double integrals carried out in the same way as (5.3), for example:
$$\begin{eqnarray}\mathit{E}\{\unicode[STIX]{x0394}\boldsymbol{R}\unicode[STIX]{x0394}W_{w}\}=\mathit{E}\left\{\int _{t}^{t+\unicode[STIX]{x0394}t}\int _{t}^{t+\unicode[STIX]{x0394}t}\,\text{d}\boldsymbol{R}(\unicode[STIX]{x1D70F}_{1})\,\text{d}W_{w}(\unicode[STIX]{x1D70F}_{2})\right\}=-\frac{1}{\unicode[STIX]{x1D6FA}}\int _{t}^{t+\unicode[STIX]{x0394}t}\,\text{d}\unicode[STIX]{x1D70F}\boldsymbol{e}_{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x1D70F})b_{ww}=0.\end{eqnarray}$$
The derivation of the main drift differential
$\bar{\text{d}}\boldsymbol{R}$
is realized by integrating
$\text{d}t$
-terms in (5.1) over time, which is equivalent to averaging over the gyro-oscillations. As a result, we have the velocity components of the guiding centres motion known in drift theory:
$$\begin{eqnarray}\bar{\text{d}}\boldsymbol{R}=\left(u\langle \boldsymbol{h}\rangle +\boldsymbol{V}-\frac{1}{\unicode[STIX]{x1D6FA}}\langle \boldsymbol{h}\times \boldsymbol{a}^{(c)}\rangle \right)\,\text{d}t,\end{eqnarray}$$
where
$\boldsymbol{V}=\boldsymbol{V}_{\bot }+V_{\Vert }\boldsymbol{h}$
,
$\boldsymbol{V}_{\bot }=\boldsymbol{V}_{E}+\boldsymbol{V}_{B}$
,
$V_{\Vert }=w^{2}(\boldsymbol{h},\text{rot}\boldsymbol{h})/2\unicode[STIX]{x1D6FA}$
,
$$\begin{eqnarray}\boldsymbol{V}_{E}=\frac{e}{m\unicode[STIX]{x1D6FA}}\boldsymbol{E}\times \boldsymbol{h},\quad \boldsymbol{V}_{B}=\langle \boldsymbol{v}\times (\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{g}\rangle _{\bot }=\boldsymbol{h}\times \left\{\frac{w^{2}}{2\unicode[STIX]{x1D6FA}^{2}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FA}+\frac{u^{2}}{\unicode[STIX]{x1D6FA}}(\boldsymbol{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{h}\right\}.\end{eqnarray}$$
All the parameters in above formulas are defined as functions of the particle coordinates
$\boldsymbol{r}=\boldsymbol{R}+\unicode[STIX]{x1D746}$
but to close the drift description it is necessary to have the argument expressed at the fixed guiding centre coordinates
$\boldsymbol{R}$
. In the final formulas,
$\boldsymbol{R}$
should be substituted for
$\boldsymbol{r}$
since
$\langle \boldsymbol{h}(\boldsymbol{r})\rangle =\boldsymbol{h}(\boldsymbol{R})$
up to second-order terms with respect to
$\unicode[STIX]{x1D746}$
, so the account of difference
$\boldsymbol{r}-\boldsymbol{R}$
in (5.7), (5.9) would lead to an excess of the basic first-order approximation for the space drift motion. The contribution of the frictional force in the expression (5.9) is zero at this order, since the transverse friction term
$\boldsymbol{a}^{(c)}$
in (5.9) is proportional to the vector
$\boldsymbol{e}_{w}$
. Finally, considering (5.7), the equation for charged particle spatial motion is represented by a three-dimensional stochastic equation in vector form that includes a two-dimensional isotropic Wiener scattering process with the coefficient
$\sqrt{D_{c}}$
$$\begin{eqnarray}\text{d}\boldsymbol{R}=(u\boldsymbol{h}+\boldsymbol{V})\,\text{d}t+\sqrt{D_{c}}(\boldsymbol{n}_{1}\,\text{d}\bar{W}_{1}+\boldsymbol{n}_{2}\,\text{d}\bar{W}_{2}).\end{eqnarray}$$
We point out that general vector equation (5.11), written in the proper coordinate representation, is applicable for any equilibria when the drift approximation is proved. Below we consider the application of the algorithm to axisymmetric equilibria as a practical important example. In the case of an axially symmetric magnetic configuration of the tokomak-type devices, the magnetic field decomposes into poloidal and toroidal components,
$\boldsymbol{B}=\boldsymbol{B}_{p}+\boldsymbol{B}_{t}$
, represented in cylindrical coordinates
$\boldsymbol{R}=(R,\unicode[STIX]{x1D717},Z)$
,
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}=\boldsymbol{e}_{\unicode[STIX]{x1D717}}/R$
:
$$\begin{eqnarray}\boldsymbol{B}_{p}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}(R,Z),\quad \boldsymbol{B}_{t}=I(R,Z)\unicode[STIX]{x1D735}\unicode[STIX]{x1D717},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6F9}$
and
$I$
are poloidal magnetic and current fluxes. It is natural to introduce the positive triplet
$(\boldsymbol{h},\boldsymbol{n}_{1},\boldsymbol{n}_{2})$
by formulas:
$$\begin{eqnarray}\boldsymbol{n}_{1}=\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}|},\quad \boldsymbol{n}_{2}=\boldsymbol{h}\times \boldsymbol{n}_{1}=\frac{B_{p}}{B}\boldsymbol{e}_{\unicode[STIX]{x1D717}}-\frac{h_{t}}{B_{p}}\boldsymbol{B}_{p}.\end{eqnarray}$$
The SDE (5.11) together with equations for
$u,w$
(3.8) makes it possible to get a picture of the poloidal motion of charges in a magnetic field in the poloidal plane
$(R,Z)$
that would demonstrate the nature of charged particles neoclassical diffusion with trapping and collision scattering effects.
For this purpose, we should educe SDE (5.11) to stochastic differentials in the cylindrical coordinates
$R(\boldsymbol{R})$
,
$Z(\boldsymbol{R})$
using Ito formulas. The nonlinearity of this transformation causes the contribution of the second derivatives
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}R=\boldsymbol{e}_{\unicode[STIX]{x1D717}}\boldsymbol{e}_{\unicode[STIX]{x1D717}}/R$
in the equation for
$\text{d}R$
. This transformation results in the additional
$\text{d}t$
-term proportional to
$D_{c}$
. However, the account of the second-order quantity
$D_{c}$
would exceed accuracy of drift approximation. Therefore, the drift differentials
$\bar{\text{d}}R$
,
$\bar{\text{d}}Z$
can be obtained by the rules of ordinary analysis:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{d}R=(uh_{R}+V_{R})\,\text{d}t+\sqrt{D_{c}}(n_{1R}\,\text{d}\bar{W}_{1}+h_{t}n_{1Z}\,\text{d}\bar{W}_{2}),\\ \text{d}Z=(uh_{Z}+V_{Z})\,\text{d}t+\sqrt{D_{c}}(n_{1Z}\,\text{d}\bar{W}_{1}-h_{t}n_{1R}\,\text{d}\bar{W}_{2}),\end{array}\right\}\end{eqnarray}$$
where
$h_{t}=I/BR$
. Equation (5.14) are applicable for both gyro-tropic and isotropic plasmas but require the poloidal distributions of
$\unicode[STIX]{x1D6F9}(R,Z)$
,
$I(R,Z)$
. It is of interest to derive SDEs for particle motion in the axisymmetric plasma in terms of pseudo-cylindrical coordinates
$\boldsymbol{R}=(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D717})$
defined by transformation
$R=R_{0}-r\cos \unicode[STIX]{x1D703}$
,
$Z=r\sin \unicode[STIX]{x1D703}$
, assuming the dependence of
$\unicode[STIX]{x1D6F9}$
only on
$r$
. In this case we have
$\boldsymbol{h}=[\unicode[STIX]{x1D6F9}^{\prime }\boldsymbol{e}_{\unicode[STIX]{x1D703}}+I\boldsymbol{e}_{\unicode[STIX]{x1D717}}]/BR$
,
$\boldsymbol{n}_{1}=\boldsymbol{e}_{r}$
,
$\boldsymbol{n}_{2}=[I\boldsymbol{e}_{\unicode[STIX]{x1D703}}-\unicode[STIX]{x1D6F9}^{\prime }\boldsymbol{e}_{\unicode[STIX]{x1D717}}]/BR$
, and the poloidal model (5.14) reduces to the following equations:
$$\begin{eqnarray}\text{d}r=V_{\bot r}\,\text{d}t+\sqrt{D_{c}}\,\text{d}\bar{W}_{1},\quad r\,\text{d}\unicode[STIX]{x1D703}=(uh_{\unicode[STIX]{x1D703}}+V_{\unicode[STIX]{x1D703}})\,\text{d}t+h_{t}\sqrt{D_{c}}\,\text{d}\bar{W}_{2}.\end{eqnarray}$$
Equation (5.15) must be solved together with (3.8) or (3.9), (3.10), determining the velocity parameters
$u=\unicode[STIX]{x1D70E}v\sqrt{1-\unicode[STIX]{x1D706}}$
,
$w=v\sqrt{\unicode[STIX]{x1D706}}$
,
$v=\sqrt{2\unicode[STIX]{x1D700}}$
,
$\unicode[STIX]{x1D70E}=\pm 1$
. Such a simplified model is sufficient to elucidate the features of the spatial trajectories of the bulk plasma particles as well as the additional spices (injected and ICRF (ion cyclotron resonant frequency) heated ions or fusion products).
6 Conclusions and discussion
The above analysis provide the method for the simulation of drift stochastic motion of charged particles in plasmas via the set of four SDEs, i.e. two equations for the orbit invariant scattering in the velocity space (3.9), (3.10) or (4.6), (4.11), and two SDEs (5.14) and/or (5.15) for spatial guiding centre convection and cross-field diffusion. The simplest formulation is possible in the case when the Rosenbluth potentials are determined in isotropic plasma using the axisymmetric pseudo-cylindrical representation of SDEs presented by the set of (4.6), (4.11), (5.15), that describe the 4-D diffusion process in the phase space
$(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D700},\unicode[STIX]{x1D706})$
. All coefficients in these equations can be explicitly calculated in the approximation of Maxwellian background plasma. Each of these four equations includes its own Wiener noise independently, that also facilitates the numerical solution of the SDEs (Kloeden & Platen Reference Kloeden and Platen1995).
It should be noted that the theory developed above determines the presence of Wiener terms in the equations for spatial diffusion (5.11), (5.14), (5.15), that are usually absent in all up-to-date modelling where the authors restrict themselves to the definition of the drift velocity scattering due to collisions. The Wiener differentials, being formally defined in the same first order of the drift theory in addition to the
$\text{d}t$
-terms, dominate in the expression for the stochastic differentials, since their order is estimated by the quantity
$\sqrt{\text{d}t}$
. It is obvious that a correct description of the random trajectories by the Monte Carlo method is possible with strict account of all components of the stochastic differentials.
The SDEs obtained above, as well as a result of Hirvijoki et al. (Reference Hirvijoki, Brizard, Snicker and Kurki-Suoni2013); describe the motion of particles in the main drift time scale, which makes it possible to exclude only the particle gyration, while the orbital motion of the guiding centres and their scattering by collisions are taken into account together. Thus, the spatial neoclassical diffusion is introduced by the classical diffusion process. Stochastic solutions in a long-time scale must contain scattering results of passing or trapped trajectories or their fragments which occur within the total time of particle motion. These SDEs equations are free from the limitations of the axial symmetry magnetic configuration, which is the condition for applicability of the Monte Carlo operators for modelling orbit-averaged particle diffusion according to the models (Xu & Rosenbluth Reference Xu and Rosenbluth1991; Eriksson & Helander Reference Eriksson and Helander1994; Tessarotto et al. Reference Tessarotto, White and Zheng1994; Boozer Reference Boozer2002). The advantage of the presented SDEs description of all the processes affecting the particle motion is the correct account of the Coulomb collisions and particle trajectories, with possible relative increase of the computational time. Thus one can expect that a random walk along the orbit fragments would reflect the diffusion process more correctly.
The proposed poloidal diffusion model serves to visualize the projection of random trajectories motion on the poloidal plane of the plasma column cross-section. This approach also makes it possible to determine the azimuthal distributions of particle density and energy loading to the plasma facing elements, depending on the initial distribution of the ensemble of test particles. Such results are hardly may be achieved using the bounce-averaged Monte Carlo operators utilized in the codes (Xu & Rosenbluth Reference Xu and Rosenbluth1991; Eriksson & Helander Reference Eriksson and Helander1994; Tessarotto et al.
Reference Tessarotto, White and Zheng1994; Boozer Reference Boozer2002). The poloidal diffusion model proposed in (5.10), (5.13), (5.14) being especially attractive for applications enables the enhanced diffusion coefficients estimated locally for various azimuthal and radial coordinates by the mean values
$D_{\text{neo}}=\mathit{E}\{(\unicode[STIX]{x0394}r)^{2}\}/\unicode[STIX]{x0394}t$
for the sample path ensembles within the drift time intervals
$\unicode[STIX]{x0394}t$
. The numerical model of the stochastic diffusion induced by classical processes presented here are strictly consistent with the kinetic theory of Coulomb collisions and can serve as a standard gauge for the analytical models of the neoclassical diffusion theory and for the verification of the appropriate numerical approaches.
