Notes on single particle motion II

                Equipo de Propulsión Espacial y Plasmas, Mario Merino, 2021

Variable fields

We will consider only cases with:

L;v/Lωc;E/t/Eωc;B/t/Bωc\ell \ll L; \quad v_\parallel / L \ll \omega_c; \quad |\partial \bm E/\partial t|/E \ll \omega_c; \quad |\partial \bm B/\partial t|/B \ll \omega_c

where LL is the gradient length of changes in E\bm E and B\bm B (LE/EB/BL\sim E/|\nabla \bm E|\sim B/|\nabla \bm B|).
This ensures that the particle does not "see" large changes in E\bm E, B\bm B within an orbit.

We introduce the small parameter ε=/L\varepsilon = \ell / L and rewrite m/qεm^/q^m/q \to \varepsilon \hat m/\hat q, where m^,q^\hat m,\hat q are arbitrary normalization constants:

εdvdt=±q^m^(E+v×B);drdt=v\varepsilon\frac{d \bm v}{dt} = \pm \frac{\hat q}{\hat m}(\bm E + \bm v \times \bm B); \quad \frac{d \bm r}{d t} = \bm v

  • We will study the limit ε0+\varepsilon \to 0^+ and drop the hat symbols for simplicity

Two-scale averaging method

There are two distinct timescales in the problem: τc=1/ωc\tau_c = 1/\omega_c for the fast changes of the gyrophase φ\varphi, and τ\tau for the slow changes of the position of the gyrocenter:

τcττcτε1\tau_c \ll \tau \quad\Rightarrow\quad \frac{\tau_c}{\tau} \sim \varepsilon \ll 1

  • While the gyrophase is a function of time, φ=φ(t)\varphi = \varphi(t), we will write all variables as functions of tt and φ\varphi as if they were two independent variables:

r(t)r(t,φ(t))r(t,φ)\bm r(t) \quad \rightarrow \quad \bm r(t,\varphi(t)) \quad \rightarrow \quad \bm r(t,\varphi)

  • We will compute derivatives as: drdt=rt+φ˙rφ\dfrac{d \bm r}{dt} = \dfrac{\partial \bm r}{\partial t} + \dot\varphi \dfrac{\partial \bm r}{\partial \varphi}, with φ˙ωc=O(1/ε)\dot\varphi \approx \omega_c = O(1/\varepsilon)

Next, we express our solution in the following form:

r(t,φ)=R(t)+ερ(t,φ);v(t,φ)=U(t)+u(t,φ)\bm r(t,\varphi) = \bm R(t) + \varepsilon \bm \rho(t,\varphi); \quad \bm v(t,\varphi) = \bm U(t)+ \bm u(t,\varphi)

  • R\bm R and U=U+U1b\bm U = \bm U_\perp + U_\parallel \bm 1_b describe the gyrocenter, with U=dR/dt\bm U = d\bm R/dt
  • ρ,u\bm \rho,\bm u represent the fast gyromotion:
    • They are perpendicular to B(R,t)\bm B(\bm R,t), i.e., ρ1b=0\bm \rho \cdot \bm 1_b =0, u1b=0\bm u \cdot \bm 1_b =0
    • They are 2π2\pi-periodic in φ\varphi for fixed tt
    • They have zero φ\varphi-average:

<ρ>=12πρdφ=0;<u>=12πudφ=0\left<\bm \rho\right> = \frac{1}{2\pi}\oint \bm \rho d\varphi = \bm 0; \quad \left<\bm u\right> = \frac{1}{2\pi}\oint \bm u d\varphi = \bm 0

The vector basis {1b,1r,1φ}\{\bm 1_b, \bm 1_r, \bm 1_\varphi\} is defined at (R,t)(\bm R,t).

Finally, we expand:

