-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplasma_writeup.tex
1948 lines (1795 loc) · 103 KB
/
plasma_writeup.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
%% LyX 2.3.4.2 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[11pt,english,nofootinbib, tightenlines, notitlepage]{revtex4-1}
\usepackage{tgpagella}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\setcounter{secnumdepth}{3}
\setlength{\parskip}{\smallskipamount}
\setlength{\parindent}{0pt}
\usepackage{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{cancel}
\usepackage{setspace}
\usepackage{esint}
\usepackage{microtype}
\usepackage[unicode=true,
bookmarks=true,bookmarksnumbered=false,bookmarksopen=false,
breaklinks=false,pdfborder={0 0 1},backref=false,colorlinks=false]
{hyperref}
\hypersetup{pdftitle={Plasma Physics},
pdfauthor={Christopher Albert}}
\begin{document}
\global\long\def\tht{\vartheta}%
\global\long\def\ph{\varphi}%
\global\long\def\balpha{\boldsymbol{\alpha}}%
\global\long\def\btheta{\boldsymbol{\theta}}%
\global\long\def\bJ{\boldsymbol{J}}%
\global\long\def\bGamma{\boldsymbol{\Gamma}}%
\global\long\def\bOmega{\boldsymbol{\Omega}}%
\global\long\def\d{\mathrm{d}}%
\global\long\def\t#1{\text{#1}}%
\global\long\def\m{\text{m}}%
\global\long\def\v#1{\boldsymbol{#1}}%
\global\long\def\t#1{\mathbf{#1}}%
\global\long\def\tT{\boldsymbol{\mathsf{T}}}%
\global\long\def\tI{\boldsymbol{\mathsf{I}}}%
\global\long\def\tP{\boldsymbol{\mathsf{P}}}%
\global\long\def\vv{\boldsymbol{v}}%
\title{Plasma Physics}
\author{Christopher Albert}
\date{\today}
\maketitle
This document is a write-up for the plasma physics course at TU Graz
with the goal of becoming the new lecture notes. While derivations
in the original lecture notes are kept short, they are listed more
extensively here for better understanding. Before going to the actual
problems, we summarize some general concepts that are required for
treating them. Most of those concepts have a much wider application
than plasma physics, so it's a good idea to be familiar with them.
\tableofcontents{}
\section*{Common Concepts}
\subsection*{Concept 1: Tensor calculus}
First we summarize some rules of tensor calculus that are required
for certain tasks within plasma physics. A summary of rules of vector
and tensor calculus is available in the NRL plasma formulary. Some
rules are
BAC-CAB rule:
\begin{equation}
\v a\times(\v b\times\v c)=(\v c\times\v b)\times\v a=\v b(\v a\cdot\v c)-\v c(\v a\cdot\v b)
\end{equation}
The scalar product is not commutative for tensors:
\begin{equation}
\v a\cdot\nabla\v b=(\v a\cdot\nabla)\v b=\v a\cdot(\nabla\v b)\neq(\nabla\v b)\cdot\v a
\end{equation}
Cross-product with a curl:
\begin{equation}
\v a\times(\nabla\times\v b)=(\nabla\v b)\cdot\v a-\v a\cdot\nabla\v b\label{eq:eq3}
\end{equation}
Divergence of a dyad:
\begin{equation}
\nabla\cdot(\v a\v b)=(\v a\cdot\nabla)\v b+(\nabla\cdot\v a)\v b\label{eq:diaddiv}
\end{equation}
Those rules can be deduced from using index notation for the respective
operations. For application in research and engineering it is important
to note that tensor calculus can be formulated in a flexible and coordinate-independent
way within the co-contravariant formalism. More detailed derivations
on this very convenient way of treating curvilinear coordinates can
be found in the first chapters of the book of Callen/d'haeseleer \citep{dhaeseleer1991}.
\subsubsection*{Taylor expansion of vector fields\label{subsec:Taylor-expansion-of}}
The second order expansion of a vector field $\v A(\v r)$ around
a point $\v R$ is given by
\begin{align}
\v A(\v R+\v{\rho}) & =\v A(\v R)+(\v{\rho}\cdot\nabla)\v A(\v R)+(\v{\rho}\v{\rho}:\nabla\nabla))\v A(\v R)+\mathcal{O}(\rho^{3}).
\end{align}
Here $\v{\rho}\v{\rho}$ and $\nabla\nabla$ are two dyades combined
by a scalar product. where differentiation doesn't act on $\v{\rho}$.
In index notation this means
\begin{equation}
A_{i}(\v R+\v{\rho})=A_{i}(\v R)+\rho_{j}\partial_{j}A_{i}(\v R)+(\rho_{k}\rho_{l}\partial_{k}\partial_{l})A_{i}(\v R)+\mathcal{O}(\rho^{3}).
\end{equation}
The convention is that a sum is taken automatically over indexes appearing
twice, and we skip the sum sign. If we consider the dyade
\begin{equation}
(\nabla\v A)_{kl}=\partial_{k}A_{l}.
\end{equation}
We see that the first order term can be written without brackets as
\begin{equation}
(\v{\rho}\cdot\nabla)\v A=\v{\rho}\cdot\nabla\v A.
\end{equation}
For the second order term we combine
\begin{equation}
\rho_{k}\rho_{l}\partial_{k}\partial_{l}=(\rho_{k}\partial_{k})(\rho_{l}\partial_{l}).
\end{equation}
Thus it can also be written as a dyad itself via
\begin{equation}
\v{\rho}\v{\rho}:\nabla\nabla=(\v{\rho}\cdot\nabla)(\v{\rho}\cdot\nabla),
\end{equation}
where derivatives don't act on $\v{\rho}$. As an example we can look
at Cartesian coordinates $(x^{1},x^{2})$ in 2D and a vector component
\begin{align*}
A_{i}(\v R+\v{\rho}) & =A_{i}(\v R)\\
& +\rho^{1}\partial_{1}A_{i}(\v R)+\rho^{2}\partial_{2}A_{i}(\v R)\\
& +(\rho^{1})^{2}\partial_{11}A_{i}(\v R)+2\rho^{1}\rho^{2}\partial_{12}A_{i}(\v R)+(\rho^{2})^{2}\partial_{22}A_{i}(\v R)\\
& +\mathcal{O}(\rho^{3}).
\end{align*}
Here we used the shorthand notation $\partial_{k}\partial_{l}\equiv\partial_{kl}$.
Since $\rho^{1}\rho^{2}\partial_{12}=\rho^{2}\rho^{1}\partial_{21}$
a coefficient $2$ appears in front of this term in the second order
term of the Taylor expansion. Going to higher orders, such factors
will be binomial coefficients in 2D and multinomial coefficients in
$N$-D.
\subsection*{Concept 2: Fluid mechanics}
Let's look at the general form of fluid equations. The derivation
relies on the intuitive picture what a ``flux'' is. A rigorous derivation
of fluid equations can be performed from individual particles via
kinetic theory.
We start with the conservation law for mass density $\rho_{m}(\v r,t)$
known as the continuity equation. In a fluid moving with velocity
$\v v(\v r,t)$, mass is transported by the mass flux
\begin{equation}
\boldsymbol{\Gamma}_{m}(\v r,t)=\rho_{m}(\v r,t)\v v(\v r,t).
\end{equation}
The mass flux is a vector field parallel to $\v v(\v r,t)$ and of
scales with the mass available to move in $\v x$ quantified by $\rho_{m}(\v r,t)$.
Now we consider a volume $\Omega$ fixed to the lab frame.\footnote{This is the Eulerian picture, as opposed to the Lagrangian picture,
where a volume is moved and distorted together with the fluid flow. } We compute the total mass inside $\Omega$ via a volume integral
over $\rho_{m}$. If the system's total mass should be conserved,
a change in mass in the volume $\Omega$ can only be caused by an
in-/outflux of mass via $\v{\Gamma}_{m}$ over its boundary surface
$\partial\Omega$. Mathematically we formulate this as a relation
between a partial time derivative over the total mass and a surface
integral over $\v{\Gamma}_{m}$ and obtain the law of\textbf{}\\
\textbf{Mass conservation (integral)}
\begin{equation}
\frac{\partial}{\partial t}\int_{\Omega}\rho_{m}(\v r,t)\d V+\int_{\partial\Omega}\boldsymbol{\Gamma}_{m}(\v r,t)\cdot\d\v S=0.
\end{equation}
Note that a positive flux $\boldsymbol{\Gamma}_{m}$ points outwards,
such that it needs to be compensated by a reduction of integral mass
via its partial time derivative. Now we can exchange derivative and
integral, as $\Omega$ is fixed during time and use the divergence
theorem to obtain
\begin{equation}
\int_{\Omega}\frac{\partial\rho_{m}(\v r,t)}{\partial t}+\nabla\cdot\boldsymbol{\Gamma}_{m}(\v r,t)\d V=0.
\end{equation}
Since this relation is true independent on the choice of $\Omega$,
it must hold for the integrand. Thus we obtain the law of\textbf{\smallskip{}
}\\
\textbf{Mass conservation (differential) / continuity equation}
\begin{equation}
\frac{\partial\rho_{m}(\v r,t)}{\partial t}+\nabla\cdot\boldsymbol{\Gamma}_{m}(\v r,t)=0.
\end{equation}
or explicitly
\begin{equation}
\frac{\partial\rho_{m}(\v r,t)}{\partial t}+\nabla\cdot(\rho_{m}(\v r,t)\v v(\v r,t))=0.
\end{equation}
For a single species of particles of mass $m$ we can write $\rho_{m}(\v r,t)=m\,n(\v r,t)$
via the number density $n(\v x,t)$ and write the law of\textbf{\smallskip{}
}\\
\textbf{Number density conservation}
\begin{equation}
\frac{\partial n(\v r,t)}{\partial t}+\nabla\cdot(\v v(\v r,t)n(\v r,t))=0.
\end{equation}
The second conserved quantity we are interested in is the momentum
density $\boldsymbol{p}$ in absence of external forces. Since momentum
is a vector, its flux $\boldsymbol{\Gamma}_{\v p}$ must be a rank-2
tensor. The result is still the same as for mass or any other conservation
law, yielding the\textbf{\smallskip{}
}\\
\textbf{Total momentum conservation} \textbf{law (differential)}
\begin{equation}
\frac{\partial\boldsymbol{p}(\v r,t)}{\partial t}+\nabla\cdot\boldsymbol{\Gamma}_{\v p}(\v r,t)=0.
\end{equation}
In case of external forces, they would appear as a \emph{source} on
the right-hand side as force densities $\boldsymbol{f}$. The overall
momentum density consists of \emph{kinematic fluid} $\v p^{\text{fl}}=\rho_{m}\v v$
and \emph{electromagnetic} momentum density $\v p^{\text{EM}}=\boldsymbol{S}/c^{2}$,
where $\boldsymbol{S}$ is the Poynting vector. The momentum flux
$\boldsymbol{\Gamma}_{\v p}$ consists of a fluid part\footnote{Here $\t T=\boldsymbol{v}\v w$ means the outer product $\boldsymbol{v}\otimes\boldsymbol{w}$,
being a rank-2 tensor with Cartesian components $T_{ij}=v_{i}w_{j}$.} $\boldsymbol{\Gamma}_{\v p}^{\text{fl}}=\v v\v p_{\text{fl}}$, but
is also influenced by \emph{thermodynamic} effects, namely the width
of the distribution around the average velocity $\boldsymbol{v}$.
We separate out \emph{electrodynamic} and \emph{thermodynamic} effects
in the total momentum conservation law, such that they appear as forces
on the right-hand side, leading to the\textbf{\smallskip{}
}\\
\textbf{Kinematic fluid momentum law}
\begin{equation}
\frac{\partial\v p^{\mathrm{fl}}(\v r,t)}{\partial t}+\nabla\cdot(\v v(\v r,t)\v p^{\mathrm{fl}}(\v r,t))=-\nabla p(\v r,t)+q\left(\v E(\v r,t)+\v v(\v r,t)\times\v B(\v r,t)\right).
\end{equation}
In Cartesian coordinates the divergence of a rank-2 tensor $\t T$
is a column vector whose components are the divergence of tensor columns,
so
\begin{equation}
\nabla\cdot\t T=\left(\begin{array}{c}
\nabla\cdot\boldsymbol{T}_{\cdot1}\\
\nabla\cdot\boldsymbol{T}_{\cdot2}\\
\nabla\cdot\boldsymbol{T}_{\cdot3}
\end{array}\right)=\left(\begin{array}{c}
\frac{\partial}{\partial x^{1}}T_{11}+\frac{\partial}{\partial x^{2}}T_{21}+\frac{\partial}{\partial x^{3}}T_{31}\\
\frac{\partial}{\partial x^{1}}T_{12}+\frac{\partial}{\partial x^{2}}T_{22}+\frac{\partial}{\partial x^{3}}T_{32}\\
\frac{\partial}{\partial x^{1}}T_{13}+\frac{\partial}{\partial x^{2}}T_{23}+\frac{\partial}{\partial x^{3}}T_{33}
\end{array}\right).
\end{equation}
Now we can further modify this using mass continuity by using a rule
from tensor calculus to write the flux vector $\nabla\cdot\boldsymbol{\Gamma}_{\v p}^{\text{fl}}$
as
\begin{align}
\nabla\cdot\boldsymbol{\Gamma}_{\v p}^{\mathrm{fl}}=\nabla\cdot(\v v\v p^{\mathrm{fl}}) & =\v p^{\mathrm{fl}}(\nabla\cdot\v v)+(\v v\cdot\nabla)\v p^{\mathrm{fl}}.
\end{align}
Here, $\boldsymbol{v}\cdot\nabla$ performs a convective derivative,
i.e.
\begin{equation}
(\v v\cdot\nabla)\boldsymbol{w}=v_{1}\frac{\partial}{\partial x^{1}}\boldsymbol{w}+v_{2}\frac{\partial}{\partial x^{2}}\boldsymbol{w}+v_{3}\frac{\partial}{\partial x^{3}}\boldsymbol{w}.
\end{equation}
Now let's look back and see that \emph{incidentally }(well, not entirely)
the \emph{fluid momentum density} $\v p^{\mathrm{fl}}=\rho_{m}\v v=\boldsymbol{\Gamma}_{m}$
matches the \emph{mass flux} defined before.
\section{Basics}
\subsection*{1.1 Criteria for a plasma}
\subsection*{Problem 1: Work needed to introduce charge density fluctuation}
see lecture notes
\subsection*{Problem 2: Debye shielding - the easy way}
see lecture notes
\subsection*{Problem 3: Debye shielding - the hard way}
see lecture notes
Helmholtz equation
\begin{equation}
\Delta\Phi(\v x)-\frac{1}{\lambda^{2}}\Phi(\v x)=-\frac{Q}{\varepsilon_{0}}\delta(\v x).
\end{equation}
Linear problem and infinite or periodic space: Use spatial Fourier
transform. Be careful about normalization with $1/\sqrt{2\pi}$ (symmetric,
unitary) or $1/(2\pi)$ (non-unitary) on one of either transform or
inverse transform.
1-dimensional Fourier transform of $\delta(x)$:
\begin{align}
\mathcal{F}\delta(x) & =\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\d x\,\delta(x)e^{-ikx}\nonumber \\
& =\frac{e^{-ik0}}{\sqrt{2\pi}}=\frac{1}{\sqrt{2\pi}}.
\end{align}
$N$-dimensional Fourier transform of $\delta(\v x)$:
\begin{align}
\mathcal{F}_{N}\delta(\v x) & =\frac{1}{(2\pi)^{N/2}}\int_{-\infty}^{\infty}\d^{N}x\,\delta(\v x)e^{-i\v k\cdot\v x}\nonumber \\
& =\frac{e^{-i\v k\cdot\v 0}}{(2\pi)^{N/2}}=\frac{1}{(2\pi)^{N/2}}.
\end{align}
So in our 3D equation:
\begin{equation}
-k^{2}\Phi(\v k)-\frac{1}{\lambda^{2}}\Phi(\v k)=-\frac{Q}{(2\pi)^{3/2}\varepsilon_{0}}.
\end{equation}
Solution in $\v k$-space:
\begin{equation}
\Phi(\v k)=\frac{Q}{(2\pi)^{3/2}\varepsilon_{0}(k^{2}+1/\lambda_{D}^{\,2})}.
\end{equation}
Solution in real space:
\begin{align}
\Phi(\v x)=\mathcal{F}_{3}^{\,-1}\Phi(\v k) & =\frac{1}{(2\pi)^{3/2}}\int_{-\infty}^{\infty}\d^{3}k\,\Phi(\v k)\\
& =\frac{1}{(2\pi)^{3/2}}\int_{-\infty}^{\infty}\d^{3}k\,\frac{Q}{(2\pi)^{3/2}\varepsilon_{0}(k^{2}+1/\lambda_{D}^{\,2})}e^{i\v k\cdot\v x}.
\end{align}
How do we represent solve this? Represent scalar product by
\begin{equation}
\v k\cdot\v x=|\v k|\,|\v x|\cos\tht_{(k)}=kr\cos\tht_{kx}
\end{equation}
where $k=|\v k|$, $r=|\v x|$ and $\tht_{kx}$ is the angle between
vectors $\v k$ and $\v x$. For the integral we first rotate the
coordinate system of $\v k$-space such that the $z$-axis coincides
with the vector $\v x$. The advantage of this coordinate frame is
that
\begin{equation}
k_{z}=\v k\cdot\frac{\v x}{|\v x|}.
\end{equation}
Now we switch from $(k_{x},k_{y},k_{z})$ to spherical coordinates
$(k,\tht_{(k)},\ph_{(k)})$ where
\begin{align}
k & =|\v k|,\\
k_{z} & =k\cos\tht_{(k)}.
\end{align}
We observe that the a
\begin{equation}
\v k\cdot\v x=kr\cos\tht_{(k)}=rk_{z}.
\end{equation}
With this trick, we can thus represent the term appearing in the complex
exponent as
\begin{equation}
e^{i\v k\cdot\v x}=e^{irk_{z}}=e^{ikr\cos\tht_{(k)}}.
\end{equation}
Now back to our original problem. We compute
\begin{align}
\Phi(\v x) & =\frac{Q}{(2\pi)^{3}\varepsilon_{0}}\int_{0}^{\infty}\d k\int_{0}^{\pi}\d\tht_{(k)}\int_{0}^{2\pi}\d\ph_{(k)}\,k^{2}\sin\tht\frac{k^{2}}{k^{2}+1/\lambda_{D}^{\,2}}e^{ikr\cos\tht_{(k)}}.
\end{align}
The next trick is to substitute
\begin{align}
u= & \cos\tht_{(k)},\\
\d u= & -\sin\tht_{(k)}\d\tht_{(k)}.
\end{align}
We cancel the negative sign in $\d u$ by swapping the integration
boundaries, which become $(-1,1)$ for $u$ where $\tht_{(k)}=(\pi,0)$.
As there are no dependencies on $\varphi_{(k)}$ integration just
adds a factor of $2\pi$.
\begin{align}
\Phi(\v x) & =\frac{Q}{(2\pi)^{2}\varepsilon_{0}}\int_{0}^{\infty}\d k\frac{k^{2}}{k^{2}+1/\lambda_{D}^{\,2}}\int_{-1}^{1}\d ue^{ikru}\\
& =\frac{Q}{(2\pi)^{2}\varepsilon_{0}}\int_{0}^{\infty}\d k\frac{k^{2}}{k^{2}+1/\lambda_{D}^{\,2}}\frac{1}{ikr}(e^{ikr}-e^{-ikr})
\end{align}
The next trick is to replace the integral over the second term $e^{-ikr}$
by setting $k\rightarrow-k$, and include it by extending the remaining
integral in the range $(-\infty,0)$,
\begin{equation}
\Phi(\v x)=\frac{Q}{ir(2\pi)^{2}\varepsilon_{0}}\int_{-\infty}^{\infty}\d k\frac{k}{k^{2}+1/\lambda_{D}^{\,2}}e^{ikr}.
\end{equation}
We see that the integral is now an inverse 1D Fourier transform of
an expression of the form $k/(k^{2}+a^{2})$ that we can look up to
be
\begin{equation}
\mathcal{F}^{(-1)}\frac{a}{k^{2}+a^{2}}=\sqrt{\frac{\pi}{2}}e^{-a|x|},
\end{equation}
so taking one derivative
\begin{equation}
\mathcal{F}^{(-1)}\frac{ik}{k^{2}+a^{2}}=\sqrt{\frac{\pi}{2}}\frac{1}{a}\frac{\d}{\d x}e^{-a|x|}=-\sqrt{\frac{\pi}{2}}e^{-a|x|}.
\end{equation}
Finally, with $a=1/\lambda_{D}$, and $x=r>0$ we obtain
\begin{align}
\Phi(\v x) & =\frac{Q}{ir(2\pi)^{2}\varepsilon_{0}}\frac{\sqrt{2\pi}}{i}\mathcal{F}^{(-1)}\frac{ik}{k^{2}+1/\lambda_{D}^{\,2}}\nonumber \\
& =\frac{Q}{4\pi\varepsilon_{0}r}e^{-r/\lambda_{D}}.
\end{align}
\begin{itemize}
\item Summary:
\begin{itemize}
\item We computed \emph{potential }of \emph{point charge} inside \emph{electron
fluid} with a certain \emph{temperature }where electrons are allowed
to \emph{move} and thus \emph{shield} the potential.
\item First we used vortex-free \emph{momentum law} (Euler/Navier-Stokes)
to compute \emph{local} \emph{Boltzmann law} for \emph{electron density}
according to $V$, $\Phi$ and $T$.
\item Then we inserted convection-free ($V=0$) case into \emph{Poisson
equation}\textbf{\emph{ }}and \emph{linearized} to obtain linear \emph{Helmholtz
equation}, where \emph{Debye length }$\lambda_{D}$ appeared.
\item We solved the resulting linear equation via a spatial \emph{Fourier
transform} and a few tricks.
\item The result is the usual potential of a point charge weighted by an
\emph{exponential decaying} term on the scale of $\lambda_{D}$.
\item The interpretation is, that outside of $\lambda_{D}$, perturbations
in charge density become \emph{invisible} $\rightarrow$ \emph{Debye
shielding}.
\end{itemize}
\item Remember: Fourier \emph{transform} works only in an infinite domain
$\mathbb{R}^{N}$ where no boundary conditions are imposed, since
we need to take integrals up to infinity. Fourier \emph{series} are
similar and work with periodic boundary conditions (lattice) or topology
(cylinder, torus, sphere). Both methods work only well for \emph{linear}
equations $\rightarrow$ no coupling between different Fourier harmonics.
\item Remark: Here we could see that notation introduces quite some ambiguity.
We often identify the coordinates of a point $\v x$ with the coordinate
system itself. If we have another point $\v k$ and switch to spherical
coordinates $(r,\tht,\ph)$ we cannot distinguish $r$ from $|\v x|$,
even if it could be the radius $|\v y|$ of another point $\v y$.
To resolve this, we gave different names to our coordinates $(k,\tht_{(k)},\ph_{(k)})$
in $k$-space. We could also have introduced the notation $r_{\v x}$
and $r_{\v k}$, which is non-standard and maybe even more confusing.
\end{itemize}
\subsection*{Problem 4: Electron plasma oscillations - the easier way}
see lecture notes
\subsection*{Problem 5: Electron plasma oscillations - the harder way}
We look at oscillations/waves in an unmagnetized plasma due to electrostatic
forces.
Approximations:
\begin{enumerate}
\item Ions are resting compared to fast electron dynamics.
\item Electrons modeled as fluid.
\item Electrostatic forces are stronger than thermodynamics forces (``cold
plasma'').
\item Linear expansion around equilibrium with subsonic flow and quasineutrality.
\end{enumerate}
Take mass continuity and momentum equation for electron fluid with
particles of mass $m_{e}$ and charge $q_{e}=-e$, as well as Poisson
equation for electric potential and constant and static ion density
$n_{i}$,
\begin{align}
m_{e}\frac{\partial n_{e}(\v x,t)}{\partial t}+\nabla\cdot\v p_{\text{fl}}(\v x,t) & =0,\label{eq:conti}\\
\frac{\partial\v p_{\text{fl}}(\v x,t)}{\partial t}+m_{e}\nabla\cdot(\v v(\v x,t)\v p_{\text{fl}}(\v x,t)) & =-\nabla p(\v x,t)+e\nabla\Phi(\v x,t),\\
\Delta\Phi(\v x,t) & =-\frac{e}{\varepsilon_{0}}(n_{i}-n_{e}(\v x,t)).\label{eq:poiss}
\end{align}
Here we could pull out $m_{e}$ as a constant, and replaced mass flux
$\v{\Gamma}_{m}$ by the identical kinematic momentum density $\v p^{fl}$.
In order to combine all three equations we now take the time derivative
of the continuity equation, and the divergence of the momentum equation
to obtain
\begin{align}
m_{e}\frac{\partial^{2}n_{e}(\v x,t)}{\partial t^{2}}+\frac{\partial}{\partial t}\nabla\cdot\v p_{\text{fl}}(\v x,t) & =0,\label{eq:conti-1}\\
\nabla\cdot\frac{\partial\v p_{\text{fl}}(\v x,t)}{\partial t}+\nabla\cdot(\nabla\cdot(\v v(\v x,t)\v p_{\text{fl}}(\v x,t))) & =-\Delta p(\v x,t)+en_{e}(\v x,t)\Delta\Phi(\v x,t),\\
\Delta\Phi(\v x,t) & =-\frac{e}{\varepsilon_{0}}(n_{i}-n_{e}(\v x,t)).\label{eq:poiss-1}
\end{align}
Exchanging divergence and time derivative in front of $\v p^{fl}$
we can eliminate the respective terms by subtracting the momentum
equation from the continuity equation. The Laplacian $\Delta\Phi$
is substituted from Poisson's equation to yield a single equation
\begin{align}
\frac{\partial^{2}n_{e}(\v x,t)}{\partial t^{2}}+\frac{e^{2}n_{e}(\v x,t)}{m_{e}\varepsilon_{0}}(n_{e}(\v x,t)-n_{i}(\v x))-\nabla\cdot(\nabla\cdot(\v v(\v x,t)\v v(\v x,t))) & =\frac{\Delta p(\v x,t)}{m_{e}}.
\end{align}
Now we start with a stationary equilibrium without time derivatives,
and flow velocity small enough (``subsonic'') such that the term
$\nabla\cdot(\nabla\cdot(\v v(\v x,t)\v v(\v x,t)))$ can be neglected.
The equilibrium should fulfill quasineutrality,
\begin{equation}
n_{0}(\v x)=n_{e}^{0}(\v x)=n_{i}^{0}(\v x),
\end{equation}
so $\Delta p^{0}=0$ in the subsonic approximation.
Now we introduce a perturbation in $n_{e}(\v x,t)=n_{e}^{0}(\v x)+\delta n_{e}(\v x,t)$
and $\v v=\v v^{0}+\delta\v v(\v x,t)$, but neglect changes in $p$
(``cold plasma''). Keeping only first order terms (remember, $n_{e}^{0}-n_{i}^{0}=0$),
we obtain
\begin{align}
\frac{\partial^{2}\delta n_{e}(\v x,t)}{\partial t^{2}}+\frac{e^{2}n_{0}(\v x)}{m_{e}\varepsilon_{0}}\delta n_{e}(\v x,t) & =0.
\end{align}
Here the electron plasma frequency appears locally at position $\v x$
as
\begin{equation}
\omega_{\text{pe}}^{\,2}(\v x)=\frac{e^{2}n_{0}(\v x)}{m_{e}\varepsilon_{0}}.
\end{equation}
\section{Collisions}
\subsection*{2.1 Binary Coulomb Collisions}
see lecture notes
\section{Charged Particle Dynamics}
\subsection*{3.2 Motion in static and homogeneous fields}
Our \textbf{goal} is to solve
\begin{align}
m\ddot{\v r}(t) & =\frac{e}{c}\dot{\v r}(t)\times\v B+e\v E.\label{eq:Lorentz}
\end{align}
or written in terms of a differential operator describing Lorentz
force:
\begin{equation}
\hat{L}\v r\equiv\left(\frac{\d^{2}}{\d t^{2}}+\left(\text{\ensuremath{\frac{e\v B}{mc}\frac{\d}{\d t}}}\right)\times\right)\v r(t)=\frac{e}{m}\v E\label{eq:LorentzOperator}
\end{equation}
and identify\emph{ gyromotion} circling perpendicular around $\v B$
with $v_{\perp}$,
\emph{parallel }motion along $\v B$ with $v_{\parallel}$ and
average \emph{drift} motion perpendicular to $\v B$ with $\v v_{D}$.\\
Our \textbf{strategy} is to separate out gyromotion in $\v B$-field
alone and solve its equations first. Then we could solve parallel
and drift motion with the constraint that $\v v_{D}$ remains constant.
In the combined solution we verified there are 6 free parameters that
can reproduce any initial conditions.\\
~
The\textbf{ assumption}\textbf{\emph{ }}is that we considere only
static and homogeneous fields. Under this simplification the solution
describes particle motion exactly, as $\hat{L}$ in Eq.~(\ref{eq:LorentzOperator})
is \emph{linear} in \textbf{$\v r(t)$}, so the overall solution can
be described as a linear superposition of two solutions.
Splitting Eq.~(\ref{eq:Lorentz}) into a sum of two solutions $\v R+\v{\rho}$
with $\v r=\v R+\v{\rho}$ yields
\begin{align}
m(\ddot{\v R}+\ddot{\v{\rho}}) & =e\v E+\frac{e}{c}(\dot{\v R}+\dot{\v{\rho}})\times\v B.\label{eq:split}
\end{align}
Since we introduced 3 more variables, we can also introduce 3 more
equations, as long as they are not in conflict with Eq. (\ref{eq:split}).
We choose
\begin{align}
m\ddot{\v{\rho}} & =\frac{e}{c}\dot{\v{\rho}}\times\v B,\\
m\ddot{\v R} & =e\v E+\frac{e}{c}\dot{\v R}\times\v B.
\end{align}
This way, the gyromotion in $\v{\rho}$ is only affected by the magnetic,
but not the electric field.
\subsubsection*{Solution for gyration}
First we solve $m\ddot{\v{\rho}}=\frac{e}{c}\dot{\v{\rho}}\times\v B$.
We choose our coordinate system such that field $\v B=B\v b=B\v e_{1}\times\v e_{2}$.
(or: $\v e_{3}=\v b$) Then
\begin{align}
m\ddot{\v{\rho}} & =\frac{e}{c}\dot{\v{\rho}}\times(B\v b)=\frac{eB}{c}(\dot{\rho}_{1}\v e_{1}+\dot{\rho}_{2}\v e_{2})\times\v e_{3}\nonumber \\
& =\frac{eB}{c}(-\dot{\rho}_{1}\v e_{2}+\dot{\rho}_{2}\v e_{1}).
\end{align}
We combine this to
\begin{align}
\ddot{\rho}_{1} & =\frac{eB}{mc}\dot{\rho}_{2}=\omega_{c}\dot{\rho}_{2},\\
\ddot{\rho}_{2} & =-\frac{eB}{mc}\dot{\rho}_{1}=-\omega_{c}\dot{\rho}_{1}.
\end{align}
Taking another time derivative and denoting $v_{1}\equiv\dot{\rho}_{1}$,
$v_{2}\equiv\dot{\rho}_{2}$ we obtain
\begin{equation}
\ddot{v}_{1}=\omega_{c}\dot{v}_{2}=-\omega_{c}^{\,2}v_{1}.
\end{equation}
In this second order equation we have six free constants. We choose
one to be the \textbf{perpendicular velocity} $v_{\perp}$, which
is the amplitude of the gyration velocity and another one as a phase-shift
$\phi_{0}$. Then we obtain the solution of position, velocity, and
acceleration
\begin{align}
\v{\rho} & =\frac{v_{\perp}}{\omega_{c}}(\sin(\omega_{c}t+\phi_{0})\,\v e_{1}+\cos(\omega_{c}t+\phi_{0})\,\v e_{2}),\\
\dot{\v{\rho}} & =v_{\perp}(\cos(\omega_{c}t+\phi_{0})\,\v e_{1}-\sin(\omega_{c}t+\phi_{0})\,\v e_{2}),\\
\ddot{\v{\rho}} & =-v_{\perp}\omega_{c}(\sin(\omega_{c}t+\phi_{0})\,\v e_{1}+\cos(\omega_{c}t+\phi_{0})\,\v e_{2}).
\end{align}
From here we define the \textbf{gyroradius} $\rho_{c}=v_{\perp}/\omega_{c}$
and the \textbf{gyrophase} $\phi=\omega_{c}t+\phi_{0}$.\\
\\
Additional four integration constants are a translational shift by
$\v{\rho}_{0}$ plus a constant parallel velocity $\dot{\v{\rho}}_{0}\cdot\v b$
which we both set to zero! This means that we didn't include any motion
parallel to the magnetic field, or drift across the field intro $\v{\rho}$
and need to account for them in $\v R$.\\
\\
(Have $6\times2$ integration constants in $\v{\rho},\v R$ but need
only $6$ of them for overall problem in $\v r$).
\subsubsection*{Solution for parallel and drift motion}
Now we split
\begin{equation}
\dot{\v R}=v_{\parallel}\v b+\v v_{D},\qquad\v v_{D}\cdot\v b=0.
\end{equation}
We need to solve
\begin{equation}
m\ddot{\v R}=\dot{v}_{\parallel}\v b+\dot{\v v}_{D}=e\v E+\frac{e}{c}\dot{\v R}\times\v B.
\end{equation}
For parallel motion, obtain immediately
\begin{equation}
m\dot{v}_{\parallel}=eE_{\parallel}\Rightarrow v_{\parallel}(t)=v_{\parallel0}+eE_{\parallel}t.
\end{equation}
For perpendicular drift use
\begin{equation}
m\dot{\v v}_{D}=e\v E_{\perp}+\frac{e}{c}\v v_{D}\times\v B
\end{equation}
take $\times\v B$ and use bac-cab rule (NRL)
\begin{equation}
m\dot{\v v}_{D}\times\v B=e\v E\times\v B-\frac{e}{c}B^{2}v_{D}.
\end{equation}
Consistent solution is possible if $\dot{\v v}_{D}\times\v B$ remains
zero. This way we do not include any gyration in the perpendicular
motion. Solution for this drift motion (\textbf{ExB drift}):
\begin{equation}
\v v_{D}=c\frac{\v E\times\v B}{B^{2}}.
\end{equation}
\emph{... independent of charge and mass! }This also means setting
initial condition $v_{D0}=c\frac{\v E\times\v B}{B^{2}}$ which fixes
another integration constant.
The solution for parallel and drift motion together is thus
\begin{equation}
\v R(t)=\v R_{0}+c\frac{\v E\times\v B}{B^{2}}t+\left(v_{\parallel0}t+\frac{eE_{\parallel}}{2}t^{2}\right)\v b.
\end{equation}
\subsubsection*{Combined solution}
The combined exact solution for particle motion in static, homogeneous
$\v E$ and $\v B$ is
\begin{align}
\v r(t) & =\v R(t)+\v{\rho}(t)\nonumber \\
& =\v R_{0}+c\frac{\v E\times\v B}{B^{2}}t+\left(v_{\parallel0}t+\frac{eE_{\parallel}}{2}t^{2}\right)\v b+\frac{v_{\perp}}{\omega_{c}}\left(\sin(\omega_{c}t+\phi_{0})\,\v e_{1}+\cos(\omega_{c}t+\phi_{0})\,\v e_{2}\right).
\end{align}
~
Let's count integration constants:
\begin{itemize}
\item 1,2,3: $\v R_{0}$ contains 3 constants of initial position of Larmor
circle center (\textbf{gyrocenter}).
\item 4: Constant $v_{\parallel0}$ sets initial gyrocenter parallel velocity.
\item 5,6: Constants $v_{\perp}$ and $\phi_{0}$ fix the gyromotion ($v_{\perp}$
can optionally be replaced by $\rho_{c}$).\\
\end{itemize}
Thus we could find a consistent solution that separates \emph{gyromotion}
from \emph{parallel} and \emph{drift} motion in \emph{static} and
\emph{homogeneous} fields. In more complicated fields we will rely
on methods of averaging that work as long as gyromotion i.e. $\omega_{c}$
is fast, and $\rho_{c}$ is small, compared to field scales.
\subsection*{3.4 Motion in weakly inhomogeneous fields\label{subsec:Inhomo}}
Now the \textbf{goal} is to solve
\begin{align}
m\ddot{\v r} & =\frac{e}{c}\dot{\v r}\times\v B(\v r)+e\v E(\v r).
\end{align}
or written in terms of a \emph{non-linear} Lorentz operator:
\begin{equation}
\hat{L}\v r\equiv\left(\frac{\d^{2}}{\d t^{2}}+\frac{e\v B(\v r)}{mc}\frac{\d}{\d t}\times\right)\v r=\frac{e}{m}\v E(\v r).
\end{equation}
\textbf{Strategy:} Instead of an exact splitting as before, separate
into \emph{fast} gyromotion for $\v{\rho}$ and \emph{gyro-averaged}
motion in $\v R$. The latter means to take time averages over\emph{
gyroperiod} $\tau_{c}=2\pi/\omega_{c}$ and spatial averages within
\emph{gyroradius $\rho_{c}=v_{\perp}/\omega_{c}$. }Perform a series
expansion of fields around average $\v R$ in direction of $\v{\rho}$.\\
~
We\textbf{ assume}\textbf{\emph{ }}that relative variations of $\v E$
and $\v B$ in time and space are small over $\tau_{c}$ and within
$\rho_{c}$. This means that the plasma must be \emph{strongly magnetized
}to get reasonably large $\omega_{c}=eB/(mc)$. Due to the expansion
the solution is only \textbf{approximate}\emph{.}
\subsubsection*{Calculation}
A detailed calculation can be found in the lecture notes, chapter
3.4. For details how to do an expansion of tensor fields, see p.~\pageref{subsec:Taylor-expansion-of}.
\subsubsection*{Result}
The orbital motion is split into the sum of the gyrocenter position
$\v R$ and the Larmor gyration vector $\v{\rho}$, where
\begin{align}
\v{\rho} & =\frac{v_{\perp}}{\omega_{c}(\v R)}(\sin(\phi(\v R,t))\,\v e_{1}(\v R)+\cos(\phi(\v R,t))\,\v e_{2}(\v R)),\\
m\ddot{\v R} & =e\left(1+\frac{\rho_{c}^{2}}{4}\nabla_{\perp}^{\,2}\right)\v E(\v R)+\frac{e}{c}\dot{\v R}\times\v B(\v R)-\mu\nabla B(\v R).\label{eq:R}
\end{align}
Fields are evaluated in the gyrocenter position $\v R$ that can be
computed ignoring $\v{\rho}$. The magnetic moment $\mu$ is a conserved
quantity and related to $v_{\perp}$ via
\begin{equation}
\mu=\frac{mv_{\perp}^{\,2}}{2B(\v R)}=\frac{ev_{\perp}}{2\pi\rho_{c}}\pi\rho_{c}^{\,2}.\label{eq:magmoment}
\end{equation}
This is the ring current over one gyration period and $\mu=I\cdot A$
over the gyration circle with surface area $A=\pi\rho_{c}^{2}$. If
$E_{\parallel}=\v b\cdot\v E$ then $\mu$ is approximately conserved
-- an \textbf{adiabatic invariant}.
\begin{align}
\frac{\d\mu}{\d t} & =\frac{1}{B}\frac{\d}{\d t}\left(\frac{mv_{\perp}^{\,2}}{2}\right)-\frac{mv_{\perp}^{\,2}}{2B^{2}}\dot{\v R}\cdot\nabla B\nonumber \\
& \approx\frac{1}{B}\frac{\d}{\d t}\left(\frac{mv_{\perp}^{\,2}}{2}\right)-\frac{\mu}{B}v_{\parallel}\v b\cdot\nabla B.
\end{align}
Compare this with total energy conservation
\begin{align}
0 & =\frac{\d}{\d t}\left(\frac{mv_{\parallel}^{\,2}}{2}+\frac{mv_{\perp}^{\,2}}{2}+e\Phi\right)\nonumber \\
& \approx\frac{\d}{\d t}\left(\frac{mv_{\perp}^{\,2}}{2}\right)-v_{\parallel}\mu\v b\cdot\nabla B.
\end{align}
Here $v_{\parallel}\v b\cdot\nabla\Phi$ has vanished due to $E_{\parallel}=0$
and $\dot{v}_{\parallel}$ has been replaced by its expression from
the equations of motion (\ref{eq:R}) with $\dot{\v R}\approx v_{\parallel}\v b$.
Together we see that the total time derivative of $\mu$ vanishes
in this first order calculation. A subtle point that should concern
us is the actual nature of $v_{\perp}$ and $\phi$. If we assume
$\mu$ to be exactly conserved, Eq.~(\ref{eq:magmoment}) becomes
a definition for $v_{\perp}=v_{\perp}(\v R)$ rather than a quantity
resulting from it. While the motion of the gyrocenter $\v R$ can
be computed from Eq.~(\ref{eq:R}) without referring to the gyromotion,
the definition what $\v R$ actually describes is still somewhat dubious.
For this reason we will perform a more ``watertight'' treatment
of the problem using Lagrangian mechanics.
\section{Guiding-Center Dynamics}
The modern way to treat guiding-center motion via a phase-space Lagrangian
has been introduced by Littlejohn~\citep{littlejohn1983} in the
1980s. Besides that there exists a review article on guiding-center
motion by Cary and Brizard \citep{cary2009} containing background,
main derivation steps and further literature references. A detailed
derivation of the guiding-center Lagrangian can also be found in Lectures
3-4 of 2017 and https://itp.tugraz.at/\textasciitilde ert/assets/internal/plasma/2017/A1\_Guiding\_Center\_Lagrangian.pdf
. More information on asymptotic expansions with a small parameter
$\varepsilon$ can be found in the lecture of Stefan Possanner at
TUM. Here we treat the case where fields are not time-dependent.
We transform from the \textbf{particle} coordinates $\v z_{P}=(\v r,\v v)$
in phase-space to new \textbf{guiding-center} coordinates $\v z=(\v R,\phi,v_{\parallel},v_{\perp})$
related to particle position $\v r$ and velocity $\v v$ via
\begin{align}
\v r(\v z) & =\v R+\v{\rho}(\v R,\phi,v_{\perp}),\text{ where }\v{\rho}(\v R,\phi,v_{\perp})\equiv\frac{mcv_{\perp}}{eB(\v R)}\hat{\v{\rho}}(\v R,\phi)\equiv\rho_{c}(v_{\perp},\v R)\hat{\v{\rho}}(\v R,\phi),\label{eq:rz}\\
\v v(\v z) & =v_{\parallel}\hat{\v b}(\v R)+v_{\perp}\hat{\v n}(\v R,\phi).\label{eq:vz}
\end{align}
Here $\hat{\v b}(\v R)=\v B(\v R)/B(\v R)$ and forms a tripod with
some pair of \textbf{fixed }but possibly space-dependent unit vectors
$\hat{\v e}_{1}(\v R)$ and $\hat{\v e}_{2}(\v R)$, and definitions
for \textbf{dynamic} (gyrating) unit vectors $\hat{\v{\rho}}$ and
$\hat{\v n}$ are
\begin{align}
\hat{\v{\rho}}(\v R,\phi) & =\cos(\phi)\,\hat{\v e}_{1}(\v R)-\sin(\phi)\,\hat{\v e}_{2}(\v R)=-\frac{\partial\hat{\v n}(\v R,\phi)}{\partial\phi},\\
\hat{\v n}(\v R,\phi) & =\hat{\v{\rho}}(\v R,\phi)\times\hat{\v b}(\v R)=-\sin(\phi)\,\hat{\v e}_{1}(\v R)-\cos(\phi)\,\hat{\v e}_{2}(\v R)=\frac{\partial\hat{\v{\rho}}(\v R,\phi)}{\partial\phi}.
\end{align}
The transformation of Eqs.~(\ref{eq:rz}-\ref{eq:vz}) is still exact!
We call $\v R$ the guiding-center position. Both, $(\hat{\v e}_{1},\hat{\v e}_{2},\hat{\v b})$
and $(\hat{\v n},\hat{\v{\rho}},\hat{\v b})$ form an orthonormal
basis, where $\hat{\v e}_{1},\hat{\v e}_{2}$ depend only on $\boldsymbol{R}$,
but $\hat{\v n},\hat{\v{\rho}}$ also on $\phi$. Variables have suspicious
names $\phi$, $v_{\parallel}$ and $v_{\perp}$. Indeed, for the
case of a homogeneous $\v B$-field those variables are equal to gyrophase,
parallel and perpendicular velocity of the particle. If the field
becomes weakly inhomogeneous, they still approximate the particle
values averaged over one gyro-period. Similarly, $\v R$ describes
the center of the Larmor circle only in the homogeneous limiting case,
and remains close to it in a field with sufficiently weak gradients.
To sum up one could say that $(\phi,v_{\parallel},v_{\perp})$ are
values we would get if we pretend that $\v B$ is evaluated in the
guiding-center position $\v R$ rather than the particle position
$\v R+\v{\rho}$. In contrast to the treatment in section~\ref{subsec:Inhomo},
all variables in $\v z=(\v R,\phi,v_{\parallel},v_{\perp})$ have
a well-defined relation to the particle position $\v r$ and velocity
$\v v$ if Eqs.~(\ref{eq:rz}-\ref{eq:vz}) are solved implicitly
in $\v z$.
The phase-space Lagrangian of a particle in an electromagnetic field
in terms of canonical $(\v r,\v p$) is given by
\begin{equation}
L(\v r,\dot{\v r},\v p,\cancel{\dot{\v p}})=\v p\cdot\dot{\v r}-\left(\frac{(\v p-\frac{e}{c}\v A(\v r))^{2}}{2m}+e\Phi(\v r)\right).\label{eq:Lag1-1}
\end{equation}
Using non-canonical variables $\v v=\frac{\v p}{m}-\frac{e}{mc}\v A$
instead of the momenta $\v p$ the particle phase-space Lagrangian
is
\begin{equation}
L(\v r,\dot{\v r},\v v,\cancel{\dot{\v v}})=\left(m\v v+\frac{e}{c}\v A(\v r)\right)\cdot\dot{\v r}-\left(\frac{m\v v^{2}}{2}+e\Phi(\v r)\right).\label{eq:Lag1-1-1}
\end{equation}
In terms of $\v z=(\v R,\phi,v_{\parallel},v_{\perp})$ we have
\begin{equation}
L(\v R,\dot{\v R},\phi,\dot{\phi},v_{\parallel},v_{\perp})=\left(mv_{\parallel}\hat{\v b}(\v R)+mv_{\perp}\hat{\v n}(\v R,\phi)+\frac{e}{c}\v A(\v R+\v{\rho})\right)\cdot(\dot{\v R}+\dot{\v{\rho}})-\left(\frac{m(v_{\parallel}^{\,2}+v_{\perp}^{\,2})}{2}+e\Phi(\v R+\v{\rho})\right),\label{eq:Lag1-1-1-1}
\end{equation}
with gyration vector $\v{\rho}=\v{\rho}(\v R,\phi,v_{\perp})=\v r-\v R$
from Eq.~(\ref{eq:rz}). This Lagrangian is still exact, as are transformations
to guiding-center variables. We will now make approximations to the
Lagrangian, affecting the \emph{dynamics} of motion in guiding-center
variables. As in the previous chapter we use a Taylor expansion which
we break at the first order. For higher order computations, Lie transform
methods are better suited (see Littlejohn, 1983). Here the small parameter
compared to typical length scales is the gyroradius
\begin{equation}
\rho_{c}=|\v{\rho}|=\frac{mcv_{\perp}}{eB(\v R)}\sim\varepsilon.
\end{equation}
Furthermore we define the cyclotron frequency at the guiding-center
position
\begin{equation}
\omega_{c}=\frac{eB(\v R)}{mc}\sim\varepsilon^{-1},
\end{equation}
as a large parameter compared to typical time scales. This corresponds
to the physical conditions where the guiding-center approximation
is applicable, leaves the perpendicular velocity $v_{\perp}$ being
a quantity of interest with order $\varepsilon^{0}=1$ and thus prevents
it from being swallowed by a perturbation expansion. The parameter
$\varepsilon$ is called an ordering parameter that allows us to collect
terms of similar magnitude. Vector and scalar fields are expanded
in the following way,
\begin{align}
\v A(\v r) & =\v A(\v R)+\varepsilon\v{\rho}\cdot\nabla\v A(\v R)+\mathcal{O}(\varepsilon^{2}),\\
\Phi(\v r) & =\Phi(\v R)+\varepsilon\v{\rho}\cdot\nabla\Phi(\v R)+\mathcal{O}(\varepsilon^{2}),
\end{align}
where $\varepsilon$ acts as a marker for intermediate calculations
and is set to one only in the final result after neglecting all terms
of order $\varepsilon$ or higher. We mark each small quantity with
an $\varepsilon$ in Eq.~\ref{eq:Lag1-1-1-1} to obtain
\begin{equation}
L(\v R,\dot{\v R},\phi,\dot{\phi},v_{\parallel},v_{\perp})=\left(mv_{\parallel}\hat{\v b}(\v R)+mv_{\perp}\hat{\v n}(\v R,\phi)+\varepsilon^{-1}\frac{e}{c}\v A(\v R+\varepsilon\v{\rho})\right)\cdot(\dot{\v R}+\dot{\v{\rho}})-\left(\frac{m(v_{\parallel}^{\,2}+v_{\perp}^{\,2})}{2}+e\Phi(\v R+\varepsilon\v{\rho})\right).\label{eq:Lag1-1-1-1-1}
\end{equation}
We cannot simply write $\varepsilon$ in front of $\dot{\boldsymbol{\rho}}$
just because the magnitude of $\v{\rho}$ is small, since gyration
is fast. This becomes apparent by using the chain rule,
\begin{align}
\dot{\v{\rho}} & =\frac{\d}{\d t}\v{\rho}(\v R,\phi,v_{\perp})=\dot{\v R}\cdot\nabla\v{\rho}+\dot{\phi}\frac{\partial\v{\rho}}{\partial\phi}+\dot{v}_{\perp}\frac{\partial\v{\rho}}{\partial v_{\perp}}\nonumber \\
& =\varepsilon\dot{\v R}\cdot\nabla\v{\rho}+\varepsilon^{-1}\dot{\phi}\varepsilon\rho_{c}\hat{\v n}+\varepsilon\frac{\dot{v}_{\perp}}{v_{\perp}}\v{\rho}=\dot{\phi}\rho_{c}\hat{\v n}+\varepsilon(\dot{\v R}\cdot\nabla\v{\rho}+\frac{\dot{v}_{\perp}}{v_{\perp}}\v{\rho}).
\end{align}
Here in the first term the fast gyration $\dot{\phi}$ and the small
$\rho_{c}$ cancel to order $\varepsilon^{0}$, and the last term
should be even smaller than order $\varepsilon$ as long as $v_{\perp}$
doesn't change too fast due to strongly inhomogeneous fields. Now
we look at the first term in~(\ref{eq:Lag1-1-1-1-1}) up to order
$\varepsilon^{0}$,
\begin{align}
\left(mv_{\parallel}\hat{\v b}(\v R)+mv_{\perp}\hat{\v n}(\v R,\phi)+\varepsilon^{-1}\frac{e}{c}\v A(\v R+\v{\rho})\right)\cdot(\dot{\v R}+\dot{\v{\rho}}) & =\left(mv_{\parallel}\hat{\v b}+mv_{\perp}\hat{\v n}+\varepsilon^{-1}\frac{e}{c}\v A+\frac{e}{c}\v{\rho}\cdot\nabla\v A\right)\cdot\dot{\v R}\nonumber \\
& +mv_{\perp}\dot{\phi}\rho_{c}+\varepsilon^{-1}\frac{e}{c}\v A\cdot\dot{\v{\rho}}+\frac{e}{c}(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v{\rho}}+\mathcal{O}(\varepsilon).
\end{align}
Here the term with $\hat{\v b}\cdot\hat{\v n}$ vanished and we took
a scalar product $\hat{\v n}\cdot\hat{\v n}=1$ to obtain the first
term in the second line. Now we use the fact that we can add/subtract
an arbitrary total time derivative to the Lagrangian without changing
its dynamics. We will use this trick two times in a row.
First we take
\begin{align}
\frac{\d}{\d t}(\v{\rho}\cdot\v A(\v R)) & =\dot{\v{\rho}}\cdot\v A+\v{\rho}\cdot\dot{\v A}\nonumber \\
& =\dot{\v{\rho}}\cdot\v A+\v{\rho}\cdot(\dot{\v R}\cdot\nabla\v A).\label{eq:dt}
\end{align}
This we can combine with the term $(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v R}=\v{\rho}\cdot((\nabla\v A)\cdot\dot{\v R})$
via tensor algebra rule (\ref{eq:eq3}),
\begin{equation}
(\nabla\v A)\cdot\dot{\v R}-\dot{\v R}\cdot\nabla\v A=\dot{\v R}\times(\nabla\times\v A)=\dot{\v R}\times\boldsymbol{B}.\label{eq:eq3-1}
\end{equation}
So
\begin{align}
(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v R}-\v{\rho}\cdot(\dot{\v R}\cdot\nabla\v A) & =\v{\rho}\cdot(\dot{\v R}\times\v B).
\end{align}
and together with~(\ref{eq:dt}) we obtain
\begin{equation}
(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v R}=\v{\rho}\cdot(\dot{\v R}\times\v B+\dot{\v R}\cdot\nabla\v A)=\v{\rho}\cdot(\dot{\v R}\times\v B)+\frac{\d}{\d t}(\v{\rho}\cdot\v A)-\dot{\v{\rho}}\cdot\v A.
\end{equation}
Inserting everything into $L$ cancels the term $\v A\cdot\dot{\v{\rho}}$
and we obtain
\begin{align}
L-\frac{e}{c}\frac{\d}{\d t}(\v{\rho}\cdot\v A)+\mathcal{O}(\varepsilon) & =\left(mv_{\parallel}\hat{\v b}+mv_{\perp}\hat{\v n}+\varepsilon^{-1}\frac{e}{c}\v A\right)\cdot\dot{\v R}\nonumber \\
& +\frac{e}{c}\v{\rho}\cdot(\dot{\v R}\times\v B)+mv_{\perp}\dot{\phi}\rho_{c}+\frac{e}{c}(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v{\rho}}\nonumber \\
& -\left(\frac{m(v_{\parallel}^{\,2}+v_{\perp}^{\,2})}{2}+e\Phi\right),
\end{align}
where $\v A,\v B,$ and $\Phi$ are now evaluated in the guiding-center
position $\v R$ instead of the particle position $\v r$. Two terms
cancel since by rotating the triple-product and extracting basis vectors,
\begin{equation}
\frac{e}{c}\v{\rho}\cdot(\dot{\v R}\times\v B)=\frac{e}{c}\rho_{c}B\dot{\boldsymbol{R}}\cdot(\hat{\v b}\times\hat{\v{\rho}})=-mv_{\perp}\dot{\boldsymbol{R}}\cdot\hat{\v n}.
\end{equation}
To eliminate the term $(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v{\rho}}=\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}}$
we use another total time derivative (see also review paper of Cary
and Brizard, 2009),
\begin{align}
\frac{1}{2}\frac{\d}{\d t}\left(\varepsilon\v{\rho}\cdot\varepsilon^{-1}\frac{e}{c}\nabla\v A\cdot\varepsilon\v{\rho}\right) & =\frac{1}{2}\left(\dot{\v{\rho}}\cdot\nabla\v A\cdot\v{\rho}+\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}}\right)\nonumber \\
& +\frac{\varepsilon}{2}\v{\rho}\cdot\left(\frac{\d\nabla\v A}{\d t}\right)\cdot\v{\rho}.
\end{align}
Despite being a total time derivative one could also argue that this
term is of order $\varepsilon$ for the leading order approximation.
Since it contains terms of higher order due to the time derivative
and fast gyration one should be careful with this. Anyway, we can
use the expression for the term
\begin{align}
(\v{\rho}\times\dot{\v{\rho}})\cdot\v B & =(\v{\rho}\times\dot{\v{\rho}})\cdot(\nabla\times\v A)\nonumber \\
& =(\varepsilon_{ijk}\rho_{j}\dot{\rho}_{k})(\varepsilon_{imn}\partial_{m}A_{n})\nonumber \\
& =(\delta_{jm}\delta_{kn}-\delta_{jn}\delta_{km})(\rho_{j}\dot{\rho}_{k}\partial_{m}A_{n})\nonumber \\
& =\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}}-\dot{\v{\rho}}\cdot\nabla\v A\cdot\v{\rho}.
\end{align}
Thus,
\begin{align}
\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}} & =\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}}=-\dot{\v{\rho}}\cdot\nabla\v A\cdot\v{\rho}\nonumber \\
& =\frac{1}{2}\underbrace{\left(\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}}-\dot{\v{\rho}}\cdot\nabla\v A\cdot\v{\rho}\right)}_{(\v{\rho}\times\dot{\v{\rho}})\cdot\v B}+\frac{1}{2}\underbrace{\left(\dot{\v{\rho}}\cdot\nabla\v A\cdot\v{\rho}+\v{\rho}\cdot\nabla\v A\cdot\dot{\v{\rho}}\right)}_{\mathcal{O}(\varepsilon)-\frac{1}{2}\frac{\d}{\d t}\dots}.
\end{align}
Further,
\begin{equation}
(\v{\rho}\times\dot{\v{\rho}})\cdot\v B=\left(\rho_{c}\hat{\v{\rho}}\times\dot{\phi}\rho_{c}\hat{\v n}\right)\cdot\boldsymbol{B}=-\dot{\phi}\rho_{c}^{\,2}B=-\frac{mcv_{\perp}}{e}\rho_{c}\dot{\phi}.
\end{equation}
and we can replace $(\v{\rho}\cdot\nabla\v A)\cdot\dot{\v{\rho}}$
in the Lagrangian to obtain to order $\varepsilon$ (formally setting
$\varepsilon=1$)
\begin{align}
L-\frac{e}{c}\frac{\d}{\d t}\left(\v{\rho}\cdot\v A-\frac{1}{2}\v{\rho}\cdot\nabla\v A\cdot\v{\rho}\right)+\mathcal{O}(\varepsilon) & =\left(mv_{\parallel}\hat{\v b}+\frac{e}{c}\v A\right)\cdot\dot{\v R}\nonumber \\
& +mv_{\perp}\dot{\phi}\rho_{c}-\frac{mv_{\perp}}{2}\rho_{c}\dot{\phi}\nonumber \\
& -\left(\frac{m(v_{\parallel}^{\,2}+v_{\perp}^{\,2})}{2}+e\Phi\right).
\end{align}
In the bracket under the total time derivative we can see that there
is some system to the way the gauge transform in phase-space is done.
More generally the technique is known as the \textbf{Lie-transform
method}. To leading order in $\varepsilon$ the \textbf{guiding-center
Lagrangian} follows as
\[
L_{\mathrm{gc}}(\v R,\dot{\v R},\cancel{\phi},\dot{\phi},v_{\parallel},v_{\perp})=\left(mv_{\parallel}\hat{\v b}(\v R)+\frac{e}{c}\v A(\v R)\right)\cdot\dot{\v R}+\frac{m^{2}cv_{\perp}^{\,2}}{2eB(\v R)}\dot{\phi}-\left(\frac{m(v_{\parallel}^{\,2}+v_{\perp}^{\,2})}{2}+e\Phi(\v R)\right).
\]
Here all fields are evaluated in $\v R$ and dependencies on $\phi$
have vanished. Thus we can identify a conserved generalized momentum
\begin{equation}
\frac{\partial L_{\mathrm{gc}}}{\partial\dot{\phi}}=\frac{m^{2}cv_{\perp}^{\,2}}{2eB(\v R)}\equiv J_{\perp}=\frac{mc}{e}\mu,
\end{equation}
known as the perpendicular invariant with $\dot{J}_{\perp}=0$ along
the trajectory. The perpendicular invariant is proportional to the
magnetic moment $\mu$ of the guiding-center orbit, which is thus
also conserved. Again, $\mu$ is not the exact magnetic moment of
the particle, but rather an approximation up to order $\varepsilon$.
If we replace $v_{\perp}$ by the variable $\mu$ we can write
\begin{equation}
L_{\mathrm{gc}}(\v R,\dot{\v R},\dot{\phi},v_{\parallel},\mu)=\left(mv_{\parallel}\hat{\v b}(\v R)+\frac{e}{c}\v A(\v R)\right)\cdot\dot{\v R}-\frac{mc}{e}\mu\dot{\phi}-\left(\frac{mv_{\parallel}^{\,2}}{2}+\mu B(\v R)+e\Phi(\v R)\right).
\end{equation}
The Hamiltonian is the last part in brackets,
\begin{equation}
H=\frac{mv_{\parallel}^{\,2}}{2}+\mu B(\v R)+e\Phi(\v R),
\end{equation}
and is equal to the total energy. For equations of motion, it is of
limited use, since we use non-canonical variables. We rather need
Euler-Lagrange equations of the guiding-center Lagrangian that follow
as
\begin{align}
\dot{\v R} & =\frac{v_{\parallel}\v B^{\star}(\v R,v_{\parallel})+c\v E^{\star}(\v R,\mu)\times\hat{\v b}(\v R)}{\hat{\v b}(\v R)\cdot\v B^{\star}(\v R,v_{\parallel})},\quad\dot{v}_{\parallel}=\frac{e}{m}\frac{\v E^{\star}(\v R,\mu)\cdot\v B^{\star}(\v R,v_{\parallel})}{\hat{\v b}(\v R)\cdot\v B^{\star}(\v R,v_{\parallel})},\label{eq:rdot-1}\\
\dot{\phi} & =-\omega_{c}(\v R)\equiv\frac{e\v B(\v R)}{mc},\quad\dot{\mu}=0.\label{eq:phidot}
\end{align}
where $\hat{\v b}(\v R)=\v B(\v R)/B(\v R)$ and
\begin{align}
\v E^{\star}(\v R,\mu) & =-\frac{1}{e}\nabla H(\v R,v_{\parallel},\mu)=-\frac{\mu}{e}\nabla B(\v R)-\nabla\Phi(\v R),\\
\v B^{\star}(\v R,v_{\parallel}) & =\nabla\times\v A^{\star}(\v R,v_{\parallel}),\\
\v A^{\star}(\v R,v_{\parallel}) & =\v A(\v R)+\frac{mc}{e}v_{\parallel}\hat{\v b}(\v R).
\end{align}
\textbf{Derivation of equations of motion}
From now on we will call the guiding-center Lagrangian $L_{\mathrm{gc}}$
just $L$ to reduce notational clutter. Take
\begin{equation}
L(\v R,v_{\parallel},\mu,\dot{\v R},\dot{\phi})=\left(mv_{\parallel}\hat{\v b}(\v R)+\frac{e}{c}\v A(\v R)\right)\cdot\dot{\v R}-\frac{mc}{e}\mu\dot{\phi}-\left(\frac{mv_{\parallel}^{\,2}}{2}+\mu B(\v R)+e\Phi(\v R)\right).
\end{equation}
Derivatives are
\begin{align}
\frac{\partial L}{\partial\v R} & =mv_{\parallel}(\nabla\hat{\v b})\cdot\dot{\v R}+\frac{e}{c}(\nabla\v A)\cdot\dot{\v R}-\left(\mu\nabla B+e\nabla\Phi\right)\\
\frac{\partial L}{\partial\phi} & =0\\
\frac{\partial L}{\partial v_{\parallel}} & =m\hat{\v b}\cdot\dot{\v R}-mv_{\parallel}\label{eq:vpar}\\
\frac{\partial L}{\partial\mu} & =-\frac{mc}{e}\dot{\phi}-B\label{eq:mu}
\end{align}
Generalized velocities:
\begin{align}
\frac{\d}{\d t}\frac{\partial L}{\partial\dot{\v R}} & =\frac{\d}{\d t}(mv_{\parallel}\hat{\v b}(\v R(t))+\frac{e}{c}\v A(\v R(t)))=m\dot{v}_{\parallel}\hat{\v b}+mv_{\parallel}\dot{\v R}\cdot\nabla\hat{\v b}+\frac{e}{c}\dot{\v R}\cdot\nabla\v A.\\
\frac{\d}{\d t}\frac{\partial L}{\partial\dot{\phi}} & =\frac{\d}{\d t}(-\frac{mc}{e}\mu)=-\frac{mc}{e}\dot{\mu}
\end{align}
Euler-Lagrange euqations in $\v R$ are
\begin{align}
0=\frac{\d}{\d t}\frac{\partial L}{\partial\dot{\v R}}-\frac{\partial L}{\partial\v R} & =m\dot{v}_{\parallel}\hat{\v b}+mv_{\parallel}(\nabla\hat{\v b})\cdot\dot{\v R}+\frac{e}{c}(\nabla\v A)\cdot\dot{\v R}\nonumber \\
& -\mu\nabla B-e\nabla\Phi-mv_{\parallel}\dot{\v R}\cdot\nabla\hat{\v b}-\frac{e}{c}\dot{\v R}\cdot\nabla\v A\nonumber \\
& =m\dot{v}_{\parallel}\hat{\v b}+mv_{\parallel}\dot{\v R}\times(\nabla\times\hat{\v b})+\frac{e}{c}\dot{\v R}\times(\nabla\times\v A)\nonumber \\
& -\mu\nabla B+e\v E-mv_{\parallel}\dot{\v R}\cdot\nabla\hat{\v b}\nonumber \\
& =m\dot{v}_{\parallel}\hat{\v b}+\frac{e}{c}\dot{\v R}\times(\nabla\times\v A^{\star})-\mu\nabla B+e\v E-mv_{\parallel}\dot{\v R}\cdot\nabla\hat{\v b}\nonumber \\
& =m\dot{v}_{\parallel}\hat{\v b}+\frac{e}{c}\dot{\v R}\times\v B^{\star}-\mu\nabla B+e\v E-mv_{\parallel}\dot{\v R}\cdot\nabla\hat{\v b}.
\end{align}
Here we have used $\v a\times(\nabla\times\v b)=(\nabla\v a)\cdot\v b-\v b\cdot\nabla\v a$
from (\ref{eq:eq3}) and defined
\begin{align}
\v A^{\star}(\text{\ensuremath{\v R}},v_{\parallel}) & =\v A(\text{\ensuremath{\v R}})+\frac{mc}{e}v_{\parallel}\hat{\v b}(\text{\ensuremath{\v R}}),\\
\v B^{\star}(\text{\ensuremath{\v R}},v_{\parallel}) & =\nabla\times\v A^{\star}(\text{\ensuremath{\v R}},v_{\parallel})
\end{align}
as a modified vector potential and magnetic field. The latter depend
on the parallel velocity $v_{\parallel}$ that can take both, positive
or negative sign and vary during the orbit. The Euler-Lagrange equation
in the ignorable $\phi$ yield conservation of magnetic moment $\dot{\mu}=0$.
The remaining equations have no generalized velocity terms in $L_{\text{gc}}$
being a phase-space Lagrangian. Thus we can set terms (\ref{eq:vpar}-\ref{eq:mu})
to zero right away. To sum up our equations of motion become
\begin{align}
m\dot{v}_{\parallel}\hat{\v b}(\v R)+\frac{e}{c}\dot{\v R}\times\v B^{\star}(\text{\ensuremath{\v R}},v_{\parallel}) & =\mu\nabla B(\v R)-e\v E(\v R),\label{eq:elg1}\\
\hat{\v b}(\v R)\cdot\dot{\v R} & =v_{\parallel},\label{eq:elg2}\\
\dot{\phi} & =-\frac{eB(\v R)}{mc}\equiv-\omega_{c}(\v R).\label{eq:elg3}
\end{align}
Eq.~(\ref{eq:elg1}) describes drifts due to $\v E$ and gradient
and curvature of $\v B$. Eq.~(\ref{eq:elg2}) confirms the meaning
of $v_{\parallel}$ as the parallel velocity, and Eq.~(\ref{eq:elg3})
of $\phi$ as the gyrophase with regard to the guiding-center position
$\v R$. Now we would like to decouple Eqs.~(\ref{eq:elg1}-\ref{eq:elg2})
to get explicit expressions for $\dot{\v R}$ and $\dot{v}_{\parallel}$.
Taking a cross product of $\hat{\v b}$ with Eq.~\ref{eq:elg1} removes
the first term and together with the BAC-CAB rule and $\v b\cdot\dot{\v R}=v_{\parallel}$
from Eq.~(\ref{eq:elg2}) yields
\begin{equation}
\frac{e}{c}((\hat{\v b}\cdot\v B^{\star})\dot{\v R}-v_{\parallel}\v B^{\star})=\hat{\v b}\times(\mu\nabla B-e\v E).
\end{equation}
Finally we obtain
\begin{align}
\dot{\v R} & =\frac{v_{\parallel}\v B^{\star}+c\v E^{\star}\times\hat{\v b}}{\hat{\v b}\cdot\v B^{\star}},\label{eq:rdot}\\
\dot{v}_{\parallel} & =\frac{e}{m}\frac{\v B^{\star}\cdot\v E^{\star}}{\hat{\v b}\cdot\v B^{\star}}.
\end{align}
Here (convention taken from \citep{lanthaler2017}) we defined
\begin{equation}