\documentclass{article}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\beq}{\begin{eqnarray}}
\newcommand{\eeq}{\end{eqnarray}}
\newcommand{\ds}{\displaystyle}
\newcommand{\ts}{\textstyle}
\newcommand{\li}{\limits}
\renewcommand{\abstractname}{}

\topmargin=-10mm \oddsidemargin=-10mm \evensidemargin=-10mm
\textwidth=170mm \textheight=240mm

\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{graphics}
\usepackage [koi8-r] {inputenc}
\usepackage [T2A] {fontenc}
\usepackage [english,russian] {babel}

\usepackage {epsfig}
\usepackage {floatflt}

\usepackage {psfrag}

\usepackage [math]{pscyr}

%% use PSCyr fonts:
\newcommand\fixfonts{%
  \renewcommand {\rmdefault}{faq}%
  \renewcommand {\sfdefault}{ftx}%
  \renewcommand {\ttdefault}{fco}%
}
\fixfonts


\psfrag {phase}[][b][0.8]{phase}
\psfrag {U}[][b][0.8]{\(u\)}
\psfrag {sqrU}[][b][0.8]{\(u^2\)}
\begin{document}

\fixfonts

\title{Dynamic theory of relativistic electrons stochastic heating
by whistler mode waves with application to the Earth magnetosphere
 }

\date{}
\maketitle
\section{Introduction}

\indent It is suggested that whistler mode waves are responsible for
electron precipitation caused by the interaction of these waves
and trapped radiation belt electrons. Research has focused
mainly on the interaction of waves whose Doppler-shifted frequency
matches that of the electron cyclotron frequency (Walker, 1993).

\noindent It is known (e.g. Baker et al. 1986; 1994) that the flux of
outer radiation zone relativistic (\(E \sim 1\)~MeV) electrons
exhibits considerable variability during geomagnetic storms. Horne and
Thorne (1998) have identified potential wave modes that are capable of
resonating with electrons over the important energy range from 100~keV
to a few Mev in different regions of Earth's magnetosphere, 
principally whistler modes. The basic concept of energy diffusion of
relativistic electrons resulting from resonant interaction with
whistlers in the magnetosphere has been discussed, for example, by Li
et al. (1997) and Summers et al. (1998).

While these theories have been successful in explaining some of the
events, there is experimental evidence that other processes are also
important (Goldberg et al., 1987). Moreover it is known that one given
mechanism of resonance is in fact, unlikely, because of chaotic, particle
motion caused by nonlinear wave-particle interaction as
well as by the group velocity dispersion (Lichtenberg and Lieberman,
1983; Zaslavsky, 1985).  Attracting interest for many years,
stochastic dynamics of
charged particles in the field of a wave packet is one of the fundamental
problem in the theory of plasma physics; e.g., the stochastic dynamics 
stochastic dynamics of relativistic
electrons in the time-like wave packet has been discussed by Chernikov
et al. (1989). A theory of dynamic chaos has also been employed to
deal with wave-particle interactions of the high-frequency heating of
electrons in the space-like wave packet (Klimov and Tel'nikhin,
1995). In these papers respectively, the relativistic generalization of
the standard maps introduced by Chirikov (1979) and Chernikov and
Zaslavsky (1985) are derived. Possible applications of these results
were discussed by Karimbady and Angelopous (1989), and Krotov and
Tel'nikhin (1998).Nagornykh and Tel'nikhin (2002) have announced
some results of theoretical study concerning stochastic particle
motion in the magnetized 0.  In addition, electron precipitation 
caused by chaotic
motion due to coupling the bounce motion of trapped electrons with
a large-amplitude whistler wave has been studied by Faith et al
(1997). The interaction between two oppositely directed whistler waves
and non-relativistic electrons, and the appearance of chaos in this
system has been considered by Matsoukis et al (1998). Wykes et
al. (2001) have generalized this system to include relativistic
effects. With applications to Jupiter in a given particle 
energy approximation, the degree of pitch
angle scattering as a function of electron energy and initial pitch
angle has been calculated.

The research presented here, systematically studies stochastic motions
of high-energy electrons due to a spectrum of whistler mode waves with
applications to the Earth magnetosphere. The paper is organized as
follows:
\begin{enumerate}
\item We derive the canonical equation of motion in terms of the 
action-angle variables.

\item In section two, we discuss the nonlinear resonance 
wave particle interaction.

\item In section three, it is shown that the stochastic motion of 
relativistic electrons can be described by a closed set of nonlinear 
difference equations.  The solutions of these equations are obtained
both analytically and numerically.

\item We derive equations in section four that describe the dynamics of
resonant non-relativistic electrons in the time-like wave packet.  
We show also that solutions of these equations allow chaotic motion.

\item The Fokker-Planck-Kolmogorov equation is derived in section five 
and the affects
associated with particle stochastic heating are studied.  

\item Application
of our results to high-energy electrons observed in Earth's election
radiation belts is described in section six.

\item In section seven, we discuss the
conclusions of our studies.
\end{enumerate}

\section{Basic equations}


Let us consider a relativistic particle of charge \(|e|\) and mass
 \(m\) in the wave packet of extraordinary electromagnetic waves
 propagating along an external uniform magnetic field of strength
 \(B\). The Hamiltonian corresponding to the problem is

\be 
H(\bar{r}, \bar{p};t)= \sqrt{m^2+ (\bar{p}+\bar{A})^2},
\label{1} 
\ee 
and the canonical equations of motion are 
\be
\dot{\bar{p}} = [\bar{p},H], \ \ \ \ \dot{\bar{r}} = [\bar{r},H],
\label{2} 
\ee 
where \(\bar{p}\) is the particle momentum, \(\bar{r}\) is the
position vector, \(\bar{A}\) the vector potential divided into two
parts, \(\bar{A}=\bar{A}^w +\bar{A}^{ext}\), superscripts \(w \) and
\(ext\) denote both the wave and external fields, and [\ \ ,\ \ ]
stand for the Poisson brackets.\\
We have employed here and throughout this paper, the frame of reference in which the
speed of light \(c=1\) and charge \(|e|=1.\)

\noindent We denote by \(R\) the set of all real numbers. Then
\(\bar{p}\in R^3, \bar{r}\in R^3,\) and the smooth manifold \(M =
R^6\) will be a canonical space of this dynamic system, and \(R^6=R^3\times
R^3 \) is a direct product space.

In order to write down equations of the particle motion one must
specify a coordinate system. We have chosen a Cartesian spatial
coordinates system whose \(z\) axis is directed along the external
magnetic field, the plane perpendicular to this direction is spanned
by the orthogonal coordinates \(x\) and \(y\).

\noindent Making use of the connecting relations, 
\beq
\bar{B}=[\bar{\nabla},\bar{A}], \ \ \ \bar{E} =
-\partial{\bar{A}}/\partial{t}, \label{3}\eeq we have in the
coordinate representation \beq \bar{r} =(x,y,z), \ \ \
\bar{B}^{ext}=(0,0,B), \label{4} 
\eeq 
\beq \bar{A}^w=\left(\sum_k
A_k \sin\varphi,\ \ \sum_kA_k\cos\varphi , 0\right),\ \  \varphi =
zk -t\omega_k, \label{5} 
\eeq 
\beq \bar{A}^{ext}
=\frac{1}{2}(-By,Bx,0). \label{6}
\eeq 
Here the expression for \(\bar{A}^{ext}\) is written in the axial
gauge, \(A_k\) is the amplitude of mode in the wave packet, \(k\) is
the wave number, and \(\omega_k\) is the dispersion equation.

The dispersion relations for electron branch of the whistler mode
waves in the cold magnetoplasma is written as 
\beq
\frac{k^2}{\omega^2}= 1+\frac{\omega_p^2}{\omega(\omega_B
-\omega)}, \label{7}
\eeq
where \(\omega_B\) and \(\omega_p\) are the gyrofrequency and electron
plasma frequency, respectively.

Equation (\ref{7}) is a long-wavelength approximation for the low
density plasma, viz., the condition \((\omega_B \omega/\omega_p^2)
\ll 1\) is valid,  reduces to 
\beq  v_{ph}^2=
\frac{\omega(\omega_B-\omega)}{\omega_p^2}, \ \ v_{ph}^2\ll 1.
\label{9} 
\eeq

\noindent We now take into account axial symmetry of the
non-perturbative problem, and  introduce the new variables, an
action\((I)\), and an angle \((\theta )\), by a canonical
transformation \((x, p_x; y, p_y) \rightarrow (\theta, I):\) 
\beq
x=r\cos\theta, && p_x=
-\frac{1}{2}mr\omega_B\sin\theta,\label{10} \\
y=r\sin\theta, && p_y =\frac{1}{2}mr\omega_B\cos\theta;
\nonumber \\
r = \sqrt{2m\omega_B I}/m\omega_B, && \omega_B =B/m, \label{11}
\eeq 
where \(r\) is the gyroradius.

The Hamiltonian (\ref{1}) in this representation becomes 
\beq
H(z,p_z; \theta,I;t)= H_0(p,I)+ \frac{\sqrt{2m\omega_B
I}}{H_0}\sum_kA_k\cos(zk+\theta -t\omega_k),\label{12}
\\
H_0(p,I)=\sqrt{m^2+p^2+2m\omega_B I}.
\phantom{H(z,p_z; \theta,I;t)= H_0(p,I)+aaa}
\label{13} 
\eeq

\noindent Here we have assumed that the ratio \(\mu = A/m \ \ll 1\) is
the small parameter of the problem, and one retains in (\ref{12}) only
the leading terms.