R=R0+εR1+ε2R2+U=U0+εU1+ε2U2+ρ=ρ0+ερ1+ε2ρ2+u=u0+εu1+ε2u2+\bm R = \bm R_0 +\varepsilon\bm R_1 +\varepsilon^2\bm R_2 + \ldots \\ \bm U = \bm U_0 +\varepsilon\bm U_1 +\varepsilon^2\bm U_2 + \ldots \\ \bm \rho = \bm \rho_0 +\varepsilon\bm \rho_1 +\varepsilon^2\bm \rho_2 + \ldots \\ \bm u = \bm u_0 +\varepsilon\bm u_1 +\varepsilon^2\bm u_2 + \ldots

The gyrophase evolves as:

dφdt=1εΩ1+Ω0+εΩ1+\frac{d \varphi}{dt} = \frac{1}{\varepsilon}\Omega_{-1} + \Omega_0 + \varepsilon \Omega_1 + \ldots

Magnetic field gradient:

Consider B(R,t)=B01z\bm B (\bm R,t) = B_0 \bm 1_z. The gradient of B\bm B has different terms:

B=[Bx/xBy/xBz/xBx/yBy/yBz/yBx/zBy/zBz/z]\nabla \bm B = \left[ \begin{array}{cc|c} \partial B_x/\partial x & \partial B_y/\partial x & \partial B_z/\partial x \\ \partial B_x/\partial y & \partial B_y/\partial y & \partial B_z/\partial y \\ \hline \partial B_x/\partial z & \partial B_y/\partial z & \partial B_z/\partial z \\ \end{array} \right]

  • Diagonal terms are the divergence terms. Observe that always B=0\nabla \cdot \bm B = 0
  • The variation of BzB_z with xx and yy are the gradient terms
  • Variations of BxB_x, ByB_y with zz are the curvature terms
  • Terms Bx/y\partial B_x/\partial y and By/x\partial B_y/\partial x are shear terms

We write:

B(r,t)B(R,t)+ε(ρ)B(R,t)+E(r,t)E(R,t)+ε(ρ)E(R,t)+\bm B(\bm r,t)\simeq \bm B(\bm R,t) + \varepsilon(\bm \rho\cdot \nabla) \bm B (\bm R,t) + \ldots \\ \bm E(\bm r,t)\simeq \bm E(\bm R,t) + \varepsilon(\bm \rho\cdot \nabla) \bm E (\bm R,t) + \ldots

and expand the magnetic force as:

±qv×B=±qv×B(R,t)±εqv×(ρ)B(R,t)\pm q\bm v \times \bm B = \pm q\bm v \times \bm B(\bm R,t) \pm \varepsilon q\bm v \times (\bm \rho\cdot\nabla)\bm B(\bm R,t)

  • The second term on the RHS can be thought of as a small "extra force" F\bm F

Furthermore, we require that EE_\parallel be small:

E(R,t)=E(R,t)+εE(R,t)1b\bm E(\bm R,t) = \bm E_\perp(\bm R,t) + \varepsilon E_\parallel(\bm R,t)\bm 1_b

Leading order solution

To order O(1/ε)O(1/\varepsilon), the equation of motion reads:

Ω1du0dφqmu0×B(R,t)=±qmE(R,t)±qmU0×B(R,t)\Omega_{-1}\frac{d \bm u_0}{d\varphi} \mp \frac{q}{m}\bm u_0 \times \bm B(\bm R,t) = \pm\frac{q}{m}\bm E_\perp(\bm R,t) \pm \frac{q}{m}\bm U_{\perp 0} \times \bm B(\bm R,t)

Averaging in φ\varphi:

0=E(R,t)+U0×B(R,t)\bm 0 = \bm E_\perp(\bm R,t) + \bm U_{\perp 0} \times \bm B(\bm R,t)

This yields

U0=E(R,t)×B(R,t)B2(R,t)\bm U_{\perp 0} = \frac{\bm E_\perp(\bm R,t)\times \bm B(\bm R,t)}{B^2(\bm R,t)}

Subtracting this from the equation:

Ω1du0dφqmu0×B(R,t)=0;Ω1dρ0dφ=u0\Omega_{-1}\frac{d \bm u_0}{d\varphi} \mp \frac{q}{m}\bm u_0 \times \bm B(\bm R,t) = \bm 0; \quad \Omega_{-1}\frac{d \bm \rho_0}{d \varphi} = \bm u_0

