1 Introduction
Most astrophysical plasmas are collision free even on very large scales due to the extremely low particle densities. The motion of charged particles is, therefore, not governed by collisions, but by interactions with the ubiquitous turbulent magnetic fields (Bell Reference Bell1978; Schlickeiser Reference Schlickeiser1989). It is a common assumption that this interaction can be described in terms of wave–particle interaction. This process is relevant in astro- and space physics, particle acceleration, plasma heating and the physics of fusion devices. From the interaction of waves and particles all relevant transport parameters (i.e. Fokker–Planck coefficients, mean free path) can be calculated in a typically very tedious procedure.
As wave–particle interaction is a fundamental process, the subject has been frequently revisited during the past decades, oftentimes focusing on relatively simple toy models to demonstrate basic features of more complex physical processes. While this is generally not the case in real world plasmas, models may consider only a single or a limited number of individual plasma waves in order to investigate particle transport. For the simple case of a single, circularly polarized electromagnetic wave propagating along a background magnetic field, the equations of motion of a charged particle are integrable and exact solutions can be found (Roberts & Buchsbaum Reference Roberts and Buchsbaum1964; Arnold Reference Arnold1978; Kong & Liu Reference Kong and Liu2007). Furthermore it can be shown that particles can be trapped in potential wells created by the fields of the wave (Sudan & Ott Reference Sudan and Ott1971; Schreiner, Kilian & Spanier Reference Schreiner, Kilian and Spanier2017a ).
However, this is only a special case of a much wider field of problems. Obliquely propagating waves or waves propagating perpendicular to the background field, for instance, introduce stochastic behaviour of the particles above a certain threshold amplitude of the wave (Palmadesso Reference Palmadesso1972; Smith & Kaufman Reference Smith and Kaufman1978; Varvoglis Reference Varvoglis1984). Chaotic particle orbits may also be encountered as soon as a second wave comes into play. Formerly trapped particles may be de-trapped by a second wave (Sakai & Kamimura Reference Sakai and Kamimura1972), and the equations of motion are in general no longer integrable. Depending on the phase space positions of the wave–particle resonances and the amplitudes of the waves, stochastic behaviour may occur in limited regions of phase space (Murakami, Nomura & Momota Reference Murakami, Nomura and Momota1982; Tran Reference Tran1982; Balakirev et al. Reference Balakirev, Buts, Tolstoluzhkii and Turkin1989; Shklyar & Zimbardo Reference Shklyar and Zimbardo2014). Besides these regions of chaotic particle motion individual islands exist in phase space, where different types of closed orbits can be found (Lehmann & Spatschek Reference Lehmann and Spatschek2010).
One might expect that at least a single, linearly polarized wave propagating along a static background field would present a simple enough problem to find exact solutions to the equations of motion of a charged particle. This is, however, not the case (Bourdier & Gond Reference Bourdier and Gond2001; Bourdier & Patin Reference Bourdier and Patin2005; Bourdier & Drouin Reference Bourdier and Drouin2009): while the equations can be integrated if the linearly polarized wave propagates in an unmagnetized medium, this no longer holds for the magnetized case.
To further add to the discussion of wave–particle interaction and whether the equations of motion of a charged particle are integrable in a specific set-up, we revisit the case of a magnetostatic, monochromatic, circularly polarized wave in a magnetized plasma. We present a new approach to describing the motion of particles in the magnetic field of the wave and derive a solution to the equations of motion which, to our best knowledge, has not been considered elsewhere. The calculations of our dynamical system are shown in detail and are easy to follow. We then make use of numerical simulations to illustrate and expand our results.
Our work is differing from previous attempts, such as Qian (Reference Qian2000) since this work focuses on the integrability of the system and not the resonances of particles and waves. With a potential application in space physics in mind, we choose the magnetostatic case instead of a scenario including electromagnetic waves (such as Balakirev et al.
Reference Balakirev, Buts, Tolstoluzhkii and Turkin1989; Kong & Liu Reference Kong and Liu2007; Bourdier & Drouin Reference Bourdier and Drouin2009; Essén & Stén Reference Essén and Stén2015). In space plasmas, particle transport is dominated by the interaction with Alfvén waves, which are often assumed to be magnetostatic in analytic models, due to their low frequency and constant velocity
$v_{\text{A}}\ll c$
, which is considerably smaller than the speed of light
$c$
.
The article is organized as follows: we present the basic idea of our model for particle transport in magnetostatic slab fluctuations in § 2. A novel set of dynamic variables is motivated in § 3, where the equations of motion are also presented. We then discuss the motion of particles in the magnetic field of a circularly polarized wave in § 4. With our new approach an exact solution to the equations of motion can be derived, which holds without any further assumptions. We present phase space trajectories which show that particles move on closed orbits and calculate equilibrium positions. In § 5 we then address the case of a linearly polarized wave. Using the same mathematical approach, we demonstrate that no analytical solution can be found and that no equilibrium points exist in this case. Numerical simulations of test particles propagating in a prescribed magnetic field set-up are used to illustrate our previous findings in § 6. First we validate our simulations by reproducing the analytical results derived for the case of a monochromatic, circularly polarized wave. As a next step we introduce a second circularly polarized wave with opposite helicity and varying amplitude, resulting in a monochromatic wave with general elliptic polarization. We find chaotic particle motion, but also islands in phase space, where closed orbits seem to be possible. Section 7 concludes the article with a brief summary and discussion.
2 Model
We consider a transverse perturbation
$\unicode[STIX]{x1D6FF}B$
on top of a constant (
$z$
-directed) magnetic field
$B_{0}$
:
$$\begin{eqnarray}\boldsymbol{B}=\unicode[STIX]{x1D6FF}B_{x}(z)\boldsymbol{e}_{x}+\unicode[STIX]{x1D6FF}B_{y}(z)\boldsymbol{e}_{y}+B_{0}\boldsymbol{e}_{z},\end{eqnarray}$$
where
$\boldsymbol{e}_{x}$
,
$\boldsymbol{e}_{y}$
, and
$\boldsymbol{e}_{z}$
are the unit vectors in
$x$
-,
$y$
-, and
$z$
-direction. This can be obtained from the vector potential
$$\begin{eqnarray}\boldsymbol{A}(x,z)=\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{y}(z^{\prime })\,\text{d}z^{\prime }\,\boldsymbol{e}_{x}+\left(B_{0}x-\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{x}(z^{\prime })\,\text{d}z^{\prime }\right)\boldsymbol{e}_{y},\end{eqnarray}$$
since
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{A} & = & \displaystyle \left|\begin{array}{@{}ccc@{}}\boldsymbol{e}_{x} & \boldsymbol{e}_{y} & \boldsymbol{e}_{z}\\ \unicode[STIX]{x2202}_{x} & \unicode[STIX]{x2202}_{y} & \unicode[STIX]{x2202}_{z}\\ A_{x} & A_{y} & 0\end{array}\right|=-\boldsymbol{e}_{\boldsymbol{x}}\frac{\unicode[STIX]{x2202}A_{y}}{\unicode[STIX]{x2202}z}+\boldsymbol{e}_{y}\frac{\unicode[STIX]{x2202}A_{x}}{\unicode[STIX]{x2202}z}+\boldsymbol{e}_{z}\frac{\unicode[STIX]{x2202}A_{y}}{\unicode[STIX]{x2202}x}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D6FF}B_{x}(z)\boldsymbol{e}_{x}+\unicode[STIX]{x1D6FF}B_{y}(z)\boldsymbol{e}_{y}+B_{0}\boldsymbol{e}_{z}.\end{eqnarray}$$
It should be noted that the choice of the vector potential is not unique. Alternatively the potential
$\tilde{\boldsymbol{A}}$
can be used:
$$\begin{eqnarray}\tilde{\boldsymbol{A}}(y,z)=\left(-B_{0}y+\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{y}(z^{\prime })\,\text{d}z^{\prime }\right)\boldsymbol{e}_{x}-\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{x}(z^{\prime })\,\text{d}z^{\prime }\,\boldsymbol{e}_{y}.\end{eqnarray}$$
The Hamiltonian for a particle with mass
$m$
and charge
$q$
in a magnetic vector potential can be written as
$$\begin{eqnarray}H(\boldsymbol{r},\boldsymbol{P})=c\sqrt{(\boldsymbol{P}-q\boldsymbol{A})^{2}+m^{2}c^{2}},\end{eqnarray}$$
with the speed of light
$c$
and the generalized momentum
$\boldsymbol{P}=\unicode[STIX]{x1D6FE}m\,\dot{\boldsymbol{x}}+q\,\boldsymbol{A}$
. For
$\boldsymbol{A}$
the Hamiltonian does not depend on
$y$
while for
$\tilde{\boldsymbol{A}}$
the Hamiltonian does not depend on
$x$
. The Hamiltonian may be written as
$$\begin{eqnarray}H(x,z,\boldsymbol{P})=c\sqrt{(P_{x}-qA_{x})^{2}+(P_{y}-qA_{y})^{2}+P_{z}^{2}+m^{2}c^{2}}.\end{eqnarray}$$
When using the vector potential
$\boldsymbol{A}$
the coordinate
$y$
becomes cyclic. Therefore
$$\begin{eqnarray}P_{y}=p_{y}+qA_{y}=p_{y}+q\left(B_{0}x-\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{x}(z^{\prime })\,\text{d}z^{\prime }\right)\end{eqnarray}$$
is a constant of motion. The Hamiltonian
$H=\unicode[STIX]{x1D6FE}mc^{2}$
does not depend on time explicitly and is thus a constant of motion. Subsequently, the particle momentum
$p$
is also a constant of motion.
Next we define the coordinate
$\boldsymbol{x}_{\boldsymbol{F}}=(x_{F},y_{F},0)$
of the base of the field line which connects the
$x$
–
$y$
-plane to the coordinate
$(x,y,z)$
. To obtain
$x_{F}$
we rewrite
$P_{y}$
using the field line equations. The equation for a field line in the
$x$
–
$z$
-plane is
$$\begin{eqnarray}\frac{\text{d}x}{B_{x}}=\frac{\text{d}z}{B_{0}}\Rightarrow x=x_{F}+\int _{0}^{z}\frac{\unicode[STIX]{x1D6FF}B_{x}(z^{\prime })}{B_{0}}\,\text{d}z^{\prime }.\end{eqnarray}$$
We rephrase the above expression to obtain
$$\begin{eqnarray}q\left(B_{0}x-\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{x}(z^{\prime })\,\text{d}z^{\prime }\right)=qB_{0}x_{F}.\end{eqnarray}$$
Note that the left-hand side is equal to
$q\,A_{y}$
, as defined in (2.2). Thus, we find the relation
$$\begin{eqnarray}P_{y}=p_{y}+qB_{0}x_{F}.\end{eqnarray}$$
For symmetry reasons, we can assume that there must be a constant of motion related to
$P_{x}$
as well, which can be used to obtain
$y_{F}$
:
$$\begin{eqnarray}\frac{\text{d}y}{B_{y}}=\frac{\text{d}z}{B_{0}}\Rightarrow y=y_{F}+\int _{0}^{z}\frac{\unicode[STIX]{x1D6FF}B_{y}(z^{\prime })}{B_{0}}\,\text{d}z^{\prime }.\end{eqnarray}$$
Using the magnetic vector potential
$\tilde{\boldsymbol{A}}$
,
$x$
becomes the cyclic variable, making
$$\begin{eqnarray}\tilde{P_{x}}=p_{x}+q\tilde{A}_{x}=p_{x}+q\left(-B_{0}y+\int _{0}^{z}\unicode[STIX]{x1D6FF}B_{y}(z^{\prime })\,\text{d}z^{\prime }\right)\end{eqnarray}$$
a constant of motion. Relating
$\tilde{P}_{x}$
to
$P_{x}$
we can then write
$$\begin{eqnarray}\displaystyle \tilde{P}_{x} & = & \displaystyle p_{x}+q\tilde{A}_{x}=P_{x}+q(\tilde{A}_{x}-A_{x})=P_{x}-qB_{0}y\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle p_{x}-qB_{0}y_{F}.\end{eqnarray}$$
The Hamiltonian with the vector potential
$\boldsymbol{A}$
then has three constants of motion:
$$\begin{eqnarray}H(x,z,\boldsymbol{P})=c\sqrt{(P_{x}-qA_{x})^{2}+(P_{y}-qA_{y})^{2}+P_{z}^{2}+m^{2}c^{2}};\,P_{y};\,\tilde{P}_{x}=P_{x}-qB_{0}y.\end{eqnarray}$$
Clearly, the latter two do not commute with each other and, thus, the system in general is not integrable. Thus, there is a prospect of chaotic orbits.
In our model we will assume that all particles are injected at
$(x_{0},y_{0},z_{0})=(0,0,0)$
at time
$t=0$
. Since we are discussing a pure slab model this does not change the generality. Their initial velocity is assumed to be
$(v_{0x},v_{0y},v_{0z})$
. Since
$\boldsymbol{A}$
as well as
$\tilde{\boldsymbol{A}}$
are 0 at the origin, we can conclude that
$$\begin{eqnarray}\displaystyle P_{y}(t=0)=p_{y}(t=0) & = & \displaystyle \unicode[STIX]{x1D6FE}mv_{0y},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tilde{P}_{x}(t=0)=p_{x}(t=0) & = & \displaystyle \unicode[STIX]{x1D6FE}mv_{0x}.\end{eqnarray}$$
As discussed earlier, both quantities are constants of motion and therefore conserved at later times.
Before solving the problem of the interaction with magnetic irregularities
$\unicode[STIX]{x1D739}\boldsymbol{B}$
, we will discuss the solution of the unperturbed system (
$\unicode[STIX]{x1D6FF}B=0$
). Using the original Hamiltonian (2.7) with the vector potential (2.2) it is easy to see that the solution is
$$\begin{eqnarray}\displaystyle x(t) & = & \displaystyle \frac{1}{B_{0}q}\left(\vphantom{-}p_{0x}\sin \left(\frac{qB_{0}}{\unicode[STIX]{x1D6FE}m}t\right)+p_{0y}\left(1-\cos \left(\frac{qB_{0}}{\unicode[STIX]{x1D6FE}m}t\right)\right)\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle y(t) & = & \displaystyle \frac{1}{B_{0}q}\left(-p_{0x}\left(1-\cos \left(\frac{qB_{0}}{\unicode[STIX]{x1D6FE}m}t\right)\right)+p_{0y}\sin \left(\frac{qB_{0}}{\unicode[STIX]{x1D6FE}m}t\right)\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle z(t) & = & \displaystyle \frac{p_{0z}t}{\unicode[STIX]{x1D6FE}m}.\end{eqnarray}$$
The straightforward solution shows that particles gyrate around the gyrocentre at
$(x,y)=(p_{0y}/(B_{0}q),-p_{0x}/(B_{0}q))$
. It should be noted that the particles’ gyrocentres do not leave their original magnetic field line. However, the field line coordinate
$(x_{F},y_{F})$
does not denote the field line of the particle’s guiding centre, but the field line currently connecting the
$z=0$
plane with the actual particle.
3 Equations of motion
From the Hamilton equations we see that
$$\begin{eqnarray}\displaystyle {\dot{x}}=v_{x} & = & \displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}P_{x}}=\frac{P_{x}-qA_{x}}{\unicode[STIX]{x1D6FE}m},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{y}}=v_{y} & = & \displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}P_{y}}=\frac{P_{y}-qA_{y}}{\unicode[STIX]{x1D6FE}m},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}}=v_{z} & = & \displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}P_{z}}=\frac{P_{z}}{\unicode[STIX]{x1D6FE}m}.\end{eqnarray}$$
We will now take three steps to derive a system with less variables:
-
(i) Since the total momentum is conserved, we may express the velocity in
$z$
-direction using the
$x$
- and
$y$
-components of the velocity vector. -
(ii) Instead of using the generalized momentum
$P_{x}$
we will switch to the alternative formulation
$\tilde{P}_{x}$
. -
(iii) The magnetic vector potential will be expressed using the field line coordinates
$(x_{F},y_{F})$
according to (2.11) and (2.15).
The first two steps are trivial, the last step requires some caution: the introduction of the coordinates
$(x_{F},y_{F})$
is first of all a mathematical construct. We have discussed earlier that the particle can change the field line which connects it to the
$z=0$
plane without leaving the field line it is ‘tied to’, i.e. the field line along which the gyrocentre of the particle is travelling.
Applying these changes yields
$$\begin{eqnarray}\displaystyle {\dot{x}}=v_{x} & = & \displaystyle \frac{\tilde{P}_{x}+qB_{0}y_{F}(y,z)}{\unicode[STIX]{x1D6FE}m},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{y}}=v_{y} & = & \displaystyle \frac{P_{y}-qB_{0}x_{F}(x,z)}{\unicode[STIX]{x1D6FE}m},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}}=v_{z} & = & \displaystyle \pm \sqrt{v^{2}-v_{x}^{2}(y,z)-v_{y}^{2}(x,z)}.\end{eqnarray}$$
Using the initial conditions introduced in the model, we find
$$\begin{eqnarray}\displaystyle {\dot{x}} & = & \displaystyle v_{x}=v_{0x}+\unicode[STIX]{x1D6FA}_{0}y_{F},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{y}} & = & \displaystyle v_{y}=v_{0y}-\unicode[STIX]{x1D6FA}_{0}x_{F},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle v_{z}=\pm \sqrt{v^{2}-v_{x}^{2}(y,z)-v_{y}^{2}(x,z)},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FA}_{0}=qB_{0}/\unicode[STIX]{x1D6FE}m$
is the cyclotron frequency of the particle.
Since there is a bijective connection between the particle coordinates
$(x,y,z)$
and the coordinates
$(x_{F},y_{F},z)$
, consisting of the field line coordinates
$(x_{F},y_{F},0)$
and the
$z$
-coordinate of the particle, we may express the equations of motion also in these coordinates:
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle \pm \sqrt{v^{2}-v_{x}^{2}(y,z)-v_{y}^{2}(x,z)}\nonumber\\ \displaystyle & = & \displaystyle \pm \sqrt{v^{2}-(v_{0x}+\unicode[STIX]{x1D6FA}_{0}y_{F})^{2}-(v_{0y}-\unicode[STIX]{x1D6FA}_{0}x_{F})^{2}}.\end{eqnarray}$$
The motion of
$x$
and
$y$
in field line coordinates can be derived from the time derivatives of (2.9) and (2.12):
$$\begin{eqnarray}\displaystyle {\dot{x}} & = & \displaystyle {\dot{x}}_{F}+{\dot{z}}\frac{\unicode[STIX]{x1D6FF}B_{x}}{B_{0}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{y}} & = & \displaystyle {\dot{y}}_{F}+{\dot{z}}\frac{\unicode[STIX]{x1D6FF}B_{y}}{B_{0}}.\end{eqnarray}$$
These may be inserted into (3.7) and (3.8):
$$\begin{eqnarray}\displaystyle {\dot{x}}_{F} & = & \displaystyle v_{0x}+\unicode[STIX]{x1D6FA}_{0}y_{F}-{\dot{z}}\frac{\unicode[STIX]{x1D6FF}B_{x}}{B_{0}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{y}}_{F} & = & \displaystyle v_{0y}-\unicode[STIX]{x1D6FA}_{0}x_{F}-{\dot{z}}\frac{\unicode[STIX]{x1D6FF}B_{y}}{B_{0}}.\end{eqnarray}$$
In the next transformation step we will move to the coordinate set
$(X_{G},Y_{G})$
, which describes the position of the field line at
$z=0$
relative to the initial guiding centre position:
$$\begin{eqnarray}\displaystyle X_{G} & = & \displaystyle x_{F}-\frac{v_{0y}}{\unicode[STIX]{x1D6FA}_{0}}=-\frac{v_{y}}{\unicode[STIX]{x1D6FA}_{0}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle Y_{G} & = & \displaystyle y_{F}+\frac{v_{0x}}{\unicode[STIX]{x1D6FA}_{0}}=\frac{v_{x}}{\unicode[STIX]{x1D6FA}_{0}}.\end{eqnarray}$$
The equations of motion in this coordinate set are then
$$\begin{eqnarray}\displaystyle {\dot{X}}_{G} & = & \displaystyle \phantom{-}\unicode[STIX]{x1D6FA}_{0}Y_{G}-{\dot{z}}\frac{\unicode[STIX]{x1D6FF}B_{x}}{B_{0}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{Y}}_{G} & = & \displaystyle -\unicode[STIX]{x1D6FA}_{0}X_{G}-{\dot{z}}\frac{\unicode[STIX]{x1D6FF}B_{y}}{B_{0}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle \pm \sqrt{v^{2}-\unicode[STIX]{x1D6FA}_{0}^{2}(X_{G}^{2}+Y_{G}^{2})}.\end{eqnarray}$$
The coordinates with index
$G$
describe the circular motion of the base of the field line. It is therefore advisable to use polar coordinates:
$$\begin{eqnarray}\displaystyle X_{G} & = & \displaystyle R_{G}\cos \unicode[STIX]{x1D711}_{G},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle Y_{G} & = & \displaystyle R_{G}\sin \unicode[STIX]{x1D711}_{G},\end{eqnarray}$$
for which the following equations of motion are derived:
$$\begin{eqnarray}\displaystyle {\dot{R}}_{G} & = & \displaystyle {\dot{X}}_{G}\cos \unicode[STIX]{x1D711}_{G}+{\dot{Y}}_{G}\sin \unicode[STIX]{x1D711}_{G},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle R_{G}\dot{\unicode[STIX]{x1D711}}_{G} & = & \displaystyle -{\dot{X}}_{G}\sin \unicode[STIX]{x1D711}_{G}+{\dot{Y}}_{G}\cos \unicode[STIX]{x1D711}_{G}.\end{eqnarray}$$
Using (3.17)–(3.19) we can write:
$$\begin{eqnarray}\displaystyle {\dot{R}}_{G} & = & \displaystyle -{\dot{z}}\left(\frac{\unicode[STIX]{x1D6FF}B_{x}(z)}{B_{0}}\cos \unicode[STIX]{x1D711}_{G}+\frac{\unicode[STIX]{x1D6FF}B_{y}(z)}{B_{0}}\sin \unicode[STIX]{x1D711}_{G}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle R\dot{\unicode[STIX]{x1D711}}_{G} & = & \displaystyle {\dot{z}}\left(\frac{\unicode[STIX]{x1D6FF}B_{x}(z)}{B_{0}}\sin \unicode[STIX]{x1D711}_{G}-\frac{\unicode[STIX]{x1D6FF}B_{y}(z)}{B_{0}}\cos \unicode[STIX]{x1D711}_{G}\right)-\unicode[STIX]{x1D6FA}_{0}R_{G},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle \pm \sqrt{v^{2}-\unicode[STIX]{x1D6FA}_{0}R_{G}^{2}}.\end{eqnarray}$$
We solve the last equation for
$R_{G}$
using
${\dot{z}}=v_{z}$
and obtain
$R_{G}=\unicode[STIX]{x1D6FA}_{0}^{-1}\sqrt{v^{2}-v_{z}^{2}}$
. The time derivative of
$R_{G}$
then reads
$$\begin{eqnarray}{\dot{R}}_{G}=-\frac{\dot{v}_{z}v_{z}}{\unicode[STIX]{x1D6FA}_{0}\sqrt{v^{2}-v_{z}^{2}}},\end{eqnarray}$$
and can be inserted into the equations of motion:
$$\begin{eqnarray}\displaystyle \dot{v}_{z} & = & \displaystyle \unicode[STIX]{x1D6FA}_{0}\sqrt{v^{2}-v_{z}^{2}}\left(\frac{\unicode[STIX]{x1D6FF}B_{x}(z)}{B_{0}}\cos \unicode[STIX]{x1D711}_{G}+\frac{\unicode[STIX]{x1D6FF}B_{y}(z)}{B_{0}}\sin \unicode[STIX]{x1D711}_{G}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D711}}_{G} & = & \displaystyle -\unicode[STIX]{x1D6FA}_{0}+\unicode[STIX]{x1D6FA}_{0}\frac{v_{z}}{\sqrt{v^{2}-v_{z}^{2}}}\left(\frac{\unicode[STIX]{x1D6FF}B_{x}(z)}{B_{0}}\sin \unicode[STIX]{x1D711}_{G}-\frac{\unicode[STIX]{x1D6FF}B_{y}(z)}{B_{0}}\cos \unicode[STIX]{x1D711}_{G}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle v_{z}.\end{eqnarray}$$
Using the definition of the pitch-angle cosine
$\cos \unicode[STIX]{x1D703}=\unicode[STIX]{x1D707}=v_{z}/v$
, we find the final set of equations:
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D707}} & = & \displaystyle \unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\left(\frac{\unicode[STIX]{x1D6FF}B_{x}(z)}{B_{0}}\cos \unicode[STIX]{x1D711}_{G}+\frac{\unicode[STIX]{x1D6FF}B_{y}(z)}{B_{0}}\sin \unicode[STIX]{x1D711}_{G}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D711}}_{G} & = & \displaystyle -\unicode[STIX]{x1D6FA}_{0}+\unicode[STIX]{x1D6FA}_{0}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\left(\frac{\unicode[STIX]{x1D6FF}B_{x}(z)}{B_{0}}\sin \unicode[STIX]{x1D711}_{G}-\frac{\unicode[STIX]{x1D6FF}B_{y}(z)}{B_{0}}\cos \unicode[STIX]{x1D711}_{G}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle v\unicode[STIX]{x1D707}.\end{eqnarray}$$
4 Circularly polarized wave
So far we have only specified that the fluctuations are magnetostatic slab fluctuations. We will now focus on the circularly polarized case
$$\begin{eqnarray}\unicode[STIX]{x1D739}\boldsymbol{B}=B_{0}\unicode[STIX]{x1D716}\left(\cos \unicode[STIX]{x1D719}(z)\,\boldsymbol{e}_{x}-\sin \unicode[STIX]{x1D719}(z)\,\boldsymbol{e}_{y}\right),\end{eqnarray}$$
with the relative amplitude
$\unicode[STIX]{x1D716}=|\unicode[STIX]{x1D739}\boldsymbol{B}|/B_{0}$
of the plasma wave. Here the equations of motion, equations (3.31)–(3.33), reduce to
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D707}} & = & \displaystyle \unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\cos (\unicode[STIX]{x1D711}_{G}+\unicode[STIX]{x1D719}(z)),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D711}}_{G} & = & \displaystyle -\unicode[STIX]{x1D6FA}_{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\sin (\unicode[STIX]{x1D711}_{G}+\unicode[STIX]{x1D719}(z)),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle v\unicode[STIX]{x1D707}.\end{eqnarray}$$
These equations are exact.
The radius of motion of the particle mapped to the
$z=0$
plane (using
$X_{G}$
and
$Y_{G}$
as defined in (3.15) and (3.16)) is
$$\begin{eqnarray}R_{G}=\frac{v}{\unicode[STIX]{x1D6FA}_{0}}\sqrt{1-\unicode[STIX]{x1D707}^{2}},\end{eqnarray}$$
as we have derived from (3.26). This expression is identical to the gyroradius in the unperturbed field
$B_{0}$
.
4.1 Monochromatic wave
In the first step we have described the polarization of the wave, but the phase relation
$\unicode[STIX]{x1D719}(z)$
has not been specified. The simplest case is the monochromatic wave with
$\unicode[STIX]{x1D719}(z)=kz$
. The phase between wave and particle is given by
$$\begin{eqnarray}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D711}_{G}+kz.\end{eqnarray}$$
Using equations (4.4) and (4.6) its time derivative is obtained as
$$\begin{eqnarray}\dot{\unicode[STIX]{x1D713}}=\dot{\unicode[STIX]{x1D711}}_{G}+k\unicode[STIX]{x1D707}v,\end{eqnarray}$$
and substituting
$\dot{\unicode[STIX]{x1D711}}_{G}$
from (4.3) gives us a set of two autonomous equations describing particle motion:
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D707}} & = & \displaystyle \unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\cos \unicode[STIX]{x1D713},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D713}} & = & \displaystyle -\unicode[STIX]{x1D6FA}_{0}+kv\unicode[STIX]{x1D707}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\sin \unicode[STIX]{x1D713}.\end{eqnarray}$$
This already shows that there are no chaotic solutions to the problem. By dividing one of the above equations by the other, a first-order ordinary differential equation for one of the variables (either
$\unicode[STIX]{x1D713}$
or
$\unicode[STIX]{x1D707}$
) as a function of the second variable can be obtained (see e.g. Prelle & Singer Reference Prelle and Singer1983; Teschl Reference Teschl2012, chap. 7). The resulting equation can be solved for a family of curves and one constant of integration, which then fixes the correct path of a particle in phase space. This can be done for all systems which depend on only two variables, while a minimum of three variables is required for chaotic behaviour.
For the case of (4.8) and (4.9) we obtain a family of curves in the
$(\unicode[STIX]{x1D713},\unicode[STIX]{x1D707})$
plane:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D713}}{\text{d}\unicode[STIX]{x1D707}}=\frac{-\unicode[STIX]{x1D6FA}_{0}+kv\unicode[STIX]{x1D707}+\unicode[STIX]{x1D6FA}_{0}\displaystyle \frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\unicode[STIX]{x1D716}\cos \unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{\text{d}\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D716}\cos \unicode[STIX]{x1D713}\frac{\text{d}\unicode[STIX]{x1D713}}{\text{d}\unicode[STIX]{x1D707}}=\frac{-\unicode[STIX]{x1D6FA}_{0}+kv\unicode[STIX]{x1D707}+\unicode[STIX]{x1D6FA}_{0}\displaystyle \frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}. & \displaystyle\end{eqnarray}$$
In order to find a relation between
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D713}$
, which allows us to draw the particles’ paths in
$(\unicode[STIX]{x1D713},\unicode[STIX]{x1D707})$
coordinates, we define a function
$F(\unicode[STIX]{x1D707})=\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}$
. We can show the following relation:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}F}{\text{d}\unicode[STIX]{x1D707}}-\frac{\unicode[STIX]{x1D707}}{1-\unicode[STIX]{x1D707}^{2}}F=\frac{kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0}}{\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D707}}(\sqrt{1-\unicode[STIX]{x1D707}^{2}}F)=\frac{kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0}}{\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}\unicode[STIX]{x1D707}}(\sqrt{1-\unicode[STIX]{x1D707}^{2}}F)=\frac{kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0}}{\unicode[STIX]{x1D6FA}_{0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \sqrt{1-\unicode[STIX]{x1D707}^{2}}F=\frac{(kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0})^{2}-C}{2kv\unicode[STIX]{x1D6FA}_{0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}=F=\frac{(kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0})^{2}-C}{2kv\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}. & \displaystyle\end{eqnarray}$$
The problem has been reduced to one integration and is now identical to the motion along
$\unicode[STIX]{x1D707}$
in a given potential. We use
$\unicode[STIX]{x1D716}\cos \unicode[STIX]{x1D713}=\sqrt{\unicode[STIX]{x1D716}^{2}-\unicode[STIX]{x1D716}^{2}\sin ^{2}\unicode[STIX]{x1D713}}$
and are thus able to reformulate equation (4.8). We insert the expression on the right-hand side of (4.16) and eliminate
$\unicode[STIX]{x1D713}$
from the equation. After separating the variables
$\unicode[STIX]{x1D707}$
and
$t$
we are finally able to write the formal integral:
$$\begin{eqnarray}t=\pm \int _{\unicode[STIX]{x1D707}_{0}}^{\unicode[STIX]{x1D707}(t)}\frac{2\,k\,v\,\text{d}\unicode[STIX]{x1D707}^{\prime }}{\sqrt{4\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FA}_{0}^{2}k^{2}v^{2}(1-{\unicode[STIX]{x1D707}^{\prime }}^{2})-[(kv\unicode[STIX]{x1D707}^{\prime }-\unicode[STIX]{x1D6FA}_{0})^{2}-C]^{2}}}.\end{eqnarray}$$
We rearrange (4.16) to find an expression for the constant of integration,
$C$
:
$$\begin{eqnarray}C=(kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0})^{2}-2\unicode[STIX]{x1D6FA}_{0}kv\sqrt{1-\unicode[STIX]{x1D707}^{2}}\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713},\end{eqnarray}$$
which is a constant of motion. Thus, the value of
$C$
can be calculated at
$t=0$
. The contours for
$C$
give the orbits in the
$(\unicode[STIX]{x1D713},\unicode[STIX]{x1D707})$
space, as shown in figure 1. For the trajectories we may easily find the position of
$\text{d}\unicode[STIX]{x1D713}/\text{d}\unicode[STIX]{x1D707}=0$
along the curves:
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\text{d}C=\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\,\text{d}\unicode[STIX]{x1D707}+\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\text{d}\unicode[STIX]{x1D713}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\frac{\text{d}\unicode[STIX]{x1D713}}{\text{d}\unicode[STIX]{x1D707}}=-\frac{\displaystyle \frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}}{\displaystyle \frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}}\Leftrightarrow 0=\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}=2kv(kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0})+2\unicode[STIX]{x1D6FA}_{0}kv\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\left(\frac{kv\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D6FA}_{0}}-1\right)\sqrt{1-\unicode[STIX]{x1D707}^{2}}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D713}. & \displaystyle\end{eqnarray}$$
This relation describes the curves connecting the turning points of the closed contours, as shown by the dashed lines in figure 1. At the intersections with the dashed lines the particle trajectories (represented by the solid lines) are vertical, i.e.
$\text{d}\unicode[STIX]{x1D713}=0$
.
Figure 1. Phase space trajectories (solid curves) for a circularly polarized wave with
$k=2\unicode[STIX]{x1D6FA}_{0}/v$
and
$\unicode[STIX]{x1D716}=0.3$
. The dashed curves give the track of
$\text{d}\unicode[STIX]{x1D713}=0$
along the trajectories. Four equilibrium points – three stable ones (labelled 1, 2, 3) and one unstable one (4) – can be seen at the intersection of the dashed curves and lines
$\unicode[STIX]{x1D713}=(1/2)\unicode[STIX]{x03C0}$
and
$\unicode[STIX]{x1D713}=(3/2)\unicode[STIX]{x03C0}$
.
So far we have discussed particle orbits in
$\unicode[STIX]{x1D713}$
–
$\unicode[STIX]{x1D707}$
-phase space, where
$\unicode[STIX]{x1D707}$
is the pitch-angle cosine in the unperturbed magnetic field
$\boldsymbol{B}_{0}$
. It should be kept in mind that the trajectories have different characteristic when the actual pitch angle
$\unicode[STIX]{x1D6FC}$
– with
$\cos \unicode[STIX]{x1D6FC}=\boldsymbol{v}\cdot \boldsymbol{B}/(|\boldsymbol{v}|\,|\boldsymbol{B}|)$
– is considered. We discuss the relation of
$\cos \unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D707}$
in appendix A.
4.2 Equilibrium points
The equations of motion (4.8) and (4.9) can be rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D6FA}_{0}\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\cos \unicode[STIX]{x1D713}\equiv \unicode[STIX]{x1D6FA}_{0}f(\unicode[STIX]{x1D707},\unicode[STIX]{x1D713}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D6FA}_{0}\left(-1+\unicode[STIX]{x1D705}\unicode[STIX]{x1D707}+\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\sin \unicode[STIX]{x1D713}\right)\equiv \unicode[STIX]{x1D6FA}_{0}g(\unicode[STIX]{x1D707},\unicode[STIX]{x1D713}), & \displaystyle\end{eqnarray}$$
with
$\unicode[STIX]{x1D705}=kv/\unicode[STIX]{x1D6FA}_{0}$
Footnote
1
. The equilibrium is defined by
$\dot{\unicode[STIX]{x1D707}}=0$
and
$\dot{\unicode[STIX]{x1D713}}=0$
. The first equation yields the requirement
$\cos \unicode[STIX]{x1D713}=0$
and we can conclude
$$\begin{eqnarray}\unicode[STIX]{x1D705}\unicode[STIX]{x1D707}=1-\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D713}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}},\end{eqnarray}$$
i.e.
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}\unicode[STIX]{x1D707}=\left\{\begin{array}{@{}ll@{}}1-\unicode[STIX]{x1D716}\displaystyle \frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\quad & \text{for }\unicode[STIX]{x1D713}=\displaystyle \frac{\unicode[STIX]{x03C0}}{2},\\ 1+\unicode[STIX]{x1D716}\displaystyle \frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\quad & \text{for }\unicode[STIX]{x1D713}=\displaystyle \frac{3\unicode[STIX]{x03C0}}{2}.\end{array}\right. & & \displaystyle\end{eqnarray}$$
The quartic equation
$$\begin{eqnarray}(\unicode[STIX]{x1D705}\unicode[STIX]{x1D707}-1)^{2}(1-\unicode[STIX]{x1D707}^{2})=\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D707}^{2}\end{eqnarray}$$
describes the position of the equilibria. Only solutions
$-1\leqslant \unicode[STIX]{x1D707}\leqslant 1$
are physical. There are either two or four equilibrium points for
$\unicode[STIX]{x1D707}$
depending on the values of
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D716}$
. We plot the roots of (4.26) in figure 2 for
$\unicode[STIX]{x1D716}=0.3$
. The four equilibrium points for
$\unicode[STIX]{x1D705}=2$
are also visible in figure 1.
Figure 2. Equilibrium points for
$\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D705})$
for
$\unicode[STIX]{x1D716}=0.3$
. The solid curve gives the positive root at
$\unicode[STIX]{x1D713}=(1/2)\unicode[STIX]{x03C0}$
. The dashed and dash-dotted curves give the roots at
$\unicode[STIX]{x1D713}=(3/2)\unicode[STIX]{x03C0}$
for
$\unicode[STIX]{x1D707}<0$
and
$\unicode[STIX]{x1D707}>0$
, respectively. The dotted line is the limit
$\unicode[STIX]{x1D707}=1/\unicode[STIX]{x1D705}$
at
$\unicode[STIX]{x1D716}=0$
. The red dashed line marks
$\unicode[STIX]{x1D705}=2$
and the labels at the intersections with the black lines correspond to the equilibrium points in figure 1.
The stability of the equilibrium points is determined by the eigenvalues of
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}=\left(\begin{array}{@{}cc@{}}\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}} & \displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\\ \displaystyle \frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}} & \displaystyle \frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\end{array}\right)=\left(\begin{array}{@{}cc@{}}-\unicode[STIX]{x1D716}\displaystyle \frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\cos \unicode[STIX]{x1D713} & -\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}\\ \unicode[STIX]{x1D705}+\displaystyle \frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}} & \unicode[STIX]{x1D716}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\cos \unicode[STIX]{x1D713}\end{array}\right). & & \displaystyle\end{eqnarray}$$
We employ the condition
$\cos \unicode[STIX]{x1D713}=0$
, which we have obtained from combining equation (4.22) with the requirement
$\dot{\unicode[STIX]{x1D707}}=0$
, to find the eigenvalues
$\unicode[STIX]{x1D706}$
:
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \left|\begin{array}{@{}cc@{}}-\unicode[STIX]{x1D706} & -\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}\\ \unicode[STIX]{x1D705}+\displaystyle \frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}} & -\unicode[STIX]{x1D706}\end{array}\right|\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D706}^{2}+\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}\left(\unicode[STIX]{x1D705}+\frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2} & = & \displaystyle -\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}\left(\unicode[STIX]{x1D705}+\frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle -\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}\left(\frac{\unicode[STIX]{x1D705}\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D707}}+\frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}\right).\end{eqnarray}$$
In the next step we replace
$\unicode[STIX]{x1D705}\unicode[STIX]{x1D707}$
according to (4.24) and then make use of
$\sin ^{2}\unicode[STIX]{x1D713}=1$
, since
$\sin \unicode[STIX]{x1D713}=\pm 1$
:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2} & = & \displaystyle -\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}\left(\frac{1}{\unicode[STIX]{x1D707}}\left(1-\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\right)+\frac{\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle -\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\left[\frac{\sin \unicode[STIX]{x1D713}}{\unicode[STIX]{x1D707}}+\frac{\unicode[STIX]{x1D707}^{2}\unicode[STIX]{x1D716}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}\right].\end{eqnarray}$$
We now have to consider the two cases
$\sin \unicode[STIX]{x1D713}=\pm 1$
. For
$\sin \unicode[STIX]{x1D713}=+1$
we find
$\unicode[STIX]{x1D706}^{2}<0$
(marginally stable point) for all
$\unicode[STIX]{x1D707}>0$
and
$$\begin{eqnarray}\unicode[STIX]{x1D707}<-\frac{1}{\sqrt{1+\unicode[STIX]{x1D716}^{2/3}}}.\end{eqnarray}$$
Since the only root for
$\sin \unicode[STIX]{x1D713}=+1$
has
$\unicode[STIX]{x1D707}>0$
it is marginally stable (see label ‘(1)’ in figures 1 and 2).
For
$\sin \unicode[STIX]{x1D713}=-1$
we obtain
$\unicode[STIX]{x1D706}^{2}<0$
for all
$\unicode[STIX]{x1D707}<0$
and
$$\begin{eqnarray}\unicode[STIX]{x1D707}>\frac{1}{\sqrt{1+\unicode[STIX]{x1D716}^{2/3}}}.\end{eqnarray}$$
We note that
$$\begin{eqnarray}\unicode[STIX]{x1D705}=\frac{1}{\unicode[STIX]{x1D707}}-\unicode[STIX]{x1D716}\frac{\sin \unicode[STIX]{x1D713}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}},\end{eqnarray}$$
so
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}=-\frac{1}{\unicode[STIX]{x1D707}^{2}}-\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}=-\frac{1}{\unicode[STIX]{x1D707}}\left(\frac{1}{\unicode[STIX]{x1D707}}+\frac{\unicode[STIX]{x1D707}^{2}\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}}{(1-\unicode[STIX]{x1D707}^{2})^{3/2}}\right) & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D706}^{2}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D713}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}.\end{eqnarray}$$
This yields
$\unicode[STIX]{x1D706}^{2}=0\Leftrightarrow \unicode[STIX]{x2202}\unicode[STIX]{x1D705}/\unicode[STIX]{x2202}\unicode[STIX]{x1D707}=0$
. Thus, the negative root is always stable (see label ‘(2)’ in figures 1 and 2) and the turning point
$\unicode[STIX]{x2202}|\unicode[STIX]{x1D707}|/\unicode[STIX]{x2202}\unicode[STIX]{x1D705}\rightarrow \infty$
at
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{1}{\sqrt{1+\unicode[STIX]{x1D716}^{2/3}}},\quad \unicode[STIX]{x1D705}=(1+\unicode[STIX]{x1D716}^{2/3})^{3/2}\end{eqnarray}$$
of the equilibrium point marked with the dot-dashed curve in figure 2 denotes the loss of stability, when going toward lower values of
$\unicode[STIX]{x1D707}$
. The point is (marginally) stable at
$$\begin{eqnarray}\unicode[STIX]{x1D707}>\frac{1}{\sqrt{1+\unicode[STIX]{x1D716}^{2/3}}}\end{eqnarray}$$
and unstable (a saddle) below (see labels ‘(3)’ and ‘(4)’, respectively, in figures 1 and 2).
5 Linear polarization
5.1 Analytical approach
For the case of linear polarization we may start by replacing (4.1) with the following polarization:
$$\begin{eqnarray}\unicode[STIX]{x1D739}\boldsymbol{B}=B_{0}\unicode[STIX]{x1D716}\cos \unicode[STIX]{x1D719}(z)\boldsymbol{e}_{x}.\end{eqnarray}$$
According to (3.31)–(3.33), this immediately yields the following equations:
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D707}} & = & \displaystyle \unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\unicode[STIX]{x1D716}\cos \unicode[STIX]{x1D711}_{G}\cos \unicode[STIX]{x1D719}(z),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D711}}_{G} & = & \displaystyle -\unicode[STIX]{x1D6FA}_{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\sin \unicode[STIX]{x1D711}_{G}\cos \unicode[STIX]{x1D719}(z),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}} & = & \displaystyle v\unicode[STIX]{x1D707}.\end{eqnarray}$$
Using trigonometric identities the first and second equation are easily transformed into:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\unicode[STIX]{x1D716}{\textstyle \frac{1}{2}}(\cos (\unicode[STIX]{x1D711}_{G}+\unicode[STIX]{x1D719}(z))+\cos (\unicode[STIX]{x1D711}_{G}-\unicode[STIX]{x1D719}(z))), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D711}}_{G}=-\unicode[STIX]{x1D6FA}_{0}+\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}}{2}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}(\sin (\unicode[STIX]{x1D711}_{G}+\unicode[STIX]{x1D719}(z))+\sin (\unicode[STIX]{x1D711}_{G}-\unicode[STIX]{x1D719}(z))). & \displaystyle\end{eqnarray}$$
These results are as expected and identical to the case of two waves with the same wavenumber and opposite magnetic helicity,
$\unicode[STIX]{x1D70E}=\pm 1$
.
Using the phase relation
$$\begin{eqnarray}\unicode[STIX]{x1D719}(z)=kz,\end{eqnarray}$$
we find the equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\unicode[STIX]{x1D716}{\textstyle \frac{1}{2}}(\cos (\unicode[STIX]{x1D711}_{G}+kz)+\cos (\unicode[STIX]{x1D711}_{G}-kz)), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D711}}_{G}=-\unicode[STIX]{x1D6FA}_{0}+\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}}{2}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}(\sin (\unicode[STIX]{x1D711}_{G}+kz)+\sin (\unicode[STIX]{x1D711}_{G}-kz)), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{z}}=v\unicode[STIX]{x1D707}. & \displaystyle\end{eqnarray}$$
Following the previous calculations for circularly polarized waves, we define:
$$\begin{eqnarray}\unicode[STIX]{x1D713}^{\pm }=\unicode[STIX]{x1D711}_{G}\pm kz.\end{eqnarray}$$
Using these abbreviations we can calculate the following:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\frac{\unicode[STIX]{x1D716}}{2}\left(\cos \unicode[STIX]{x1D713}^{+}+\cos \unicode[STIX]{x1D713}^{-}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D713}}^{+}=-\unicode[STIX]{x1D6FA}_{0}+kv\unicode[STIX]{x1D707}+\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}}{2}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\big(\sin \unicode[STIX]{x1D713}^{+}+\sin \unicode[STIX]{x1D713}^{-}\big), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D713}}^{-}=-\unicode[STIX]{x1D6FA}_{0}-kv\unicode[STIX]{x1D707}+\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{0}}{2}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\big(\sin \unicode[STIX]{x1D713}^{+}+\sin \unicode[STIX]{x1D713}^{-}\big). & \displaystyle\end{eqnarray}$$
These equations are a special case of the general expression for two waves, where the magnetic amplitudes are
$\unicode[STIX]{x1D716}^{+}=\unicode[STIX]{x1D716}^{-}=\unicode[STIX]{x1D716}/2$
and the wavenumbers
$k^{+}=-k^{-}=k$
.
One is tempted to follow the approach used for circular polarization. There, a function
$F(\unicode[STIX]{x1D707})=\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}$
had been defined to provide a relation between
$\unicode[STIX]{x1D713}$
and
$\unicode[STIX]{x1D707}$
. In the case of the circularly polarized wave this leads to (4.12), which is an equation for
$\unicode[STIX]{x1D713}$
only depending on
$\unicode[STIX]{x1D707}$
. Unfortunately, both
$\unicode[STIX]{x1D713}^{\pm }$
depend on the other
$\unicode[STIX]{x1D713}$
. This makes the definition
$F^{\pm }(\unicode[STIX]{x1D707})=\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}^{\pm }$
impractical.
We have also tried
$F(\unicode[STIX]{x1D707})=\unicode[STIX]{x1D716}(\sin \unicode[STIX]{x1D713}^{+}+\sin \unicode[STIX]{x1D713}^{-})$
, but we find that:
$$\begin{eqnarray}\displaystyle \frac{\text{d}F}{\text{d}\unicode[STIX]{x1D707}} & = & \displaystyle \unicode[STIX]{x1D716}\left(\frac{\text{d}\unicode[STIX]{x1D713}^{+}}{\text{d}\unicode[STIX]{x1D707}}\cos \unicode[STIX]{x1D713}^{+}+\frac{\text{d}\unicode[STIX]{x1D713}^{-}}{\text{d}\unicode[STIX]{x1D707}}\cos \unicode[STIX]{x1D713}^{-}\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{-2}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}+\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}}{1-\unicode[STIX]{x1D707}^{2}}(\sin \unicode[STIX]{x1D713}^{+}+\sin \unicode[STIX]{x1D713}^{-})+\frac{2kv\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\frac{\cos \unicode[STIX]{x1D713}^{+}-\cos \unicode[STIX]{x1D713}^{-}}{\cos \unicode[STIX]{x1D713}^{+}+\cos \unicode[STIX]{x1D713}^{-}}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Unfortunately, no combination of
$\unicode[STIX]{x1D713}^{\pm }$
and
$\sin \unicode[STIX]{x1D713}^{\pm }$
decouples
$\unicode[STIX]{x1D713}^{\pm }$
from
$\unicode[STIX]{x1D707}$
.
We may – on the other hand – still try to find stationary points:
$\text{d}\unicode[STIX]{x1D713}^{\pm }/\text{d}t=0$
and
$\text{d}\unicode[STIX]{x1D707}/\text{d}t=0$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=0=\unicode[STIX]{x1D6FA}_{0}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\frac{\unicode[STIX]{x1D716}}{2}(\cos \unicode[STIX]{x1D713}^{+}+\cos \unicode[STIX]{x1D713}^{-}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D713}}^{\pm }=0=-\unicode[STIX]{x1D6FA}_{0}\pm kv\unicode[STIX]{x1D707}+\frac{\unicode[STIX]{x1D716}}{2}\frac{\unicode[STIX]{x1D707}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}(\sin \unicode[STIX]{x1D713}^{+}+\sin \unicode[STIX]{x1D713}^{-}). & \displaystyle\end{eqnarray}$$
The two stability conditions for
$\unicode[STIX]{x1D713}^{\pm }$
can never be fulfilled simultaneously as can easily be seen when adding and subtracting them. Thus, there is no stationary point in the
$(\unicode[STIX]{x1D707},\unicode[STIX]{x1D713}^{+},\unicode[STIX]{x1D713}^{-})$
phase space.
This result can also be illustrated by looking at the definition of
$\unicode[STIX]{x1D713}^{\pm }$
in (5.11). For an equilibrium point we define the change
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{+}$
of the relative phase between a particle and the first wave as
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{+}=\unicode[STIX]{x0394}(\unicode[STIX]{x1D711}_{G}+k\,z)=\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{G}+k\,\unicode[STIX]{x0394}z=0.\end{eqnarray}$$
At the same time we require
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{-}=0$
for the relative phase between the particle and the second wave. However, this cannot be achieved due to (5.19), since
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{-}=\unicode[STIX]{x0394}(\unicode[STIX]{x1D711}_{G}-k\,z)=\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{G}-k\,\unicode[STIX]{x0394}z=2\,\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{G}.\end{eqnarray}$$
These considerations show that a particle in resonance with one wave (i.e.
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{+}=0$
) cannot be in resonance with the other wave at the same time. Formally it is possible to satisfy both
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{+}=0$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{-}=0$
, but according to (5.20) this requires
$\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{G}=0$
– a condition which implies that the particle is not gyrating since
$\unicode[STIX]{x1D711}_{G}$
is constant.
5.2 Poincaré sections
We have shown that no stationary points in three-dimensional phase space can exist when a particle is propagating in the fields of a linearly polarized wave. However, we might still be able to derive some more information about the behaviour of the particles. In § 6 we will make use of Poincaré sections to analyse particle motion. A Poincaré section can be defined as a
$(n-1)$
-dimensional cut through an otherwise
$n$
-dimensional phase space, which then allows us to reveal certain characteristics of a dynamical system. Exact mathematical definitions can be found in standard textbooks, such as e.g. Teschl (Reference Teschl2012, chaps 6 and 12).
By means of the Poincaré section a continuous dynamical system is transformed into a discrete system. For the case of particle trajectories this means that the continuous path of a particle in phase space is reduced to a set of discrete points, marking the crossings of the Poincaré section. Thus, periodic or quasi-periodic trajectories in
$n$
-dimensional space can be revealed, since they are represented by either a single point or an accumulation of nearby points in the
$(n-1)$
-dimensional section. Such analysis is a typical method used to study chaotic systems and has been applied to the problem of charged particle transport before (e.g. Murakami et al.
Reference Murakami, Nomura and Momota1982; Bouquet & Bourdier Reference Bouquet and Bourdier1998; Bourdier & Drouin Reference Bourdier and Drouin2009; Lehmann & Spatschek Reference Lehmann and Spatschek2010).
We define our Poincaré sections by the condition
$\unicode[STIX]{x1D713}^{-}=\unicode[STIX]{x1D713}^{+}+n2\unicode[STIX]{x03C0}$
,
$n\in \{0,1,2,\ldots \!\}$
. In order to find the position of a particle trajectory within the Poincaré section, one has to integrate over the particle’s orbit, which requires knowledge of the actual solution of the equation of motion. For the case of periodic solutions the Poincaré section can show equilibria even though there is no stationary point in the full three-dimensional phase space. These periodic (or quasi-periodic) solutions are marked by accumulations of points (or recurring patterns) in the Poincaré section.
6 Magnetostatic particle simulations
We perform numerical test particle simulations to expand our results beyond the case described by the analytical model. As a first step, however, we reproduce the analytical result for the case of a single, circular wave to validate our numerical approach. The magnetostatic test particle code first presented by Schreiner et al. (Reference Schreiner, Kilian and Spanier2017a ) is used. This code initializes a prescribed, static magnetic field and then uses the particle pusher of the particle-in-cell code ACRONYM (Kilian, Burkart & Spanier Reference Kilian, Burkart, Spanier, Nagel, Krner and Resch2012) to propagate the test particles.
6.1 Circular polarization
The test particle code allows us to initialize a background magnetic field
$B_{0}$
along
$z$
-direction together with the static magnetic field
$\unicode[STIX]{x1D739}\boldsymbol{B}(z)$
of a wave with left-handed magnetic helicity, as described by (4.1) for
$k>0$
. Test particles (in this case: protons) are injected as a mono-energetic population at the origin of the coordinate system, i.e.
$(x=0,\,y=0,\,z=0)$
, in accordance with the assumptions made in § 2. The orientation of the velocity vectors of these particles is chosen at random, so that the particles spread out into different directions once the simulation commences.
In certain output intervals the position and velocity of each test particle are saved to disk. This allows us to track the individual particles and to compute their pitch angles
$\unicode[STIX]{x1D703}$
(or the corresponding cosines
$\unicode[STIX]{x1D707}$
) and the phase
$\unicode[STIX]{x1D713}$
according to (4.6). Using (3.15) and (3.16), as well as the definition of the pitch angle in the unperturbed case the two variables can be expressed by:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713} & = & \displaystyle \arccos \left(\frac{-v_{y}}{\sqrt{v_{x}^{2}+v_{y}^{2}}}+kz\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703} & = & \displaystyle \arccos \left(\frac{v_{z}}{v}\right).\end{eqnarray}$$
The required parameters can be easily extracted from the simulation output.
Figure 3. Phase space trajectories of test particles in a magnetostatic simulation. Each coloured line represents the motion of a single particle. The simulation employs a single, circularly polarized wave with
$k=2\unicode[STIX]{x1D6FA}_{0}/v$
and
$\unicode[STIX]{x1D716}=0.3$
. The particles follow the expected trajectories, as predicted by the analytic model (see figure 1). For verification we show the turning points
$\text{d}\unicode[STIX]{x1D713}/\text{d}\unicode[STIX]{x1D707}=0$
(black dashed lines) and positions of the equilibria (labelled with Arabic numbers), as discussed in § 4. The four orange curves (labelled with Roman numbers) are discussed in appendix B.
The simulation uses the same values
$\unicode[STIX]{x1D716}=0.3$
and
$\unicode[STIX]{x1D705}=2$
as previously chosen for figure 1. The three-dimensional simulation box with periodic boundary conditions consists of
$128\times 128\times 2048$
grid cells in the
$x$
-,
$y$
- and
$z$
-directions, which yield a high resolution of the field fluctuations along the
$z$
-direction caused by the plasma wave. Perpendicular to
$\boldsymbol{B}_{0}$
the spatial resolution is not important, since the field configuration does not change along the
$x$
- or
$y$
-axis. The chosen extent of the simulation box allows us to recover the gyromotion of the particles in space in order to obtain some more detailed information about the particles, which can be used for diagnostic purposes. However, since (6.1) and (6.2) only require information about the particles’ velocity vectors and
$z$
-coordinates, a one-dimensional simulation set-up would also be sufficient.
Phase space trajectories of test particles are shown in figure 3. As can be seen, the magnetostatic test particle simulation is able to accurately reproduce the analytical results (see appendix B for a discussion of numerical accuracy). All particles follow closed orbits in phase space and the positions of the equilibrium points are recovered as expected.
6.2 Other polarizations
Next we investigate the behaviour of particles in the presence of two plasma waves. We consider four set-ups (S1 through S4) which each include two circularly polarized waves with opposite magnetic helicity. The first wave has
$k^{+}>0$
and employs
$\unicode[STIX]{x1D716}^{+}=0.3$
and
$\unicode[STIX]{x1D705}^{+}=2$
, as specified in the previous section. We choose
$k^{-}=-k^{+}<0$
and thus
$\unicode[STIX]{x1D705}^{-}=-2$
. The amplitude
$\unicode[STIX]{x1D716}^{-}$
of the second wave is varied in the different set-ups as shown in table 1. Depending on the ratio
$\unicode[STIX]{x1D716}^{-}/\unicode[STIX]{x1D716}^{+}$
the polarization of the superposition of both waves changes from circular to linear.
Table 1. Ratio
$\unicode[STIX]{x1D716}^{-}/\unicode[STIX]{x1D716}^{+}$
of the amplitudes of the forward (
$+$
) and backward (
$-$
) propagating waves and the resulting polarization of their superposition.
For graphical analysis, it is worthwhile to find a suitable representation of the system employing only two variables. We therefore introduce
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}^{-}-\unicode[STIX]{x1D713}^{+}=n2\unicode[STIX]{x03C0},\end{eqnarray}$$
with
$n$
being an integer number including zero. This condition defines a Poincaré section in the three-dimensional parameter space – as described in § 5.2 – which is used as the basis for further analysis. The effective dimensionality of the problem is thus reduced to two and the particle motion can be depicted in
$\unicode[STIX]{x1D713}^{+}$
-
$\unicode[STIX]{x1D707}$
-phase space, similar to the previous analysis of the one wave case. Results from simulations S1 through S4 are shown in figure 4.
Figure 4. Coloured dots represent phase space coordinates of individual particles on the
$\unicode[STIX]{x1D713}^{+}$
–
$\unicode[STIX]{x1D703}$
-plane defined by a Poincaré section. The black dashed lines and black crosses mark the turning points
$\text{d}\unicode[STIX]{x1D713}/\text{d}\unicode[STIX]{x1D707}=0$
and the equilibrium points (labelled 1 through 4 in panel a) for the case of a single circularly polarized wave. The different panels refer to the four set-ups S1, S2, S3 and S4 as specified by the label in each panel. Each label also includes the ratio
$\unicode[STIX]{x1D716}^{-}/\unicode[STIX]{x1D716}^{+}$
of the amplitudes of the forward (
$+$
) and backward (
$-$
) propagating waves, as described in the text and in table 1. The regions in which closed trajectories are possible shrink with increasing amplitude of the second plasma wave (panels a through d). In between these ordered regions chaotic behaviour of the particles can be observed. For illustration, the dots representing the phase space position of two particles (red and blue) are connected in panels (e,f). This gives an impression of particle motion in the Poincaré section and the chaotic orbits these particles are on. Note that lines crossing
$\unicode[STIX]{x1D713}^{+}=0$
(or
$\unicode[STIX]{x1D713}^{+}=2\unicode[STIX]{x03C0}$
) are not shown.
The first set-up, S1, is the same as in the previous section. However, one can formally introduce a second wave with
$k^{-}=-k^{+}$
and zero amplitude (
$\unicode[STIX]{x1D716}^{-}=0$
) to test the analysis using a Poincaré section in the three-dimensional parameter space. The result is shown in figure 4(a) and looks very similar to figure 3. As expected, the data match well with the lines of
$\text{d}\unicode[STIX]{x1D713}/\text{d}\unicode[STIX]{x1D707}=0$
(black dashed lines) and the positions of the equilibrium points (black crosses) which have been derived in § 4. However, the phase space trajectories are now represented by a series of discrete points which do not necessarily form unbroken lines. This is most obvious near
$\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$
, where the particles hardly propagate along the background magnetic field and the condition
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=n2\unicode[STIX]{x03C0}$
, as defined in (6.3), is seldom met. Thus, data are scarce in this region.
Phase space trajectories change with increasing amplitude
$\unicode[STIX]{x1D716}^{-}$
, as can be seen in figure 4(b)–(d). Regions of chaotic particle motion can be distinguished from smaller islands where ordered particle motion is found. We will first discuss the islands, where particles seem to describe closed orbits in the Poincaré section.
As discussed in § 5.2, a stationary point in the Poincaré section indicates that particle motion is periodic in the full three-dimensional phase space. We assume that each island surrounds one such stationary point (see the light blue dot at
$\unicode[STIX]{x1D713}^{+}=(3/2)\unicode[STIX]{x03C0}$
slightly above the lowest black cross in figure 4
b), although we do not have particle data to support the existence of all of them. The ordered structures around these stationary points represent trajectories which are at least quasi-periodic in three-dimensional phase space.
For a wave with elliptical polarization (S2 and S3, see figure 4
b and c) we notice that stationary points (and the corresponding islands) appear close to the positions of the three stable equilibria found in the one wave case (see figures 1 and 4
a, labels 1, 2, 3). An additional fourth island is formed around
$\unicode[STIX]{x1D713}^{+}=(1/2)\unicode[STIX]{x03C0},\,\unicode[STIX]{x1D703}\sim (3/4)\unicode[STIX]{x03C0}$
, which is not related to one of the previously known equilibrium points. As the amplitude
$\unicode[STIX]{x1D716}^{-}$
is increased (from S2 to S3) the islands around
$\unicode[STIX]{x1D713}^{+}=(1/2)\unicode[STIX]{x03C0}$
change their shapes and become smaller, while those around
$\unicode[STIX]{x1D713}^{+}=(3/2)\unicode[STIX]{x03C0}$
keep their triangular appearance, but also shrink slightly. The positions of the stationary points move to larger (smaller)
$|\unicode[STIX]{x1D707}|$
for
$\unicode[STIX]{x1D713}=(1/2)\unicode[STIX]{x03C0}$
(
$\unicode[STIX]{x1D713}=(3/2)\unicode[STIX]{x03C0}$
).
Finally, for the linearly polarized wave with
$\unicode[STIX]{x1D716}^{-}=\unicode[STIX]{x1D716}^{+}$
the shape and position of the ordered islands is symmetric about
$\unicode[STIX]{x1D703}=(1/2)\unicode[STIX]{x03C0}$
(S4, see figure 4
d). Interestingly the islands around
$\unicode[STIX]{x1D713}^{+}=\unicode[STIX]{x03C0}/2$
are split in two, as can be seen by the accumulation of green dots near
$\unicode[STIX]{x1D713}^{+}=(1/2)\unicode[STIX]{x03C0},\,\unicode[STIX]{x1D703}\sim (3/4)\unicode[STIX]{x03C0}$
or the purple dots around
$\unicode[STIX]{x1D713}^{+}=(1/2)\unicode[STIX]{x03C0},\,\unicode[STIX]{x1D703}\sim (1/4)\unicode[STIX]{x03C0}$
, which each correspond to one particle. We assume that these features are associated with two stationary points in the Poincaré section which lie between the two green and purple ellipses, respectively.
In between the different islands of ordered particle orbits it is hard to find a clear structure. Particles seem to travel through phase space in a random motion, which hints at chaotic behaviour. Panels (e) and (f) in figure 4 show the same data as panels (c) and (d), respectively, corresponding to simulations S3 and S4. In each of these two panels, two individual particles are highlighted in colour, and the discrete points denoting the particles’ coordinates in phase space at different times are connected via straight lines to give an impression of particle motion in the Poincaré section. It can be clearly seen that these particles do not enter the ordered islands, but travel through the surrounding space only. At times it is possible that a particle orbits one of the stationary points (see e.g. the red line in panel f). However, these orbits are not stable and the particle returns to a chaotic motion.
Changes of the direction of propagation relative to the background magnetic field
$\boldsymbol{B}_{0}$
are also possible, although not frequent (see figure 4 e, f). This result is interesting from the point of view of quasi-linear theory (e.g. Schlickeiser Reference Schlickeiser1989): the standard magnetostatic approach to wave–particle interaction employs the limit of
$\unicode[STIX]{x1D716}\rightarrow 0$
and finds a resonance gap at
$\unicode[STIX]{x1D707}=0$
, i.e. no particle scattering across that line. Our system, however, shows that a resonance gap can be avoided if finite amplitude (
$\unicode[STIX]{x1D716}>0$
) waves are considered. Even in the relatively simple case of a linearly polarized wave we find chaotic orbits which enable particles to cross
$\unicode[STIX]{x1D707}=0$
.
To conclude the discussion of the numerical results it has to be kept in mind that the trajectories presented in figure 4 do not represent a universal picture. The size and shape of the ordered islands in panels (b–d) depend on the amplitude of the waves. For the case of the linearly polarized wave (figure 4
d) we have carried out additional simulations (not shown here) with varied amplitudes
$\unicode[STIX]{x1D716}^{\pm }$
. For smaller amplitudes the region of stochastic particle transport shrinks, as expected (Murakami et al.
Reference Murakami, Nomura and Momota1982). This is also consistent with the idea that the stochastic regions vanish below a threshold amplitude (Balakirev et al.
Reference Balakirev, Buts, Tolstoluzhkii and Turkin1989). For increasing amplitude, on the other hand, the ordered islands shrink and it seems possible that they can vanish completely. However, for the example discussed here this would require an amplitude which is larger than
$B_{0}$
.
7 Summary and conclusions
We have discussed the problem of particle transport in magnetostatic slab fluctuations. In § 2 we have laid out the general model of a transverse perturbation on top of a static background magnetic field
$\boldsymbol{B}_{0}=B_{0}\boldsymbol{e}_{z}$
. We have derived the Hamiltonian for a charged particle in this general field set-up and have discussed the constants of motion. We have then continued to the equations of motion of the particle in § 3, where we made use of a new set of field line coordinates. These field line coordinates allow us to map the trajectory of the particle to the
$z=0$
plane. Switching to yet another set of variables has allowed us to finally express the dynamics of the system by three differential equations for the pitch-angle cosine
$\unicode[STIX]{x1D707}$
(in the unperturbed field), the phase angle
$\unicode[STIX]{x1D711}_{G}$
in the
$z=0$
plane and the
$z$
-coordinate, equations (3.31)–(3.33).
A detailed analysis for the case of a circularly polarized, monochromatic wave was presented in § 4. We could show that the equations of motion reduce to two autonomous equations, equations (4.8) and (4.9). Thus, we were able to describe the motion of a charged particle in phase space by only two variables, namely
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D713}$
, where
$\unicode[STIX]{x1D713}$
describes the phase between the wave and the particle in the
$z=0$
plane. We have demonstrated that the particles follow closed orbits in phase space and that, therefore, no chaotic behaviour can be found. The discussion was concluded by the analysis of stable equilibrium points in phase space.
In § 5 we have applied our previous approach to the case of a linearly polarized, monochromatic wave. By representing the linearly polarized wave by a superposition of two circularly polarized waves with opposite helicity, we are able to derive a set of three differential equations for
$\unicode[STIX]{x1D707}$
,
$\unicode[STIX]{x1D713}^{+}$
and
$\unicode[STIX]{x1D713}^{-}$
. The latter two variables describe the phase between the particle and the two waves with opposite helicity. However, other than for the case of the circularly polarized wave, it is not possible to separate the variables and solve the equations analytically. We have then briefly discussed the lack of equilibrium points in three-dimensional phase space and the introduction of Poincaré sections for further analysis of the system.
By means of numerical simulations of test particles in a prescribed, static magnetic field, we have revisited particle transport in the case of circularly and linearly polarized, monochromatic waves in § 6. We were able to reproduce our analytical results for the former case, and could show that chaotic behaviour can be found in the latter case. To analyse the linearly polarized wave we have resorted to Poincaré sections through phase space, in order to map particle positions in the three-dimensional phase space to a two-dimensional plane. Besides the chaotic behaviour, we could also show that islands of ordered particle motion exist in the Poincaré section. Stationary points in the Poincaré sections could be found, indicating fully periodic motion in full three-dimensional phase space for these orbits. These do not, however, correspond to equilibrium points of the full system. The calculation of the stationary points in the Poincaré section requires either the knowledge of the trajectory in three-dimensional phase space (which cannot be obtained analytically, as we have shown) or a complicated perturbation ansatz.
While the problem of charged particle transport in magnetostatic slab mode field fluctuations is by no means new and has been treated recurrently and extensively in the literature, our work presents a novel approach. By mapping the particle trajectory onto the
$z=0$
plane we are able to derive a new and elegant set of equations of motion. For the case of a monochromatic, circularly polarized wave we can find an exact analytical solution to the problem. Our solution holds for finite amplitudes of the magnetic perturbations and requires no further approximations.
The derivation of particle trajectories in
$\unicode[STIX]{x1D707}$
-
$\unicode[STIX]{x1D713}$
-phase space clearly shows that particle motion is deterministic. This might seem counterintuitive as one expects resonant scattering processes for waves fulfilling the cyclotron resonance condition. However, a single circularly polarized wave only yields closed orbits in phase space and no effective scattering. Chaotic motion can, therefore, only exist when two or more circularly polarized waves are present.
The numerical test particle simulations presented in § 6 are primarily meant to illustrate our point. However, the set-up for our simulations could also be used as a test case for the validation of numerical codes, such as the Streamline code by Dalena et al. (Reference Dalena, Chuychai, Mace, Greco, Qin and Matthaeus2012). The set-up is easy to implement and the results can be directly compared to the analytical solutions.
For the monochromatic, circularly polarized wave our simulations show the same closed orbits in phase space as our analytical model predicts. However, for the case of a linearly polarized wave, where no analytical solution to the equations of motion can be found, these simulations can contribute to a deeper understanding of the problem. By mapping the phase space trajectories of the test particles to a two-dimensional Poincaré section through the full phase space, we could show that individual islands of stable particle orbits exist. Inside these islands particle motion might be ordered, while it is chaotic elsewhere. Although no equilibrium points in phase space can be derived analytically, the numerical results suggest that there are stationary points in the Poincaré section. Particle trajectories crossing these points are periodic in three-dimensional phase space.
A more involved study of the chaotic basins and regular orbits in the Poincaré sections will be performed in future papers, where we will also extend the study from the polarization to the wave form, i.e. consider two waves with different wavenumber, and to a number of waves greater than two. For small amplitudes we expect that islands of ordered motion in phase space will persist in the Poincaré sections. As seen for the case of two plasma waves, the position of the islands may shift and their size may decrease if the amplitudes of the waves are increased.
It might also be worthwhile to study non-slab geometry with more advanced two-dimensional simulations. However, an analytic approach is not feasible, as the equations of motion include infinite sums of Bessel functions in the case of obliquely propagating waves. In some cases approximate solutions may be found (see e.g. Palmadesso Reference Palmadesso1972) or the equations of motion can be integrated numerically (e.g. Smith & Kaufman Reference Smith and Kaufman1978). It is therefore questionable whether analytic considerations can be of general use, or if numerical test particle simulations might be the better choice straight away.
Acknowledgements
We acknowledge the use of the ACRONYM code and would like to thank the developers (Verein zur Förderung kinetischer Plasmasimulationen e.V.) for their support. Funding: this work is based upon research supported by the National Research Foundation and Department of Science and Technology. Any opinion, findings and conclusions or recommendations expressed in this material are those of the authors and therefore the NRF and DST do not accept any liability in regard thereto. C.S. would like to thank J. Büchner for having him as a guest at the MPS and the Max Planck Society for granting a stipend. R.V. acknowledges the financial support of the Academy of Finland (project 267186).
Appendix A. Change of the actual pitch angle
The pitch-angle cosine
$\unicode[STIX]{x1D707}$
discussed in §§2, 3 and 4 is the angle between the background magnetic field and the particle momentum. This is different from the actual pitch angle, especially when the magnetic field perturbations are no longer negligible. For a more detailed analysis of particle motion or for deriving transport theory it is important to keep in mind that the actual pitch angle
$\unicode[STIX]{x1D6FC}$
behaves differently from expectations based on
$\unicode[STIX]{x1D707}$
.
We find that
$$\begin{eqnarray}\displaystyle \cos \unicode[STIX]{x1D6FC} & = & \displaystyle \frac{\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{B}}{|\boldsymbol{v}||\boldsymbol{B}|}=\frac{v_{x}\unicode[STIX]{x1D6FF}B_{x}+v_{y}\unicode[STIX]{x1D6FF}B_{y}+v_{z}B_{0}}{v\sqrt{\unicode[STIX]{x1D6FF}B_{x}^{2}+\unicode[STIX]{x1D6FF}B_{y}^{2}+B_{0}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{{\dot{x}}\unicode[STIX]{x1D716}\cos kz-{\dot{y}}\unicode[STIX]{x1D716}\sin kz+v\unicode[STIX]{x1D707}}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}.\end{eqnarray}$$
We want to express
$\cos \unicode[STIX]{x1D6FC}$
using the variables
$\unicode[STIX]{x1D713},\unicode[STIX]{x1D707}$
:
$$\begin{eqnarray}\displaystyle \cos \unicode[STIX]{x1D6FC} & = & \displaystyle \frac{({\dot{x}}_{F}+v\unicode[STIX]{x1D707}\unicode[STIX]{x1D716}kz)\unicode[STIX]{x1D716}\cos kz-({\dot{y}}_{F}-v\unicode[STIX]{x1D707}\unicode[STIX]{x1D716}kz)\unicode[STIX]{x1D716}\sin kz+v\unicode[STIX]{x1D707}}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{({\dot{X}}_{G}+v\unicode[STIX]{x1D707}\unicode[STIX]{x1D716}\cos kz)\unicode[STIX]{x1D716}\cos kz-({\dot{Y}}_{G}-v\unicode[STIX]{x1D707}\unicode[STIX]{x1D716}\sin kz)\unicode[STIX]{x1D716}\sin kz+v\unicode[STIX]{x1D707}}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D716}\frac{{\dot{X}}_{G}\cos kz-{\dot{Y}}_{G}\sin kz+v\unicode[STIX]{x1D707}\unicode[STIX]{x1D716}(\cos ^{2}kz+\sin ^{2}kz)}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}+\frac{\unicode[STIX]{x1D707}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D716}\frac{({\dot{R}}_{G}\cos \unicode[STIX]{x1D711}_{G}-R_{G}\dot{\unicode[STIX]{x1D711}}_{G}\sin \unicode[STIX]{x1D711}_{G})\cos kz-({\dot{R}}_{G}\sin \unicode[STIX]{x1D711}_{G}+R_{G}\dot{\unicode[STIX]{x1D711}}_{G}\cos \unicode[STIX]{x1D711}_{G})\sin kz}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D707}\frac{1+\unicode[STIX]{x1D716}^{2}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D716}\frac{{\dot{R}}_{G}(\cos \unicode[STIX]{x1D711}_{G}\cos kz-\sin \unicode[STIX]{x1D711}_{G}\sin kz)-R_{G}\dot{\unicode[STIX]{x1D711}}_{G}(\sin \unicode[STIX]{x1D711}_{G}\cos kz+\cos \unicode[STIX]{x1D711}_{G}\sin kz)}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D707}\frac{1+\unicode[STIX]{x1D716}^{2}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D716}\frac{{\dot{R}}_{G}\cos (\unicode[STIX]{x1D711}_{G}+kz)-R_{G}\dot{\unicode[STIX]{x1D711}}_{G}\sin (\unicode[STIX]{x1D711}_{G}+kz)}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}+\unicode[STIX]{x1D707}\frac{1+\unicode[STIX]{x1D716}^{2}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D716}\frac{-\unicode[STIX]{x1D716}v\unicode[STIX]{x1D707}\cos ^{2}(\unicode[STIX]{x1D711}_{G}+kz)-\unicode[STIX]{x1D716}v\unicode[STIX]{x1D707}\sin ^{2}(\unicode[STIX]{x1D711}_{G}+kz)+\unicode[STIX]{x1D6FA}_{0}R_{G}\sin (\unicode[STIX]{x1D711}_{G}+kz)}{v\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D707}\frac{1+\unicode[STIX]{x1D716}^{2}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D707}+\unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sin \unicode[STIX]{x1D713}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}.\end{eqnarray}$$
Thus, for a constant relative phase, we get a solution of constant
$\cos \unicode[STIX]{x1D6FC}$
. Using
$F=\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D713}$
and the expression from (4.16), we get
$$\begin{eqnarray}\displaystyle \cos \unicode[STIX]{x1D6FC} & = & \displaystyle \frac{\unicode[STIX]{x1D707}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}+\sqrt{\frac{1-\unicode[STIX]{x1D707}^{2}}{1+\unicode[STIX]{x1D716}^{2}}}F\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{2\unicode[STIX]{x1D6FA}_{0}kv\unicode[STIX]{x1D707}+(kv\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FA}_{0})^{2}-C}{2\unicode[STIX]{x1D6FA}_{0}kv\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{k^{2}v^{2}\unicode[STIX]{x1D707}^{2}+\unicode[STIX]{x1D6FA}_{0}^{2}-C}{2\unicode[STIX]{x1D6FA}_{0}kv\sqrt{1+\unicode[STIX]{x1D716}^{2}}}.\end{eqnarray}$$
This result clearly shows the nonlinear relation of
$\cos \unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D707}$
.
While the relation itself is not very complicated, the behaviour of
$\cos \unicode[STIX]{x1D6FC}$
might not seem intuitive. We therefore plot the relation of
$\cos \unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D707}$
in figure 5 according to (A 13). We choose the same values for
$C$
that have been used to plot the trajectories in figure 1. The constant of integration,
$C$
, can only take values in a limited range, defined by the physical parameters
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D716}$
of the system. Thus, there are areas in figure 5 where no physical solution to (A 13) can be obtained (marked in grey).
Figure 5. Relation of the actual pitch-angle cosine
$\cos \unicode[STIX]{x1D6FC}$
in the perturbed magnetic field and
$\unicode[STIX]{x1D707}$
in the unperturbed field according to (A 13). We show
$\cos \unicode[STIX]{x1D6FC}$
as a function of
$\unicode[STIX]{x1D707}$
for
$\unicode[STIX]{x1D705}=2$
,
$\unicode[STIX]{x1D716}=0.3$
and the same values of
$C$
that have been used in figure 1. Physical solutions (black solid lines) can be found in the region between the two dashed boundary lines. The continuations of these solutions (black dotted lines) can also be obtained from (A 13), but are not realized by the particles. Note that
$C$
is limited (see (4.18)) and that no solutions can be obtained for
$C>C_{\text{max}}$
(grey areas at the bottom corners) or
$C<C_{\text{min}}$
(grey area at the top).
Furthermore it should be noted that the trajectory of any given particle does not cover the whole range
$-1<\unicode[STIX]{x1D707}<1$
(see figure 1 to recall the shape of the particle orbits in phase space). We find that the minima and maxima in
$\unicode[STIX]{x1D707}$
occur at
$\unicode[STIX]{x1D713}=(1/2)\unicode[STIX]{x03C0}$
or
$\unicode[STIX]{x1D713}=(3/2)\unicode[STIX]{x03C0}$
. Thus, the available range of the actual pitch-angle cosine is given by
$$\begin{eqnarray}\cos \unicode[STIX]{x1D6FC}=\frac{\unicode[STIX]{x1D707}\pm \unicode[STIX]{x1D716}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}{\sqrt{1+\unicode[STIX]{x1D716}^{2}}}\end{eqnarray}$$
according to (A 10). We plot the two solutions as dashed lines in figure 5. Solid curves between these boundaries show the physical solutions of (A 13) for the given values of
$C$
. The continuations of these curves outside of the boundaries (dotted curves) are not realized by the particles.
Figure 6. Phase space trajectories (solid lines) in
$(\unicode[STIX]{x1D713},\,\unicode[STIX]{x1D6FC})$
coordinates for
$\unicode[STIX]{x1D705}=2$
,
$\unicode[STIX]{x1D716}=0.3$
and the same values of
$C$
that have been used in figure 1. The lines connecting the turning points
$\text{d}\unicode[STIX]{x1D713}/\text{d}\unicode[STIX]{x1D6FC}=0$
are shown as dashed lines, the equilibrium points are marked in red and labelled according to our findings in § 4. Grey areas denote regions in phase space that are not accessible to the particles.
Finally, we transform the particle orbits from figure 1 into the coordinate system
$(\unicode[STIX]{x1D713},\,\unicode[STIX]{x1D6FC})$
and show the result in figure 6. Mapping
$\unicode[STIX]{x1D713}$
-
$\unicode[STIX]{x1D707}$
-phase space to the new coordinate system produces trajectories with slightly different shapes. The four equilibrium points (marked by the red crosses and labels) and the turning points
$\text{d}\unicode[STIX]{x1D713}/\text{d}\unicode[STIX]{x1D6FC}=0$
(dashed lines) can be recovered. Note that certain regions in the new phase space are not accessible to the particles (grey shaded areas).
Appendix B. Numerical accuracy of the test particle simulations
The question of numerical accuracy and potential errors often arises when numerical simulations are employed to model complex physical processes. One way to discuss numerical accuracy is to test the consistency of the simulation itself. For self-consistent methods, such as particle-in-cell (PiC), a variety of scenarios can be analysed to check for numerical errors (see e.g. Melzani et al. Reference Melzani, Winisdoerffer, Walder, Folini, Favre, Krastanov and Messmer2013, for some typical validation methods). The ACRONYM PiC code (Kilian et al. Reference Kilian, Burkart, Spanier, Nagel, Krner and Resch2012), which serves as the basis for the magnetostatic code used in this work, has been thoroughly tested over the years. During the process new test cases have been established, such as the analysis of wave modes (Kilian et al. Reference Kilian, Muñoz, Schreiner and Spanier2017) and their damping characteristics (Schreiner, Kilian & Spanier Reference Schreiner, Kilian and Spanier2017b ), which can be used to check for inconsistencies in the interplay of particle motion and electromagnetic fields.
For other types of codes, such as test particle codes (e.g. Dalena et al.
Reference Dalena, Chuychai, Mace, Greco, Qin and Matthaeus2012) or the hybrid-magneto-hydrodynamic (hybrid-MHD) approach (e.g. Lange et al.
Reference Lange, Spanier, Battarbee, Vainio and Laitinen2013), less rigorous testing methods might be sufficient to ensure correct behaviour of test particles. A simple test could be to simulate gyrating particles in a static magnetic field
$\boldsymbol{B}_{0}=B_{0}\,\boldsymbol{e}_{z}$
. The particles’ orbits can be analysed to obtain their Larmor radii. These in turn can be compared to the exact analytic value to check for numerical errors and to get a measure for numerical accuracy. Lange et al. (Reference Lange, Spanier, Battarbee, Vainio and Laitinen2013) use this method to automatically adjust the numerical time step at the beginning of their test particle MHD simulations. As a rule of thumb, ten time steps per gyration are typically sufficient to reduce numerical errors to negligible levels.
Basically, the simulation described in § 6.1 can also be used for code validation. Particle orbits in phase space can be drawn (see figure 3) and compared to the analytic results (see figure 1). The study can be done on a qualitative level by comparing the shapes of the particle orbits in theory and simulation and by checking the positions of equilibrium points. However, preparing the analytical results might be cumbersome and a qualitative study might not be as convincing as a quantitative error measure.
We therefore suggest a different approach: particle orbits in
$\unicode[STIX]{x1D713}$
–
$\unicode[STIX]{x1D707}$
-space are lines of constant
$C$
, as was derived in (4.18) in § 4. In the analytic model we found that
$C$
is a constant of motion and can thus be calculated at
$t=0$
. This also holds in the numerical simulations.
To check for numerical accuracy it is therefore sufficient to calculate
$C(t_{0})$
for each particle at the beginning of the simulation using its initial
$\unicode[STIX]{x1D713}$
and
$\unicode[STIX]{x1D707}$
. The calculation of
$C$
can be repeated in each consecutive (output) time step
$t$
, then using
$\unicode[STIX]{x1D713}(t)$
and
$\unicode[STIX]{x1D707}(t)$
. Deviations from
$C(t_{0})$
hint at numerical errors or inconsistencies in the code. Note that this test can be performed without any knowledge of the analytical model presented in § 4.
We show results from our simulation in figure 7. For four particles (labelled by Roman numbers, see also figure 1) we track
$C$
over the course of the simulation. As can be seen, the value of
$C$
remains constant for all four particles, except for the occurrence of vertical lines in regular patterns. These vertical lines however do not show numerical errors, but are introduced for diagnostic purposes. Each time a particle crosses
$\unicode[STIX]{x1D713}=0$
the value of
$C$
is deliberately set to
$C(\unicode[STIX]{x1D713}=0)/\unicode[STIX]{x1D6FA}_{0}^{2}=3.5$
(see dotted line in figure 7). This allows us to obtain some more details about particle motion and explains the patterns of the vertical lines: the trajectory of particle (i) is limited in
$\unicode[STIX]{x1D713}$
and crosses
$\unicode[STIX]{x1D713}=0$
twice per orbit with the second crossing being in opposite direction compared to the first one. Particle (ii) is also limited in
$\unicode[STIX]{x1D713}$
, but does not cross
$\unicode[STIX]{x1D713}=0$
. Therefore, no vertical line is obtained in figure 7. Finally, particles (iii) and (iv) traverse the entire range of
$\unicode[STIX]{x1D713}$
, as can be seen in figure 3, and cross
$\unicode[STIX]{x1D713}=0$
once per orbit and each time in the same direction.
Figure 7. The value of the constant of integration,
$C$
, from (4.18) is plotted over time for four different particles (solid, dashed, dash-dotted and dash-dot-dotted lines). The orbits of particles (i) through (iv) are marked in orange and labelled with the same Roman numbers in figure 3. Crossings of
$\unicode[STIX]{x1D713}=0$
are denoted by vertical lines, where the value of
$C$
is deliberately set to
$C(\unicode[STIX]{x1D713}=0)/\unicode[STIX]{x1D6FA}_{0}^{2}=3.5$
(dotted line) for diagnostic purposes.