Associated with (\ref{12}) the  equations of motion are
\begin{eqnarray}
\dot{p}=[p,H]=\frac{\sqrt{2m\omega_B I}}{H_0}\sum\limits_kkA_k\sin
\psi ,
&&\dot{I}=[I,H]=\frac{\sqrt{2m\omega_B I}}{H_0}\sum\limits_kA_k\sin \psi ,  \label{14}\\
\dot{z}=[z,H] =pH_0^{-1},\phantom{\dot{p_z}=a a \dot{p}=[p,H]}
&&\dot{\theta}=[\theta,H]= \omega_BmH_0^{-1}.   \label{15}
\end{eqnarray}

\noindent We omit in (\ref{15}) the terms of the order of \(\mu^2\)
and introduce the definition for the phase 
\be \psi\stackrel{def}{=}
zk+\theta-\omega_kt. \label{16} 
\ee 
It is obvious that a Hamiltonian
flow  (\ref{14}, \ref{15})   preserves the measure \(d\mu
=dp_zdIdzd\theta\) on the smooth manifold \(M = R^3\times S\), where
\(S\) is  the circle and \(\theta (\mbox{mod}
2\pi) \in S.\) Two pairs of the canonical variables \(p, z\) and \(I, \theta\) are
symplectic coordinates on \(M\) with a little algebra:
\[[e_p, e_{\theta}]= [e_p, e_I]= [e_I,e_z]= e_z,e_\theta] = 0,\]
\[[e_p, e_z]=[e_I,e_{\theta}] = 1,\]
where \([,]\) denote a skew-product, the set \((e_j)\) is the
coordinate basis, which can be given by a differential operators
 \(\partial/\partial x^j \) (Schutz, 1982).\\ 
\noindent Again with a little algebra, we now proceed
to the physical aspect of the problem. \\

First we consider the dynamics of an autonomous system with the
Hamiltonian given by (\ref{13}). In this commutative integrable
situation \((\dot{I}=0,\ \ \dot{p}=0)\), the phase curves lie on
constant energy surfaces of \(H_0= E_p=\) const:
\[\{(z,p;\theta,I):\sqrt{m^2+p^2+2m\omega_BI}= E_p\}.\]
All solutions are regular, since trajectories in phase space
cannot intersect and each phase curve correspond in configuration
space to particle gyration in an external magnetic field with a
frequency
\(\omega_BmH_0^{-1}.\)\\
To understand the physical picture of the stochastic motion, let us
discuss the simple case of the single-wave Hamiltonian. In order to
investigate a non-autonomous system we start with a particular
Hamiltonian
\be H=H_0(p,I)+\frac{\sqrt{2m\omega_B
I}}{H_0}A_{k}\cos (zk+\theta-\omega_kt), 
\ee 
where \(A_{k}\) is the amplitude of certain mode in the wave packet,
\(k\) is an arbitrary wavenumber, \(\omega_k\) is the frequency of
this mode.  \\
Thus the associated equations of motion are
\begin{eqnarray}
\dot{p}&= &k\frac{\sqrt{2m\omega_B I}}{H_0}A_{k}\sin \psi , \ \ \
\dot{I}=\frac{\sqrt{2m\omega_B
I}}{H_0}A_{k}\sin \psi ,  \label{18}\\
\dot{z}\phantom{a}&=&pH_0^{-1}, \  \ \
 \phantom{aaa\sqrt{2m\omega_B I}}\dot{\theta}= \omega_BmH_0^{-1},
 \label{19}  \\
\dot{\psi_k}&=&kpH_0^{-1}+\omega_BmH_0^{-1} -\omega_k. \label{20}
\end{eqnarray}

\noindent It is easy to see that the system possesses the first
integral, \ \(F\), \ \  \(\dot{F} = (\dot{p}-k\dot{I})=0\),  since
\([F,H]=0\). We put const.= 0 below and represent the invariant of
motion as 
\be \alpha p - 2\omega_B I =0,  \ \ \alpha =
\frac{2\omega_B}{\omega}v_{ph}.\label{21} 
\ee

We specify the problem. We shall be considering the high-energy
particle dynamics such that the expression for particle energy,
\(E_p\),
 \be
 \varepsilon =
\varepsilon_z(1+\varepsilon_z^2/2\varepsilon_z^2
+1/2\varepsilon_z^2), \label{22} 
\ee 
\be \varepsilon \equiv E_p/m,
\ \ \   \varepsilon_z \equiv p/m, \ \ \ \varepsilon_t \equiv
\sqrt{2m\omega_BI}/m,  \label{23}
 \ee
suits the problem.\\

It is well known that the stochastic motion will be mainly
determined by the nonlinear resonance wave-particle interaction as
well as the group velocity   dispersion. Consequently, we must
consider the dynamics determined by conditions \be
\dot{\psi_k}=(kpH_0^{-1}+\omega_BmH_0^{-1} -\omega_k)_{\mbox{\it
res}}=0. \label{24}
 \ee
It is useful to discuss here the structure of phase space. As may well
be imagined, provided that (\ref{21}) exists, the number of dimensions
of the system reduces to two. Let the pair \(\theta, I\) be a
fundamental pair of variables for a given problem. As a result, we obtain a
frame rotating with angular velocity \((\omega_k -
kpH_0^{-1})_{res}\).  Relation (\ref{9}) along with (\ref{24})
determines in implicit form the family of resonance tori; the phase space
of the system has a natural structure of fibre bundle \((M,B)\), where
\(M\) is the smooth manifold,  \(M=B\times S.   \ \ \ \ M \rightarrow B\)
is the base space, \(B\subset R=\{I\},\) and \( \ M \rightarrow S,
\theta(\mbox{mod} 2\pi) \in S,\) so that each fibre is assigned the
torus from the family given by (\ref{24}).\\

The invariant (\ref{21}) along with (\ref{9},\ref{24})  allows us
to represent the interval between the adjacent tori from the
family of resonance tori located on the base   as 
\be
\left(\frac{\delta p}{m}\right)_{res}  = \mbox{min} \left\{
\frac{\omega_{dis}}{\omega_B}\varepsilon_z^2 \right\}_{res},
\label{25} 
\ee 
where \(\omega_{dis} \) is a dispersion of wave frequencies,
\be \omega_{dis} =
\mbox{min}\left\{ \delta\omega, \omega_B\varepsilon ^{-1}
\right\}, \ \ \ \delta\omega= \mbox{max}\left\{ \Delta
kpH_0^{-1},\Delta \omega \right\},    \label{26} 
\ee
\[\Delta k=L^{-1}, \ \ \ \Delta\omega=T^{-1},\]
\(T\) and \(L\) are the time -- and distance scales, respectively.\\
We now specify the wave spectrum. Let a particle speed is much
larger than the  group velocity,\(V_{gr}\) of wave packet, and the
characteristic spacescale \(L\) sufficiently greater than the
gyroradius, i.e. the conditions 
\be p/H_0 \gg V_{gr} , \ \ \
\varepsilon < N , \ \ \ N\equiv [kL] \label{27} 
\ee
are valid, and \(N\) represents the characteristic number of modes in
the wave packet, \([\ \ ]\) denotes here the integer part. It should
be noted the situation described by (\ref{9}) is just the
case. Then the relation (\ref{25}) with (\ref{21}) and (\ref{27})
yields
\be
\left(\frac{\delta p}{m}\right)_{res}  =
\frac{k}{\omega_B}N^{-1}\varepsilon_z^{2} \label{28}
 \ee
Next, we consider the behavior of this system in the  vicinity of some
resonance torus. Denote by \(\xi=(p-p_{res})/m\)  the new variable
in the proximity space. The equations (\ref{18}-\ref{20}) in terms
of the variables \((\xi , \psi)\) become 
\be 
\dot{\xi}=a\sin \psi, \
\ \ \dot{\psi}=-b\xi, \label{29} 
\ee 
where \(a\) and \(b\) are the functions of parameters evaluated on
definite resonance torus, namely
\be 
a=k(A/m)\alpha^{1/2}\varepsilon^{-1/2}, \ \ \ b=
\omega_B\varepsilon_z^{-2}. \label{30} 
\ee 
Equations (\ref{29}) possess canonical structure and are governed by
the Hamiltonian, \(h(\xi,\psi),\)
 \be
h(\xi,\psi)=a\cos\psi-\frac{1}{2}b\xi^2. \label{31} 
\ee 
When \(h=-a\) appropriate to the problem, we establish an upper bound
for the set
\be \sup \{\xi\}= \left
(2\alpha^{-1/2}\varepsilon_z^{3/2}(A/m)\right)^{1/2}.
       \label{32} 
\ee 
\\
Further, if the product of sets \(\{\xi\}\)   and \(\{(\delta
p/m)_{res}\},\) \(\{\xi\}\cap\{(\delta p/m)_{res}\}\neq 0\), where
\((\delta p/m)_{res}\)  is determined by (\ref{28}), then the
condition \(\sup \{\xi\} \geq (\delta p/m)_{res} \) infers that the
union of subsets  \(\bigcup\{\xi\}\) is a connected dense set. This
criterion may be written in an explicit form as follows
 \be
 \left
(2\alpha^{-1/2}\varepsilon_z^{3/2}(A/m)\right)^{1/2} \ge
\frac{k}{\omega_B}N^{-1}\varepsilon_z^{2} \label{33} 
\ee 
Finally, we find from (\ref{33})

\be 
\varepsilon \leq \varepsilon_b =
\left(1/2\alpha^{3/2}N^2(A/m)\right)^{2/5}. \label{34} 
\ee 
This means that in the ultra-relativistic limit, when the phase
velocities in a packet differ significantly from the speed of a
particle, the energy of the particle is bounded from above. As may be seen
in (\ref {34}), there is little dependence of \(\varepsilon_b\) on the
magnitude of the wave field, \(A\).