In the local cylindrical vector basis, integrating and imposing periodicity in φ\varphi:

u0=u01φ;Ω1=qB(R,t)m=ωc(R,t);ρ0=01r\bm u_0 = \mp u_{\perp 0} \bm 1_\varphi; \quad \Omega_{-1} = \mp \frac{qB(\bm R,t)}{m} = \mp\omega_c(\bm R,t); \quad \bm \rho_0 = \ell_0 \bm 1_r

Thus, we recover the solution we already had for uniform, constant fields

Guiding center drifts

To order O(1)O(1) the equation reads:

Ω1u1φqmu1×B(R,t)+dU0dt+u0t+Ω0u0φ=±qmE(R,t)1b±qm(ρ0)E(R,t)±qmU1×B(R,t)±qmU0×(ρ0)B(R,t)±qmu0×(ρ0)B(R,t)\Omega_{-1}\frac{\partial \bm u_1}{\partial \varphi} \mp \frac{q}{m}\bm u_1 \times \bm B(\bm R,t) + \frac{d \bm U_0}{dt} + \frac{\partial \bm u_0}{\partial t} + \Omega_{0}\frac{\partial \bm u_0}{\partial \varphi} \\ = \pm \frac{q}{m}E_\parallel(\bm R,t)\bm 1_b \pm \frac{q}{m}(\bm\rho_0\cdot\nabla)\bm E(\bm R,t) \pm \frac{q}{m}\bm U_{1\perp} \times \bm B(\bm R,t) \\ \pm \frac{q}{m}\bm U_{\perp 0} \times (\bm\rho_0\cdot\nabla)\bm B(\bm R,t) \pm \frac{q}{m}\bm u_0 \times (\bm\rho_0\cdot\nabla)\bm B(\bm R,t)

Averaging in φ\varphi eliminates most of the terms:

dU0dt=±qmE(R,t)1b±qmU1×B(R,t)±qm<u0×(ρ0)B(R,t)>\frac{d \bm U_0}{dt} = \pm \frac{q}{m}E_\parallel(\bm R,t)\bm 1_b \pm \frac{q}{m}\bm U_{1\perp} \times \bm B(\bm R,t) \pm \frac{q}{m}\left<\bm u_0 \times (\bm\rho_0\cdot\nabla)\bm B(\bm R,t)\right>

±qm<u0×(ρ0)B(R,t)>=qmu00<1φ×(1r)B(R,t)>=qmu0012πdφ[1b(1r)B(R,t)1r1r(1r)B(R,t)1b]\pm \frac{q}{m}\left<\bm u_0 \times (\bm\rho_0\cdot\nabla)\bm B(\bm R,t)\right> = - \frac{q}{m} u_{\perp 0} \ell_0\left<\bm 1_\varphi \times (\bm 1_r\cdot\nabla)\bm B(\bm R,t)\right> \\ = \frac{q}{m} u_{\perp 0} \ell_0 \frac{1}{2\pi} \oint d\varphi \left[ \bm 1_b (\bm 1_r\cdot\nabla)\bm B(\bm R,t) \cdot \bm 1_r - \bm 1_r (\bm 1_r\cdot\nabla)\bm B(\bm R,t) \cdot \bm 1_b \right]

To work out this expression use a local Cartesian vector basis {1x,1y,1z}\{\bm 1_x,\bm 1_y, \bm 1_z\}, with

1r=cosφ1x+sinφ1y;1b=1z\bm 1_r = \cos\varphi\bm 1_x + \sin\varphi\bm 1_y; \quad \bm 1_b = \bm 1_z

and use:

B=Bxx+Byy+Bzz=0\nabla\cdot\bm B = \frac{\partial B_x}{\partial x} + \frac{\partial B_y}{\partial y} + \frac{\partial B_z}{\partial z} = 0

After integrating, and noting that B(R,t)=Bz(R,t)1z=B(R,t)1z\bm B(\bm R,t) = B_z(\bm R,t) \bm 1_z = B(\bm R,t) \bm 1_z:

±qm<u0×(ρ0)B(R,t)>=qmu0012B(R,t)=u022B(R,t)B(R,t)\pm \frac{q}{m}\left<\bm u_0 \times (\bm\rho_0\cdot\nabla)\bm B(\bm R,t)\right> = - \frac{q}{m} u_{\perp 0} \ell_0 \frac{1}{2} \nabla B(\bm R,t) = - \frac{u_{\perp 0}^2}{2B(\bm R,t)} \nabla B(\bm R,t)

Introducing the (unsigned) magnetic moment at order O(1)O(1):

μ0(R,t)=mu022B(R,t)\mu_0(\bm R,t) = \frac{m u_{\perp 0}^ 2}{2B(\bm R,t)}

Now, we can rewrite the averaged equation as

mdU0dt=±qE(R,t)1b±qU1×B(R,t)μ0(R,t)B(R,t)m\frac{d \bm U_0}{dt} = \pm q E_\parallel(\bm R,t)\bm 1_b \pm q\bm U_{1\perp} \times \bm B(\bm R,t) -\mu_0(\bm R,t)\nabla B(\bm R,t)

The parallel projection of this equation gives U0U_{\parallel 0}:

mdU0dt=±qE(R,t)μ0(R,t)(1b)B(R,t)m1bdU0dtm\frac{d U_{\parallel 0}}{dt} = \pm q E_\parallel(\bm R,t) -\mu_0(\bm R,t)(\bm 1_b\cdot\nabla) B(\bm R,t) - m\bm 1_b \cdot\frac{d \bm U_{\perp 0}}{dt}

Dropping the (R,t)(\bm R,t) dependency for the sake of notation:

mdU0dt=±qEμ0(1b)Bm1bdU0dtm\frac{d U_{\parallel 0}}{dt} = \pm q E_\parallel -\mu_0(\bm 1_b\cdot\nabla) B - m\bm 1_b \cdot\frac{d \bm U_{\perp 0}}{dt}

  • The second term in the RHS is the magnetic mirror force
  • It pushes particles in the direction of decreasing BB regardless of their sign

The perpendicular projection gives:

U1=±1bωc(R,t)×[dU0dt+μ0(R,t)B(R,t)]\bm U_{1\perp} = \pm \frac{\bm 1_b}{\omega_c(\bm R,t)} \times \left[ \frac{d \bm U_0}{dt} +\mu_0(\bm R,t)\nabla B(\bm R,t) \right]

Dropping the (R,t)(\bm R,t) dependency for the sake of notation:

U1=±1bωc×[dU0dt+μ0mB]\bm U_{1\perp} = \pm \frac{\bm 1_b}{\omega_c} \times \left[ \frac{d \bm U_0}{dt} +\frac{\mu_0}{m}\nabla B \right]

  • The first term is known as the inertial drift
  • The second term is the magnetic drift

Magnetic drift

±μ0mωc1b×B\pm \frac{\mu_0}{m\omega_c} \bm 1_b \times \nabla B

  • It displaces the gyrocenter at right angles with B/1\partial B/\partial \bm 1_\perp
  • The effect is the opposite for positive and negative particles
  • It depends on mm and on u0u_{\perp 0} through μ0\mu_0

Inertial drift

We split it in two parts. The derivative of U01bU_{\parallel 0}\bm 1_b gives:

±1bωc×dU01bdt=±U0ωc1b×d1bdt=±U0ωc1b×[1bt+(U0)1b+U0(1b)1b]\pm \frac{\bm 1_b}{\omega_c} \times \frac{d \bm U_{\parallel 0} \bm 1_b}{dt} = \pm \frac{U_{\parallel 0}}{\omega_c} \bm 1_b \times \frac{d \bm 1_b}{dt} \\ = \pm \frac{U_{\parallel 0}}{\omega_c} \bm 1_b \times \left[ \frac{\partial\bm 1_b}{\partial t} + (\bm U_{\perp 0} \cdot \nabla) \bm 1_b + U_{\parallel 0} (\bm 1_b \cdot \nabla) \bm 1_b \right]