The expressions (\ref{28}), (\ref{34}) imply the possibility of
self-generated, chaotic behavior of the system in some region of
energetic space.  The results may be interpreted as follows. It is
easy to see that the structure of fibering mentioned above corresponds
to quantized energy space, when each fibre is assigned a stationary
state. Each stationary state has its own energy level, and all the levels
are quantized. The energy of a particle in its stationary states and
the interval between adjacent levels are given by (\ref{9}) and
(\ref{24}), and (\ref{25}), respectively. The width of each energetic
level in turn is governed by the resonance wave-particle
interaction. Thus the condition (\ref{33}) is, as a matter of fact,
some uncertainly principle and it is therefore not possible to
determine the precise energy of particles.  With regard to the particle,
only the probability of the energetic state of particle can be given
with precision.  Note also that expressions (\ref{28}) and (\ref{32})
have different dependencies on \(p,(\delta p/p) \propto p\),
\(\sup{\{\xi}\}/p\propto p^{-1/4},\) whence it follows that
limitations (\ref{34}) on stochastic heating relate directly to the
non-equidistant energy spectrum. So the results have a 
physical interpretation, 
but it is better to treat the conceptual aspects
of nonlinear motion mathematically.  We will carry on our advanced
studies of stochastic motion in the next sections.


\section{Particle dynamics in a space-like wave packet}

\indent We specify first the wave spectrum of a packet.  Let us
assume that the conditions (\ref{27}) are valid.  Thus, we will 
consider the ultra-relativistic particles motion in the space-like
wave packet in magnetized plasma (Lichtenberg and Lieberman, 1983,
Zaslavsky 1985). In this approach the equations of motion (\ref{14})
and (\ref{15}) take the open form:
\beq
\dot{p}=kA\frac{\sqrt{2m\omega_BI}}{H_0}\sin\psi\sum\limits_{n\in
Z}\delta(\zeta-n),\phantom{aaaaa} &&
\dot{I}=\frac{\sqrt{2m\omega_BI}}{H_0}A\sin\psi\sum\limits_{n\in
Z}\delta(\zeta-n), \label{35}
 \\