Using the curvature vector κ=(1b)1b\bm \kappa = (\bm 1_b \cdot \nabla) \bm 1_b, the last part is the curvature drift:

±U02ωc1b×κ\pm \frac{U_{\parallel 0}^2}{\omega_c} \bm 1_b \times \bm \kappa

  • It displaces the gyrocenter at right angles with κ\bm \kappa
  • The effect is the opposite for positive and negative particles
  • It depends on mm and on U02U_{\parallel 0}^2

The derivative of U0\bm U_{\perp 0} gives:

±1bωc×dU0dt\pm \frac{\bm 1_b}{\omega_c} \times \frac{d \bm U_{\perp 0}}{dt}

  • This is known as the polarization drift
  • The effect is the opposite for positive and negative particles

Conservation of μ0\mu_0 to O(1)O(1)

The equation of the φ\varphi-averaged kinetic energy of the particle reads:

m2ddt(u02+U02)=±qU0E±qU1E±q<(ρ0)Eu0>\frac{m}{2} \frac{d}{dt} (u_{\perp 0}^2 + U_0^2) = \pm q U_{\parallel 0}E_\parallel \pm q \bm U_{\perp 1}\cdot \bm E_\perp \pm q \left<(\bm \rho_0\cdot \nabla)\bm E \cdot \bm u_0 \right>


mU0dU0dt=±qU0E±qU0(U1×B)μ0U0Bm\bm U_0\cdot\frac{d \bm U_0}{dt} = \pm q U_{\parallel 0}E_\parallel \pm q\bm U_{\perp 0}\cdot (\bm U_{1\perp} \times \bm B) -\mu_0\bm U_0\cdot\nabla B


U0=E×BB2±qU0(U1×B)=±qU1EU_{\perp 0} = \frac{\bm E_\perp\times\bm B}{B^2} \quad\Rightarrow\quad \pm q\bm U_{\perp 0}\cdot (\bm U_{1\perp} \times \bm B)= \pm q \bm U_{\perp 1}\cdot \bm E_\perp


±q<(ρ0)Eu0>=mu022πBdφ(1r)E1φ=μ01b×E=μ0Bt\pm q \left<(\bm \rho_0\cdot \nabla)\bm E \cdot \bm u_0 \right> = \frac{m u_{\perp 0}^2}{2\pi B}\oint d\varphi (\bm 1_r \cdot \nabla)\bm E \cdot \bm 1_\varphi = - \mu_0 \bm 1_b \cdot \nabla \times \bm E = \mu_0 \frac{\partial B}{\partial t}


m2ddtu02μ0Btμ0U0B=0\frac{m}{2} \frac{d}{dt} u_{\perp 0}^2 - \mu_0 \frac{\partial B}{\partial t} -\mu_0 \bm U_0\cdot \nabla B = 0

Since μ0=u02/(2B)\mu_0 = u_{\perp 0}^2/(2B) and dB/dt=B/t+U0BdB/dt = \partial B/\partial t + \bm U_0 \cdot \nabla B in the averaged equation:

ddtμ0=0\frac{d}{dt} \mu_0 = 0

  • μ0\mu_0 is conserved to lowest order and is therefore an adiabatic invariant

See HAZE18 p 29 to learn more about Poincaré invariants and Adiabatic invariants

Magnetic mirror

In stationary fields,

H=12mv2+μB=constH = \frac{1}{2} mv_\parallel^2 + \mu B = \text{const}

μ=mv22B=const\mu = \frac{mv_\perp^2}{2B} = \text{const}

  • Particles trade parallel to perpendicular kinetic energy as they move toward a region of maximum BB
  • If BmaxB_{max} at the maximum, only particles with v2/v2<B/Bmaxv_\perp^2/v^2 < B/B_{max} can escape. This defines a loss cone in phase space

Further reading

  • I have followed mainly HAZE18 (with adaptations)
  • BITT04 gives an ad-hoc derivation for each drift, but easier to understand
  • litt83 is a beautiful derivation based on Hamiltonian dynamics—maybe a topic for another time

I follow mainly HAZE17, with adaptations