\dot{z}=pH_0^{-1}, \phantom{\dot{I}=sU\sin\psi\sum\limits_{n\in
z}\delta(\zeta a\,}&
&\dot{\theta}=\omega_BmH_0^{-1};\phantom{2m\omega_BI\delta (\tau-n)WWWW}
 \label{36}\\
H_0(p,I)=\sqrt{m^2+p^2+2m\omega_BI},\phantom{2m\omega_BIa a }
&&\psi=zk+\theta- \omega t.
\label{37}
 \eeq
\noindent Here \(A, \omega\) and \(k\) are the magnitude, frequency,
and wave number of the fundamental (characteristic) mode, 
\(\zeta=(z/L)\)(mod 1), \(L\) is the characteristic spacescale,
\(\delta_n\equiv \delta(\zeta-n), \ \ \delta (\cdot)\) is the Dirac
delta function, and \(Z\) denotes the set of all integers.  \\ These
equations are, in fact, a mathematical model of the construction mentioned
above. The evident symmetries of a turning, \(\psi \rightarrow \psi
\exp (i\psi_0), \psi_0\) is an arbitrary constant, and a translation,
\(\zeta \rightarrow \zeta + 1\), allow us to represent dynamics of the
system as a certain iterative process by identification of the planes
\(n\zeta \) and \((n+1)\zeta\). The group symmetry is realized first
as the motion invariant \(F\),
\be
F=p-kI=inv. \label{38} 
\ee
Using (\ref{23}), we write down expression (\ref{38})
as
\be 
\varepsilon_z-\alpha^{-1}\varepsilon_t^2=0, \ \ \alpha =
(2\omega_B /\omega)v_{ph},
\label{39}
 \ee where the constant of integration has been chosen to equal zero.  In
view of (\ref{39}) the number of dimensions reduces to two. Let the
variables \(z\) and \(p\) will be the represented pair.  \\
Now we find the explicit form of the iteration system. Denote by
\(TM\) the tangent fibering. There exits the map \(\bar{v}:
M\rightarrow TM\), such that the field of vectors has the form
\be
\bar{v}=k\frac{\sqrt{2m\omega_BI}}{H_0}A\sin
\psi\sum\delta_n\frac{\partial}{\partial p}+
pH_0^{-1}\frac{\partial}{\partial z}.  \label{40} 
\ee
 Then each
element of  \(\bar{v}\) is assigned the  1-form, so that the field
of 1-forms may be written as 
\beq
\tilde{v}=-k\frac{\sqrt{2m\omega_BI}}{H_0}A\sin\psi\sum\delta_n\tilde{d}z+pH_0^{-1}\tilde{d}p.
\label{41}
 \eeq
The notations \(\partial/\partial x^i \) and \(\tilde{d}x_i\) are
employed here for the orthonormal coordinate basis of vectors.  With the
algebra mentioned above and the dual basis of 1-forms, such that
\[\frac{\partial}{\partial x^i} \tilde{d}x_{i^{\prime}}=\delta_i^{i\prime}=
\left\{\begin{array}{ll}
1,&i=i^{\prime}\\
0,&i\ne i^{\prime}
\end{array}\right.,\]
where \((x^i)=(p,z)\). It is well known that the section of the field of
1-forms, \(\tilde{v}(\bar{v})=0\), on a submanifold in \(M\) is
equivalent to the solution of equations (\ref{35}, \ref{36}) ( Schutz,
1982). By convolution transform \(\tilde{v}(\bar{v})=0\), we deduce
from (\ref{40} -- \ref{41})
\beq 
\tilde{d}\varepsilon_z -
\varepsilon^{-1/2}\alpha^{1/2}N(A/m)\sin\psi\sum_{z\in Z}\delta_n
\tilde{d(\zeta)} =0, \label{42} 
\eeq 
and then equations (\ref{36}) and (\ref{37}) yield

\beq 
\tilde{d}\psi -(\omega -  \omega_B\varepsilon^{-1}
-k\varepsilon_z\varepsilon^{-1})L\frac{\varepsilon}{\varepsilon_z}\tilde{d}\zeta=0.
\label{43}
\eeq 
Making use of the invariant of motion (\ref{39}), and integrating one
by one resulting equations, we obtain the closed set of nonlinear
difference equations
\be
\begin{array}{ll}
& u_{n+1}=u_n+3/2Q\sin\psi_n,\\
g^n:&\\
& \psi_{n+1}=\psi_n+N\left( 1+\frac{1}{2}\alpha
u_{n+1}^{-2/3}-u_{n+1}^{-2/3}v_{ph}\sqrt{1+u_{n+1}^{4/3}+\alpha
u_{n+1}^{2/3}} \right) (\mbox{mod }2\pi), \label{44}
\end{array}
\ee 
expressed in terms of the variables \(u, \psi\), where \(u_{n+1}\) and
\(u_n\) are respectively the values of the normed momentum at times
\((n+1)\zeta\) and \(n\zeta\), and \(\psi_{n+1} -\psi_n\) is the phase
shift acquired by the particle. We have employed here the following notations:
\beq
 u= \varepsilon_z^{3/2}, \ \
Q =\sqrt{\alpha}\cdot Nw, \ \ \alpha =
2\frac{\omega_B}{\omega}v_{ph}, \ \ w=A/m. \label{45}
\eeq
The quantity \(w\) has a clear physical meaning: \(w\) is the
dimensionless representation of the ratio, at one wavelength, of the 
work of the wave field
to the particle rest energy. It is convenient to
represent an expression for \(w\) in some diverse form. Regarding the
relationship between fields \(A^w\) and \(B^w\), we right down
\be
w = \frac{1}{2}\alpha b, \ \ b = B^w / B.
\label{w}
\ee

It is readily seen, if the inequality \(\varepsilon \gg \alpha \) is
valid, the set of equations (\ref{44}) goes over into the following
equations: 
\be\begin{array}{ll}
& u_{n+1}=u_n+\frac{3}{2}Q\sin\psi_n,\\
g^n:&\\
& \psi_{n+1}=\psi_n+(\pi^{5/3}/Q)u_{n+1}^{-2/3} (\mbox{mod }2\pi),
\label{46}
\end{array}
\ee
 and one can write 
\beq u\rightarrow \pi\frac{u}{u_b}, \
\ Q \rightarrow \pi\frac{Q}{u_b}, \ \
u_b=\left(\frac{1}{2}\alpha^{3/2}N^2w\right)^{3/5}.\label{47}
\eeq

The map \(g^n\), or rather the family of maps depending upon the
parameter, becomes suitable for sequential analysis. So hereafter we
deal with high-frequency heating, relativistic particles, expressed in
terms of the discrete group \(g^n\), that acts on the smooth
manifold. Denote by

\be 
J=\frac{\partial
(u_{n+1},\psi_{n+1})}{\partial (u_n,\psi_n)} \label{Jacoby} 
\ee 

the Jacoby matrix of map (\ref{46}). It is important to note that the
Jacobian of (\ref{47}) is equal to one, therefore, \(g^n\) has a
structure of the differentiable area-preserving map. It stands to
reason that \(u\) and \(\psi\) are the canonical pair of variables.

\\ We now study the behavior of this dynamic system by computational
analysis. We have numerically integrated equations (\ref{46}) for
several different values of \(Q\) from \(10^{-3}\) to 0.1. Figure
\ref{f:1} shows some our results computed for map (\ref{46}) after
\((10^6\div 10^4)\) iterations for one trajectory in the
(\(u,\psi\))-phase space. The initial conditions were chosen in a
random fashion and correspond to the region of small values of
\((u,\psi)\).

\psfrag {u}[][b][0.8]{\(u\)}

\begin {figure}
\psfrag {a}[][b][0.8]{(a)}
\psfrag {b}[][b][0.8]{(b)}

\centering
\includegraphics[height=4cm]{map_r-q-0.03.eps}
\hfill%
\includegraphics[height=4cm]{map_r-q-0.021-arrows.eps}
\caption{ (a) One single trajectory of length \(10^6\)\ \ for \(Q =
0.01\pi\). The trajectory was started from the point \(u_0
=10^{-3},\ \  \psi_0 = 10^{-4}\). \ \ 
(b) The phase space of the map at \(Q = 0.007\pi \). The
identification of points is indicated by the arrows.}\label{f:1}
\end{figure}

Next we describe the structure of the point set.  Let us consider a pair
\((M,g^n)\), where \(M\) is a smooth manifold and \(g^n\) is the
differentiable area-preserving map (\ref{46}).  From (\ref{46})
follows that a pair \((M,g^n)\) is invariable because of the inversion
of a point with respect to a circle \(\psi (\mbox{mod} 2\pi)\in S \),
and the reflection of a point on a \(O\psi\) -- axis. On that basis, we
define the following equivalence rules:

\be
(-\pi ,u)\sim (\pi,u), \ \ \ (\psi,-\pi) \sim (\psi,\pi), \label{48}
\ee
where \(\sim\) stands for the equivalence sign.\\

The identification of these points is indicated by the arrows in
Fig.\ref{f:1}.

The equivalence rules (\ref{48}) allow us to represent phase space of
the system as the quotient space having the quotient topology of
a torus, \(M/\sim \simeq T^2, \ \ \ T^2=S\times S, \ \ \ \psi
(\mbox{mod} 2\pi) \in S, u\in S\)(Kosniowski, 1980).\\ 
Now we discuss the local topology considering a Jacoby matrix given by
(\ref{Jacoby}). Denote through \(\lambda_1\) and \(\lambda_2\) the
eigenvalues of the matrix \(J\). Recalling that local structure is
determined by neighbours of fixed points, we have from (\ref{Jacoby})
\be 
\det J= \lambda_1 \cdotp \lambda_2 =1, \ \ 
\ee
\be \mbox{tr}J= \lambda_1 + \lambda_2 = 2+ (\pi^5u^{-5})^{1/3},
\label{49} 
\ee 
where \( \det J\) and \(\mbox{tr}J\)  denote the
determinant and the trace of this matrix, respectively.\\ 
Arnold and Avez (1968) pointed out that the condition
\be |\,\mbox{tr}J\,|=3
\label{50} 
\ee 
implies that the manifold has the topology of a
hyperbolic torus. Thereupon from (\ref{50}) subject to (\ref{49}),
we find 
\be
 \lambda_1=\frac{3+\sqrt{5}}{2}, \ \ \  \lambda_2=\frac{3-\sqrt{5}}{2}.
\label{51} 
\ee 
Thus there exists such an exponent \(\lambda\) that
\(\lambda_1\) exceeds one, and the ratio \(\lambda_1/\lambda_2\)
is an irrational  number, therefore stochastic motion on  this
manifold is realized without fail. Present conditions (\ref{49})
and (\ref{50}) ensure that the relation 
\be
  |\,u_b \,|=\pi   \label{52}
\ee 
holds, and it determines  the upper bound of \(\{u\}\).
We know, that \(u\) depends on the parameters of this problem as
given by (\ref{45}). By that condition (\ref{52}) determines in
itself an equivalence class in parameter space.\\
Considering (\ref{45},\ref{47}) and (\ref{39}), from (\ref{52}) two
equations result 

\be
\varepsilon_{tb}=\left(\frac{1}{2}\alpha^{5/2}N^2w\right)^{1/5},
\label{53} \ee 
\be 
\varepsilon_b =
\left(\frac{1}{2}\alpha^{3/2}N^2w\right)^{2/5}= \alpha
\left(\frac{1}{4}bN^2\right)^{2/5}.
 \label{54} 
\ee

This agrees with both the numerical solutions and the results of
qualitative analysis.\\
According to Arnold and Avez (1968) the system (\ref{46}) is
classified as a \(C\)- system, and its statistical properties are
characterized by the positive Kolmogorov - Sinai entropy \( s=\ln
\lambda_1.\) The set \(\{u_n, \psi_n\}\) is generated by 
smooth map (\ref{46}); it has fractal structure, with the fractal
dimension defined by (Schuster, 1984)
\be d_F=1- \ln
\lambda_1/\ln \lambda_2. \label{55} 
\ee  
This yields \(d_F\)=2 for the  point
set of system (\ref{46}). It should be noted that
while smooth curves, piecewise smooth curves, and curves with finite
length all have fractal dimension one, the phase curves 
constructed  by iteration function systems (\ref{46}) all have fractal
dimension two.  
The visual impression of Figs.  \ref{f:1}  indicate dynamical chaos,
and numerical and analytical results confirm our intuitive
conclusion.

\section{Non-relativistic electron motion in a time-like wave packet}
We now turn to a description of the dynamics of non-relativistic
electrons in the field of the packet of whistler mode waves.\\
The Hamiltonian of the problem is: 

\be H(\bar{r}, \bar{p};t)=
\frac{1}{2m}(\bar{p}+\bar{A})^2, \label{N1} 
\ee 

where the vector potential is determined by expressions (\ref{5},
\ref{6}). We are going to consider the case, when the particle
velocity along the direction of an external magnetic field is smaller
than the group velocity of the wave packet. Therefore the wave field takes
the form of the time-like wave packet (Zaslavsky, 1985).  Then we
write down instead of (\ref{5})

\be
 \bar{A}^w= \left(A \sin\varphi \sum_{n \in Z}\delta (\tau -n),\
A \cos\varphi \sum_{n \in Z}\delta (\tau -n), \ 0\right),\ \ \varphi
= zk -t\omega, \label{N2} 
\ee 

where \(A, \omega, k \) are the
magnitude, frequency and wave number of the characteristic mode in
the wave packet, \(T\) is the timescale, \(\tau = (t/T) (\mbox{mod}
\ 1),\) \ \ \( \delta_n = \delta (\tau -n),\ \ \delta(\cdotp)\) is
the Dirac delta function. Now the Hamiltonian expressed in terms of the
action-angle coordinates given by (\ref{10}, \ref{11}) becomes 

\be
H(z,p;\, \theta, I;\, t) =
H_0(p,I)+\sqrt{\frac{2\omega_BI}{m}}A\cos\psi \sum_{n \in
Z}\delta_n,\label{N3} 
\ee 
\be 

H_0=\frac{p^2}{2m}+ \omega_BI,
\label{N4} \ee \be \psi = zk+ \theta -t\omega. \label{N5} 
\ee 

The equations of motion associated with (\ref{N3}) become

\be
\dot{p_z}= k\sqrt{\frac{2\omega_B I}{m}}A\sin\psi\sum\delta_n, \ \ \
\dot{I}= \sqrt{\frac{2\omega_B I}{m}}A\sin\psi\sum\delta_n
\label{NN1} 
\ee 
\be \dot{z} = v_z = p_z / m, \ \ \ \dot{\theta} =
\omega_B. \label{NN2} 
\ee 

We retain in (\ref{NN2}) only the leading terms on condition that
\(A/mv_t\) is a small parameter of the problem.  In order to match the
resonance conditions \(kv_z + \omega_B - \omega = 0\), electrons in
general must move in the opposite direction as waves since the whistler
wave frequency \(\omega\) is smaller than the electron cyclotron
frequency \(\omega_B\). In view of this fact and also taking into account
 the relation \(\omega_BI = mv_t^2/2\), we write the
invariant of Hamiltonian flow (\ref{NN1}, \ref{NN2} ) as
\beq 
v_t^2 = \alpha(v_0
- |v_z|), \alpha = 2 (\omega_B / \omega)v_{ph} \label{35.2} 
\eeq
where  \(v_0 ( = v_{gr})\) is the speed of an electron along an
ambient magnitude field.

Using (\ref{35.2}) and the expression \(mv_t^2 / 2 =
\omega_BI\) once more to represent the equations of motion in the
form

\beq
\dot{v_t} = \omega_B w \sin \psi \sum_{n \in Z}\delta(\tau-n),
\nonumber\\
\dot{\psi} = \omega_B  - \omega - kv_0 + (k / \alpha)v_t^2,
\label{35.4} 
\eeq
 
where we have again used the notation 
\be w = A /
m = \frac{1}{2}\alpha b, \ \ b =B^w / B. \label {35.5} 
\ee 

By integrating these equations can be transformed into a map, \(G^n\),
\be 
G^n:
\begin{array}{ll}
u_{n+1}    = u_n    + Q\sin\psi_n,\ \ \ \ \{u \in U \subset R|\ \ 0 < u \le 1\},\\
\psi_{n+1} = \psi_n + su_{n+1}^2 \ \ (\mbox{mod  }2\pi), \label{35.6}
\end{array}
\ee 
where subscript \(n\) refer s to values  taken at time \(t =
nT.\) Here one has introduced the new variable, \(u\), 
\be u = v_t
/v_b,  \ \ v_b = \sqrt {\alpha v_{gr}}, \label{35.7} 
\ee 
\(v_b \) is
all perpendicular (to \(B\))  speed of electrons,
 and the new parameters are
\be 
Q = N w / v_b, \ \ \ N = [\omega_B T], \ \ \ s = 2N(v_b /
\alpha)^2, \label{35.8} 
\ee 
where \(w \) and \(\alpha\) are given by
(\ref{35.5})  and [\ \ ] stands for the integer part.

Thus the Jacobian of the matrix of \(Q^n\), det\(J = 1\), so \(Q^n\) is
the measure-preserving map, and \(\psi\), \(u\) are the symplectic
pair on the smooth manifold, \(M: \ U \times S\).

Thereby we apply the condition (\ref{52}) to (\ref{35.6}), to find
the expression 
\be 
u_c = inf\{u\} = (2Qs)^{-1} = \alpha
\left(2N^2bv_b\right)^{-1}, \label{35.9} 
\ee 
that well determines the lower bound \(inf\{u\}\) of the stochastic set.

We have numerically  integrated equations (\ref{35.6}) for values of
\(Q\) from 0.0005 to 0.01. Our results are shown in Fig. 2.

\begin{figure}
\psfrag {a}[][b][0.8]{(a)}
\psfrag {b}[][b][0.8]{(b)}
\psfrag {c}[][b][0.8]{(c)}
\centering
\includegraphics[height=4cm]{nonr-q-1.e-3-u-1.e-2-s-5.e3.eps}
\hfill%
\includegraphics[height=4cm]{nonr-q-1.e-3-u-1.1e-1-s-5.e3.eps}\\

\includegraphics[height=4cm]{nonr-q-0.5e-3-u-2.1e-1-s-5.e3.eps}\\
\caption{ The phase space \(G^n\). (a) A single trajectory of 
  length \(10^6\) for \(Q = 10^{-3}\),
\(s = 5\cdot 10^{3}\), \((u_0,\psi_0) = (0.01, 0.0001).\)\ \ 
(b) A
phase curve for the same values of \(Q\) and \(s\) as in (a). The
trajectory was started from the point (\(0.11, 0.0001\)). (c) One
trajectory  of length \(7\cdot 10^6\) for \(Q = 0.5\cdot 10^{-3}\) \(s=5
10^{-3}\), (\(u_0, \psi_0) = (0.21, 0.0001)\)}
\label{f:2}
\end{figure}


 The
overall picture of the phase space is quite different for \(u <
u_c\) and \(u \ge u_c.\) In the first case, the motion is regular. Then
the figures indicate the existence of a threshold for the initial
particle velocity  above which the trajectory becomes chaotic and the
extent of the stochastic region as expected, increases linearly with
\(Q\). Nonlinear electron acceleration by a wave packet of the
whistler mode  waves is therefore always a stochastic process. By
just the same operational mode as above, it is easy  to see that the
\(\{u, \psi\}\) point set, which is generated by the smooth map
\(G^n\), has fractal structure with a fractal dimension equal to two,
and its topological entropy is given by the expression \(s = \ln((3
+ \sqrt{5}) / 2).\)


\section {Diffusion evolution}

\indent Particle dynamics in random electromagnetic fields is known
to be described by the quasilinear theory (QLT) (Vedenov et al.,
1962; Drummond and Pines, 1962). The QLT approach in particular,
has been employed to deal with diffusion  of electrons in the
turbulent fieldi of whistler wave packets (Kennel and Petchek, 1966).
Theoretical work has been done to study the wave-particle
interaction in rather strong random wave fields (Dupree, 1966;
Weinstock, 1966). This theory is well-known as the
Dupree-Weinstock approach or the resonance broadening theory.

In coherent electromagnetic fields, on the other hand, the particle
dynamics are not described by these theories. The nature of
diffusion in this case is described by stochastic dynamics of particles, when
the motion along stochastic trajectories gives rise to the so-called
deterministic diffusion (Lichtenberg and Lieberman, 1983). General
mathematical and physical aspects of this problem have been
discussed, for instance, in reviews (Arnold and Avez, 1969;
Chirikov, 1979).

The purpose of the present paper is to investigate electron dynamics
in coherent packets of whistler mode waves. We have employed here the
method in which the perturbation of particles orbits are included
in a Fokker-Planck-Kolmogorov (FPK) equation (Zaslavsky, 1985).

In the present section, we shall exploit the canonical Hamiltonian
structure that has been developed above. Thus, the drift kinetic
equation follows at  once since the canonical structure equations of
motion (\ref{35}, \ref{36}) have been established. The distribution
function (probability density)\(f(u;t)\) obeys the Fokker - Planck -
Kolmogorov (FPK) equation 
\be 
\frac{\partial f(u;t)}{\partial
t}=\frac{1}{2}\frac{\partial}{\partial u}D \frac{\partial
f}{\partial u}, \label {56} 
\ee
which holds if \(Q \ll 1\) (Zaslavsky,  1985).\\
 Here \(D\) is the conventional diffusion coefficient in  phase space,
\be
 D=<(u_{n+1}- u_n)^2>T^{-1}, \label{57} 
\ee 
in which \((u_{n+1}-
u_n)\) is substituted from (\ref{46}), and \(<\cdotp>\) denotes the
phase average, \(T\) is the timescale of mapping (\ref{46}). The
function \(f(u,t)\) belongs to the space of all differentiable
functions supported in \([-\pi,\pi]\), and the function 
\be
\int_{-\pi}^{\pi} f(u,t)du=1 \label{58} 
\ee
is the condition of normalization.\\
First, by means of (\ref{46}) we calculate using (\ref{57}) the
diffusion coefficient 
\be 
D=Q^2/2T.  \label{59} 
\ee 
Then, making use of the result proven above (\ref{52}) along with
(\ref{56}) and (\ref{59}), we evaluate the characteristic time for
redistribution \(u\) over the spectrum
\be 
t_d \simeq u_b^2/D= 2T(\pi/Q)^2. \label{60} 
\ee

A objective of this study is to determine  the time-independent
distribution function and the rate  of heating.  This requires a
solution of the FPK equation together  with a normalization
(\ref{58}). Thus a solution of the FPK equation along with the
boundary condition \(f(-\pi)=f(\pi) \) in the limiting case \(t\ge
t_d\) may be given in the form of the uniform distribution 
\be
 f(u)=(2\pi)^{-1}.  \label{61}
\ee 
This result is supported by the simulations of Fig. 3.
The
number of \(n\) imaging points in \(n\) elements of phase space,
\(\Delta \psi \Delta u\), represents the 
joint probability density, \(F(\psi,u)\).   We partition each axis between external points so that the cell size is identical and equal to 
\((2\pi
/10)^2\).

\begin{figure}
\centering
\includegraphics[height=5cm]{map_r-d3-q-0.03.eps}
\caption{ Joint distribution (density of \(n\) imaging points)
(\(F(\psi,u)\)) obtained numerically from (\ref{46}) with the
parameters as in Fig. \ref {f:1}. }
\label{f:3}
\end{figure}

\begin{figure}
\psfrag {a}[][b][0.8]{(a)}
\psfrag {b}[][b][0.8]{(b)}
\psfrag {n}[][b][0.8]{\(n,10^4\)}


\includegraphics[height=4cm]{sqr-energy-q-0.03.eps}
\hfill%
\includegraphics[height=4cm]{sqr-energy-q-0.021.eps}\\

\caption{ \(u^2\) as a function of \(n\) obtained numerically
from (\ref{46}). (a) \(Q = 0.01\pi\). (b) \(Q = 0.007\pi\).}
\label{f:4}
\end{figure}


Now, it becomes relevant to determine how the system evolves in time
at \(t\le t_d\). We exploit the FPK equation with \(f(u)\) and require its
derivative \(\partial f/\partial u\) to vanish at the boundary. We
introduce the moment \(<u^2>=\int_{-\pi}^{\pi}du u^2 f(u),\) multiply
equation (\ref{56}) by \(u^2\), and integrate the resulting equation
over \(u\) to obtain
\be
 \frac {d<{u}^2>}{dt}=D. \label{62}
\ee 
\(u^2\) as a function of \(n\) is shown in Fig. 4.  The mean square
velocity grows linearly with time, as predicted by (\ref{62}). The
duration of the initial regime depends on \(Q\) as described by
(\ref{60}).

\noindent In the relativistic limit it is easy to see that
\(f(\varepsilon) = f(|\varepsilon_z|)\), or rather, the difference
between \(f(\varepsilon)\) and \(f(|\varepsilon_z|)\) is as small as
\(1/(2\varepsilon^2)\).  Then equation (\ref{46}) and the
measure-preserving point transformation

 \be
f(\varepsilon)=\frac{1}{2}f(u)(du/d(\varepsilon)) \label{64} 
\ee
along with the normalization
\[\int_0^{\varepsilon_b}f(\varepsilon)\,d\varepsilon =1\]
allow us derive the time-independent distribution 
\beq
f(\varepsilon)d\varepsilon =
\frac{3}{2}\varepsilon_b^{-3/2}\varepsilon^{1/2} d\varepsilon, \ \ \
\{\varepsilon \in R_+ \ |0 \leq \varepsilon \leq \varepsilon_b \},
\label{65} 
\eeq 
\noindent  and the nonlinear FPK equation describing the
evolution of an energetic spectrum in the allowed range of particle
energies 
\beq \frac{\partial}{\partial t}f(\varepsilon, t) &  = &
              \frac{\partial}{\partial \varepsilon}
             D(\varepsilon)
         \frac{\partial}{\partial \varepsilon} f(\varepsilon, t),
\label {in-82} \\
 D(\varepsilon) & = & \frac {\alpha ^3 N^2 b^2}{16 \varepsilon}.
\label {in-83} \\
\eeq 
\noindent One may observe that there is an explicit dependence
of \(D\) on \(\varepsilon^{-1}.\) This distinct feature is one reason
that particle energy, on the average, grows slowly in time, as
\(t^{1/3}\) as predicted by the solution of (\ref {in-82}), or, in
other words, the heating rate of electron decreases with its energy
as \(\varepsilon^{-2}\). Then from (\ref{in-82}), (\ref {in-83}) and
(\ref{56}) we find the characteristic time to establish the
distribution (\ref {65}) 
\be 
t_d = 16 T \left (N / 64
b^2\right)^{2/5}. \label {in-84} 
\ee

The distribution function \(f(u,t)\) in the non-relativistic case
also obeys the FPK equation as a consequence of the fact that the
equations of motion, expressed in coordinates \((u,\psi)\) reflect
the underlying canonical  Hamiltonian structure of the equations
(\ref{35.2}). By reason of (\ref{35.9}), we write the normalization as

\be \int_{u_c}^1f(u)\,du =1. \label {new-51} 
\ee 
Then by means of (\ref{35.8}) along with (\ref{35.9}) and
(\ref{35.8}), we calculate by formulas (\ref{57}) and (\ref{60}) the
diffusion coefficient
\be D = \frac{Q^2}{2T} =
\frac{1}{2T}\frac{N^2w^2}{v_b^2}, \label {new-52} 
\ee 
and the characteristic time for redistribution \(u\) over the spectrum

\be 
t_d = 2T/Q^2 = 2Tv_b^2/N^2/w^2. \label {new-53} 
\ee 
Next using the normalization for \(f(u,t)\) in the FPK equation, it follows
when
\(t > t_d\) the distribution function  can be given by 
\be f(u)\,du
= (1 - u_c)^{-1}\,du. \label {new-54} 
\ee

\noindent This means that the random variable  \(u\) is evenly
distributed on \([u_c,\ \ 1]\), a result supported by the
simulations in Fig. 5. As will be shown below, typically  \(u_c \ll
1;\) and in the following we write simply \(f(u)\,du = 1\,du.\)


\begin{figure}
\centering
\includegraphics[height=5cm]{map_nonr-d3-q-1.e-3-u-0.11-s-5.e3.eps}
\caption{ Joint probability density in (\(u, \psi\)) phase
space. The distribution is well defined  by \(1 /2\pi(1 -u_c)\),
except for a boundary effects.}
\label{f:5}
\end{figure}



\begin{figure}
\psfrag {n}[][b][0.8]{\(n,10^5\)}

\centering
\includegraphics[height=4cm]{sqr-energy-nonr-q-1.e-3-u-0.11-s-5e3.eps}
\caption{Time series of \(u^2\) obtained numerically from
the map \(G^n\).}
\label{f:6}
\end{figure}


\noindent Now from the transformation \(f(u)\,du = f(v)\,dv\) and
the relation \(u = v_t / v_b\) we obtain 
\be 
f(v_t) = v_b^{-1}.  \label{new-55} 
\ee 
Taking account of (\ref{new-51}), by (\ref{new-55}) we derive the
quadratic mean of \(v_t\)
\be <v_t^2>
\int_0^{v_b}f(v_t)v_t^2\,dv_t = \frac{1}{3}v_b^2 =
\frac{4}{3}v_{ph}^2\frac{\omega_B - \omega}{\omega}, \label {new-56}
\ee 
the mean value of \(|v_z|\) and the mean velocity 
\be <|v_z|> =
\frac{2}{3}v_{ph}\frac{\omega_B - \omega}{\omega_B}, \label {new-57}
\ee 
\be <v_z> = 0. \label {new-58} 
\ee

Then we define by 
\be T_t = <\frac{mv_t^2}{2}>, \ \ \ T_z =
<\frac{mv_z^2}{2}>, \label {new-59} 
\ee 
to find the ratio 
\be
\frac{T_t}{T_z} = \frac{5\alpha}{8v_{gr}} = \frac{5}{8}
\frac{\omega_B}{\omega}\,\frac{\omega_B}{\omega_B - \omega}, \label
{new-60} 
\ee 

describing the anisotropy of the energy distribution (\ref
{65}) in \((v_z, v_t)\) phase space. This is an anisotropic
distribution function with its derivative 
\be \frac {dv_t}{v_z} =
-\frac {1}{2}\sqrt {\frac {\alpha}{v_{gr}
                    - |\,v_z\,|}}sgn v_z, \label {in-95}
\ee

that tends to infinity as \(|\,v_z\,|\rightarrow v_{gr}\), with 
standard deviation
\be
\sqrt{<v_z^2> -<v_z>^2} = \sqrt {8 / 15} v_{gr}, \label {in-96}
\ee
that describes the so-called pancake distribution in \((v_z,\,v_t)\)
phase space.


Again it is useful  to evaluate the energy distribution function,
\(f(\varepsilon)\). Considering (\ref{35.2}), we note that the
particle energy \(E\) is a specified function of \(v_t\). Thus the
measure-preserving transformation \(f(\varepsilon)\,d\varepsilon =
f(v_t)\,dv_t\) determines this problem completely, subject to
appropriate boundary conditions. To that end, we have to find the
lower \((\varepsilon_t)\) and upper \((\varepsilon_z)\) boundaries
of an energetic spectrum. At first, from (\ref{35.2}) we get 

\be
v_{th} = |\,v_z\,|_c = v_{gr} - v_c^2/\alpha. 
\ee 

Then \(\varepsilon_c\) and \(\varepsilon_0\) are given by
 
\be
\varepsilon_c = \frac{1}{2}\left[\left(v_{gr} -
v_c^2/\alpha\right)^2 + v_c^2\right] \simeq \frac{1}{2}v_{gr}^2,
\label {new-61} 
\ee 

\be 
\varepsilon_0 = \frac{1}{2}v_b^2 =
\frac{1}{2}\alpha v_{gr}. \label {new-62} 
\ee 

The corresponding solution for \(f(\varepsilon)\) provided by
(\ref{new-55}) can be expressed in the form
\be 
f(\varepsilon)\,d\varepsilon =
\frac{1}{2}\frac{1} {\sqrt{(\varepsilon_0 -
\varepsilon_c)(\varepsilon - \varepsilon_c)}} \,d\varepsilon, \ \ \
\varepsilon_0 > \varepsilon
> \varepsilon_c. \label {new-63} 
\ee

A solution of this form relates directly to the behavior of the system near
the order-chaos bifurcation transition.  The function
\(f(\varepsilon)\) undergoes a sudden change at \(\varepsilon =
\varepsilon_c\), and in the region of regular motion \((\varepsilon <
\varepsilon_c, \varepsilon = \mbox {const})\) has the form of the
Dirac \(\delta - \) function, \(f (\varepsilon) = 8 (\varepsilon -
\mbox {const}).\) Using the same pattern as
above, we evaluate the heating rate

\be 
<\dot{u}^2> = \frac{Q^2}{2T}, \label {new-64} 
\ee 

which in the explicit form may be written as

\be 
<\dot{\varepsilon}> =
<\frac{1}{2}v^2> = \frac{1}{8T}\alpha ^2b^2, \ \ \ \varepsilon_c
< \varepsilon < \varepsilon_0. \label {new-65} 
\ee 

This result is nontrivial, because chaotic motion in (\(u,\psi\))
- phase space leads to important and easily observable
macroscopic effects such as the stochastic heating of plasma particles.

The nonlinear difference equations (\ref {35.6}) have been solved
numerically, and \(u^2\) as a function of \(n\) is shown in Fig.6.
\(u^2\) reaches a plateau value in time, which depends on \(Q\) as
described by (\ref {new-53}). The mean square velocity grows linearly
with time, as predicted by (\ref{new-64}).

Now we describe the effects associated with stochastic heating of 
high-energy particles.
First, we discuss the pitch angle distribution in a
range of non-relativistic energies.

Denote by \(\lambda_p\) the pitch angle, and let \(\chi\) be its
complementary angle. Then we have
\be 
\tan \chi = v_z / v_t, \ \ \
\chi \in (-\pi / 2, \pi / 2). \label {new-66} 
\ee 

One can use equations (\ref {35.2}) to represent \(\tan\chi\)
as a function of \(u(u= v_t/v_b),\) namely,

\be 
\tan\chi =
\sqrt{\frac{v_{gr}}{\alpha}}\cdot \frac{1-u^2}{u}. \label {new-67}
\ee 

And again, we use the measure-preserving point transformation \(f(u)\,du =
f(\chi)\,d\chi\) to find the angle distribution

\be
f(\chi) = \sqrt{\frac{\alpha}{v_{gr}}}\left(\frac{u^2}{1 + u^2} +
\frac{v_{gr}}{\alpha}\cdot\frac{(1-u^2)^2}{1 + u^2}\right). \label
{new-68} 
\ee 

Equation (\ref {new-68}) determines the angle distribution function in
a parametric form. Nonlinear transformation (\ref {new-67}) allows us
to express (\ref {new-68}) as a function of a proper variable
\(\chi\).

Note that \(f(\chi)\) is invariable under transformation 
\be 
u \
\rightarrow - u, \ \ \ \chi\ \rightarrow - \chi; \ \ f(\chi) =
f(-\chi), \label {new-69} 
\ee 
therefore \(f(\chi)\) is a symmetric function on \((-\pi / 2, \pi /
2).\) However, the explicit expression for \(f(\chi)\) has a very
complicated form. The following asymptotic formulas are valid:

\be 
f(\chi) \approx
           \frac{1}{2}\sqrt{\frac{\alpha}{v_{gr}}}\left( 1 -
\sqrt{\frac{\alpha}{v_{gr}}}\tan\chi\right) \mbox{     as   } \chi \
\ \rightarrow 0, \label {new-70} \ee \be f(\chi) \approx \sqrt
{\frac {v_{gr}}{\alpha}}(1 + \tan^{-1}\chi), \mbox{ \hspace {1cm} as
} \chi\ \rightarrow \pm \pi / 2,  
\label {new-71} 
\ee

\be 
f(\chi) \approx  \sqrt {\frac {v_{gr}}{\alpha}} \left(1 + 2
\left(\frac {\varepsilon}{\varepsilon_c} -1 \right)\right)
 \mbox{ \hspace {1cm} as } \varepsilon
\downarrow \varepsilon_c, \ \ \chi\ \rightarrow \pm \pi / 2,  
\label {new-72} 
\ee

\be 
f(\chi) \approx \frac{1}{2}\sqrt{\frac {v_{gr}}{\alpha}}
\left((\frac {\varepsilon}{\varepsilon_c} -1\right) \mbox{ \hspace
{1cm} as   } \varepsilon \uparrow \varepsilon_0, \ \ \chi\
\rightarrow 0. \label {new-73} 
\ee 
where \( \varepsilon_c\) and
\(\varepsilon_0\) are given by (\ref{new-61}) and
(\ref{new-62}).\\

 Making use of these equations, we
conclude that the function \(f(\chi)\) governing the distribution of 
angles
is a convex symmetric function of \(\chi\), having the maximum at
\(\chi = 0  \ \ (\chi_p = \pi / 2\)), and the relation between the
scattering degrees \(f(\chi = 0)\) and \(f(\chi = \pm \pi/2)\) is
given by \(f(0)/f(\pm \pi/2) = \varepsilon_0 / 2\varepsilon_c\). The
result indicates also that particles are equally likely to be
scattered  in the direction with or against the wave. To this may be
added that the given function is an increasing function of the particle
energy on \((\varepsilon_c, \ \varepsilon_0)\), and it hinges
asymptotically on \(\varepsilon\) as \((\varepsilon -
\varepsilon_c)\) near the energetic threshold, and as \(\sqrt
{\varepsilon / \varepsilon_c} - 1\)  at \(\varepsilon > 2
\varepsilon_c.\)

Let us discuss particle scattering in a range of relativistic
energies. Conforming to \( A/m \ll 1\), we write
\be 
\tan \chi = \frac {v_z}{v_t} \approx \frac {p_z}{p_t}
\approx \frac {\varepsilon_z}{\varepsilon_t}, \ \ \ \chi \in (-\pi /
2, \ \ \pi / 2). \label {new-74} 
\ee 
at energies above 1~MeV.

\noindent Taking account of the invariant of motion (\ref{39}) we derive
the following dependence 
\be 
\tan \chi = (\varepsilon_z /
\alpha)^{1/2}. \label {new-75} 
\ee 
Using (\ref{65}) along with the measure-preserving
transformation \(f(\chi)\,d\chi = f(\varepsilon)\,d\varepsilon\) yields
the time-independent pitch angle distribution \(f(\chi)\), namely,
\beq 
f(\chi)\,d\chi &= &3 (\alpha
/ \varepsilon_b)^{3/2}\tan^4\chi d\,\chi,
\label {new-76}\\
\chi \in (-\chi_m, \ \ \chi_m), & & \tan \chi_m = (\varepsilon_b /
\alpha)^{1/2}, \nonumber, 
\eeq 
and its dependence on the particle energy
\be 
f(\chi) = 3 \frac
{\varepsilon^2}{\varepsilon_b^{3/2}\alpha^{1/2}}. \label {new-77}
\ee 
The function given by (\ref{new-76}) is the concave, symmetric
function of \(\chi\), which has a maximum at \(\chi = \pi / 2\
(\chi_p = 0)\), and derivatives that tend to infinity as  \(\chi
\rightarrow \pm \pi / 2\). This function describes the so-called 
U-like distribution.

\noindent Comparing (\ref{new-77}) and (\ref{new-72}), we note that
in the latter, dependence on \(\varepsilon\) implies that \(f(\chi)\)
changes more slowly.

Next we show that the stochastic  heating is accompanied by a radial
diffusion of particles  in space. Indeed, subject to (\ref{11}) from
(\ref{9}), which relates directly to the local spatial coordinates,
the following expression results 
\be 
\frac{d\,<{r}^2>}{d\,t} =
\omega^{-2}_B <\dot{v}_t^2>, \label {new-78} 
\ee 
where \(<\dot{v}_t^2>\) is given by (\ref{new-65}).

\noindent Then we evaluate the radial drift velocity by 
\be v_d =
<\dot{r}^2> / <r>, \ \ <r> = <v_t>\omega^{-1}_B, \label {new-79} 
\ee
where \(<r>\) is the mean gyroradius, and \(<v_t> = \frac{1}{2}v_b\)
is the mean perpendicular velocity.

\noindent Via (\ref{new-78}), we write down the expression for the
drift velocity 
\be
 v_d = N^2w^2 / v_b\omega_BT, \label {new-80} 
\ee
which is true at \(t < t_d\).



\section {Application}
So, we point out that our result can be quite useful to the
understanding of wave phenomena in the Earth's magnetosphere. At first
the parameters of the problem are chosen to reflect a typical
whistler wave propagating along the \(L_p = 3\) shell, where the 
electron cyclotron frequency \(\omega_B = 2 \cdot 10^5\)~s\(^{-1}\),
the electron plasma frequency \(\omega_p = 6 \cdot 10^5\)~s\(^{-1}\),
and the wave frequency \(\omega = 9 \cdot 10^4\)~s\(^{-1}\). Note that
without doubt, the dispersion relation (\ref{9}) corresponds to the
problem. By that and through (\ref{9}) we find the phase velocity,
\(v_{ph} \simeq 0.17,\) the wavenumber \(k = \omega / v_{ph} (\simeq
1.8\cdot 10^{-5}\)~cm\(^{-1})\) and the group velocity, \(v_{gr} = 2
v_{ph} (\omega_B - \omega / \omega_B)\), \(v_{gr} \simeq 0.24\).
Then we employ expression (\ref{35.7}) to compute the speed of
particle, \(v_b = \sqrt{\alpha v_{gr}} \simeq 0.45\).  Since the
zone of intense wave activity is spatially localized, we have
accepted that the lengthscale \(L\) and Earth's radius \(R_E (\simeq
10^9\)~cm) are of the some order of magnitude.  It should be noted
that the geometry of an ambient magnetic field does not, in essence, play
a role in the given problem, thus the lengthscale  of an interaction
region is typically smaller than the travelling path, i.e. an
ambient magnetic field is always locally an uniform field. It is
clear, on the other hand, that the conditions \(nL \gg 1\) and
\(\omega_B T \gg 1\) must be justified. To this may be added the
analytic model describing the motion of a particle out of the interaction
region.  This motion  must include the radial diffusion of electrons, 
thus conserving the
first and second adiabatic invariants.   Obviously, the
timescale \(T\) of the problem equals the transit time of the
particle through the packet, i.e. \(T = L / v_{gr}\).  By that, the
value of \(T\) approximately equals \(1.4\cdot 10^{-1}\)~s. We
will use this value to obtain a quantitative estimate of the quantity
\(N = [\omega_B T], \ \ N = 3\cdot 10^4\). Remembering that \(A^w\)
is related to \(B^w\) by (\ref{3}), we derive \( w = B^w / km\). To
calculate this parameter we take \(B^w = 10^{-6}\)~G, to find \(w =
10^{-7}.\) Then by means of (\ref{35.9}) we calculate \(v_c \simeq
2.5\cdot 10^{-3}\). Applying this result to (\ref{35.2})
yields an estimate of \(v_{th}\), close to \(v_{gr}\). This implies that
when the speed of an electron exceeds the threshold velocity \(v_{th}\),
the electron demonstrates chaotic motion. Again, one uses the 
relation (\ref{9}) and formula (\ref{new-57}) to calculate the
characteristic values of the means \(<|v_z|>= 0.22, \ \ <v_t> = 0.32,\)
the degree of anisotropy \(T_t/T_z \simeq 2.5\), and the energy flow
along the direction magnetic field, \(Q_T = <v_z><E_p>, \ \ Q_T
\simeq 2.4 \cdot 10^{14}\)~eV\,cm\,s\(^{-1}\). At last, by formula
(\ref{new-80}) we evaluate the value of radial drift velocity, \(v_d
= 3\cdot (10-10^3)\)~cm\,s\(^{-1}\), and by  (\ref{60}) the
characteristic time of the diffusion process, \(t_d \simeq 1\) min.
 Now
we apply the result of our 0 to the magnetosphere of the Earth near the
plasmapause at \(L_p = 4.\) and use the following values:
 \(\omega_B = 6\cdot 10^4\)~s\(^{-1}\),
\(\omega_p = 12 \cdot 10^4\)~s\(^{-1}\), \(\omega = 2.4\cdot
10^4\)~s\(^{-1}\). Thereafter we obtain \(v_{ph} \simeq 0.25, \ \
v_{gr} \simeq 0.33, \ \ k \simeq 3.2\cdot 10^{-6}\)~ cm\(^{-1}\) and
\(v_b \simeq 0.62.\) Then
we calculate \(N = 2 \cdot 10^3\),
\(w \simeq 10^{-6}\), \(v_c \simeq 3 \cdot 10^{-3}\) and \(v_{th}
\leq v_{gr} (\simeq 0.33)\). Finally, the upper \((E_0)\) and lower
\((E_{c})\) values of the particle energy are estimated respectively,
at \(v = v_b\) and \(v = v_{th}\) for non-relativistic electrons to
be \(E_{c} \simeq 30\)~keV and \(E_{0} \simeq 120\)~keV.  Thus the
motion along stochastic trajectories gives rise to significant
energy diffusion, which extends to electron energies near 100~keV.
This allows whistler waves to contribute to the stochastic
acceleration of electrons in the high-energy tail of the plasma
population. Acceleration of the relativistic electrons is described
by equations (\ref{44}) and
(\ref{46}) in the ultra-relativistic limit. Solutions of these
equations show that the energy of 0 particles is bounded from above as
given by (\ref{54}). Now let us estimate the wave amplitude,
\(B_c^w\), required for strong diffusion scattering of relativistic
electrons. Considering (\ref{54}), we obtain 
\beq 
b_c = 4\alpha^{-5/2} N^{-2}, \label{Ap:1} 
\eeq 
where \(N = [kL]\) is the characteristic number of modes in a packet,
\(k \simeq 1.8 \cdot 10^{-5}\)~cm\(^{-1}\), \(\alpha \simeq 0.85\),
and the lengthscale \(L\) is the scale for the region of interaction,
which is approximately the Earth radius, \(R_E(\simeq 10^9\)~cm).

\noindent Utilizing this results in (\ref{Ap:1}) and we obtain an
estimate of \(B_c^w \simeq 7 \cdot 10^{-8}\)~G. The level
of whistler waves is fairly low, because from observational
information on wave activity we know that typical wave amplitudes
of different whistlers are in the range of \((10^{-6} \div
10^{-5})\)~G (Nagano et al., 1996; Parrot and Gaye, 1994).  Thus, the
upper value of the energy spectrum is calculated from (\ref{54}) and 
approximately equals 5~MeV provided \(b = 10^{-5}\).
 Now we estimate the
characteristic time of the redistribution over the spectrum in a
relativistic particle energy range, \(t_d \simeq (1-10)\)~hours.
The duration of this process significantly exceeds that derived from 
non-relativistic particle energies. Observations show that the power of 
whistler-mode
waves is enhanced at substorm injection and decays over a
period of hours or over a period of a few days during the storm
recovery (Horne  et al., 2003). \noindent Now to see that our
approach correctly interprets the facts, it is necessary to find
the range of validity of our results. So, in the non-relativistic
case the particle energy given by \(E_p = mv_b^2 / 2\) is 
small as compared to the rest energy of electron. According to
(\ref{9}) and (\ref{new-62}), this ratio may be evaluated as
\beq 
\frac {E_0}{m} = \frac {\omega_B}{\omega} \cdot v_{ph}\cdot
v_{gr} = 2\left (\frac{(\omega_B - \omega)}{\omega_p}\right )^2.
\label{Ap:2} 
\eeq 
Because the local electron gyrofrequency \(\omega_B\) in Earth's
magnetosphere never exceeds the corresponding plasma frequency
\(\omega_p\) except for in the low density auroral cavity, the right
hand side of (\ref{Ap:2}) is sufficiently small compared with unity,
therefore, \((E_p / m) \ll 1\) is valid.

The equations of motion (\ref{44}) in the region of relativistic
energies are obtained in the 0 representing the ultra-relativistic
approximation, \(\varepsilon^2 \gg 1\), which has been employed in
deriving (\ref{46}). We must also establish that
expression (\ref{54}) actually corresponds to the problem,  if
requiring \((\varepsilon_t / \varepsilon_z)^2 \ll 1\) holds.
According to (\ref{39}), this inequality may be written in the form
\beq 
2\frac {\omega_B}{\omega} v_{ph}\varepsilon^{-1} \ll 1. 
\eeq
Owing to (\ref{54}) the latter can be transformed into the requiring
\beq 
B^w / B \gg 4N^{-2}. \label{Ap:3} 
\eeq 
Recalling direct
amplitude measurements ( Nagano et al.,  1996) and the extremes of
spectral density measurements (Parrot and Gaye, 1994) yield \((B^w /
B) =  (10^{-4} \div 10^{-3})\), we conclude that (\ref{Ap:3}) is
trivially satisfied, thus \(\varepsilon \ll 1 \) is a fortiori valid;
therefore, the upper limit for the energetic spectrum given by
(\ref{54}) is correct.

These results demonstrate the possibility of electron
heating over a wide energy range between \(\geq 10\)~keV and a few  MeV.
This is in reasonable agreement with the experimental observations.

\section {Summary}
A canonical Hamiltonian approach has been employed to deal with the
whistler-wave, electron interaction and the stochastic heating of
high-energy electrons in magnetized plasmas.  Group symmetries
of the problem allow us to represent the equations of motion as 
measure-preserving  maps,  which describe particle motion
transmitted through a coherent packet of whistler waves. The main
properties of stochastic sets generated by the mapping equations
are studied in detail, analytically and numerically.    The
fractal dimension is defined to characterize these sets and measure
their roughness.

\noindent We show that under conditions typical for the wave
electron acceleration mechanism, the heating region is determined by
dynamic chaos in phase space and the the boundaries of that region. The
heating rate is governed by the nature of the kinetics, which in
turn depends on the canonical variables. Theory has shown that an
energetic threshold exists below which  electron motion is
regular.  When  the initial energy  is greater than the energetic 
threshold,  an electron moves
stochastically. Thus, nonlinear electron acceleration by a whistler
wave packet is always a stochastic process.  In the
ultra-relativistic limit when the phase velocities of the waves in
the packet differ significantly from the speed of the particle, the energy
of the  particle is bounded from above. These upper and lower values
of particle energy are calculated and their dependence on  
parameters of problem is determined.

Chaotic dynamics of particles in a coherent packet of whistler waves
is investigated by the FPK equation for a distribution  function  in
the corresponding phase space. In the present work, we have also
described the effect of stochastic  heating on the pitch angle
distribution and the transverse (radial) drift of high-energy electrons.

\noindent The above results afford a basis for conclusions. So, if the
condition \(E > E_th\) is valid, an electron with small speeds (\(
v_z  < v_{gr}\)) is accelerated, and its energy increases as
predicted by (\ref{new-65}). Then, when the particle velocity
exceeds \(v_{gr}\), a process such as stochastic heating is described by
(\ref{in-82}).  As a result, any distribution is established 
in an allowed range of
particle energies.  The energy spectrum in a range of
non-relativistic electron energies has the form \(f(\varepsilon)
\propto (\varepsilon - \varepsilon_c)^{-1}\) with the condition of 
the bifurcation transition between order and chaos. For relativistic 
electrons the
energy distribution is described  by a power law with the
spectral  index 1/2.

\noindent The characteristic time required for a spectrum to establish itself
hinges on the magnitude of the wave field as \(b^{-2}\) and \(b^{-0.8}\)
respectively, in the non-relativistic and ultra-relativistic range of
electron energy. Mixing over the energy spectrum is more effective
in the non-relativistic particle energy range, thus the heating rate
of relativistic electrons rapidly decreases with its energy (as
\(E^{-2}\)). The distribution of electrons with energies from near
10~keV up to 100~keV in \((v_z, v_t)\) phase space has the
shape of the pancake distribution.

 Stochasticity, the
pitch-angle distribution for non-relativistic electrons appears
sharply peaked near  \(\pi / 2.\)   The pitch-angle
distribution of relativistic electrons however, is described by the
U-like function; the intensity of the pitch-angle spectrum is
dominant along a direction parallel to the external magnetic field.

\noindent The pitch-angle distribution function depends on the
particle energy as \((\varepsilon - \varepsilon_c)^{1/2}\) and is
proportional to  \(\varepsilon^2\) in the range of non-relativistic
and relativistic energies, respectively.

\noindent The high-energy electrons with chaotic motion are equally
likely to scatter in the direction with or against the wave, and
there exists a certain probability that some particles will  be
scattered into a loss cone and will precipitate into the polar regions
along magnetic field lines. The electron transverse (radial)
drift allows precipitation of some electrons into the equatorial region
in addition to the loss cone.

\newpage
\begin{center}
\textbf{REFERENCES}
\end{center}


\noindent Arnold V.,and  Avez A. (1968) Ergodic Problems of
Classical Mechanics.
 Benjamin, N.I.\\
 Chernikov A.A., Tel T., Vattay G., Zaslavsky G.M.(1989)
Phys. Rev. A40, 4072.\\

Karimbadi H. and Angelopoulos V.  (1989) Phys. Rev. A62, 2342.\\
 Klimov O.V., and Tel'nikhin A.A. (1995) Plasma Phys. Rep. 21,
  356.\\
 Kosniowski C. (1980) A First Course in Algebraic Topology. Cambr.
 Univ. Press,
 London.\\
 Krotov A.S., and  Tel'nikhin A.A. (1998) Plasma Phys. Rep. 24, 767.\\
 Li et al. (1997) J. Geophys. Res. 102, 14123.\\
Lichtenberg  A. and Lieberman M. (1983) Regular and Stochastic
Motion. Springer, N.I.\\ Nagornykh E., and Tel'nikhin A.(2002)
Proc.  2 th Int. Congress on Plasma  Phys.: Sydney, Australia (
ed.  J.S. Falkoner, R.L. Dewar, and J.Khachan).
Amer. Inst. of Phys., p. 824.\\
Nagano I. et al. (1996) J. Geomagn. Geoelectr. 48, 299.\\
Parrot M., Gaye C.A. (1994) J. Geophys.Res. 23, 2463.\\
    Schuster H. (1984) Deterministic Chaos.  Physik Verlag. Weinheim. \\
  Schutz B.(1982) Geometrical methods of Mathematical Physics.
Cambr. Univ. Press,  Cambrige, London.\\
Summers D. et al. (1998) J. Geophys.Res. 103, 20487.\\
Wykes W.J., Chapman S.C., Rowlands G. (2001) Planet.Space Sci.,
49, 395. \\
 Zaslavsky G.M. (1985)
Chaos in Dynamical Systems.   Harwood, N.Y.\\
Vedenov A.A., Velikhov E.B. and Sagdeev R.Z. (1962) Nucl. Fusion,
Suppl.
Pt. 2, 465. \\
Drumond W.E. and Pines D. (1982) Nucl. Fusion, Suppl. Pt. 3, 1049.\\
Dupree T.H. (1966) Phys.Fluods, 9, 1773.\\
Weinstock (1969) Phys.Fluods, 12, 1045.\\
Kennel C.F. and Petchek H.E. (1966) Geophys.Res. 71,1.\\
Li et al. (1997) J. Geophys. Res. 102, 14123.\\
Walker A.D.M. (1993) Plasma Waves in the Magnetosphere. Springer,
N.Y. \\
Baker D.N. et al. (1986). J. Geophys.Res. 91, 4265.\\
Horne R.B. and Thorne R.M. (1998) J. Geophys.Res. 99, 17259.\\
Horne R.B, Thorne R.M., N.P. Meredith, and R.R. Anderson (2003) J.
Geophys.Res. 108, 1290. \\
Goldberg R.A., S.A. Curtis and J.R. Barcus (1987) J. Geophys.Res.,
92,  2505. \\
Faith J., S.Kuo and J. Huang (1997) J. Geophys.Res. 102, 2233.\\
Matsoukis K.S., S.Chapman and G. Rowlands (1998) Geophys.Res. Lett.,
25, 265 \\
\end{document}
