# Action-angle coordinates for motion in a straight magnetic field with
constant gradient

Alain J. Brizard  Department of Physics, Saint Michael’s College, Colchester,
VT 05439, USA

###### Abstract

The motion of a charged particle in a straight magnetic field ${\bf
B}=B(y)\,\widehat{\sf z}$ with a constant perpendicular gradient is solved
exactly in terms of elliptic functions and integrals. The motion can be
decomposed in terms of a periodic motion along the $y$-axis and a drift motion
along the $x$-axis. The periodic motion can be described as a particle trapped
in a symmetric quartic potential in $y$. The canonical transformation from the
canonical coordinates $(y,P_{y})$ to the action-angle coordinates $(J,\theta)$
is solved explicitly in terms of a generating function $S(\theta,J)$ that is
expressed in terms of Jacobi elliptic functions. The presence of a weak
constant electric field ${\bf E}=E_{0}\,\widehat{\sf y}$ introduces an
asymmetric component to the quartic potential, and the associated periodic
motion is solved perturbatively up to second order.

##  I Introduction

The motion of a charged particle in a straight magnetic field ${\bf
B}=B(y)\,\widehat{\sf z}$ with constant perpendicular gradient was considered
recently Brizard_2017  in order to explore the validity of the guiding-center
approximation Cary_Brizard_2009  in the present of strong gradients. There, a
single orbit in the $(x,y)$ plane was solved exactly in terms of elliptic
functions and integrals footnote , and the orbit-averaged position $\langle
y\rangle$ and drift velocity $\langle\dot{x}\rangle$ were compared,
respectively, with the polarization shift and magnetic-drift velocity
predicted by guiding-center Hamiltonian theory Cary_Brizard_2009 . An
excellent agreement was found between exact analytical results and guiding-
center predictions, even in the presence of strong magnetic gradients.

The purpose of the present paper is to consider all orbit types associated
with the straight magnetic field ${\bf B}=B_{0}\,(1-y/L)\,\widehat{\sf z}$
considered in Ref. Brizard_2017 , where $B_{0}$ is the field magnitude at
$y=0$ and $|\nabla\ln B(y)|_{y=0}\equiv 1/L$ defines a constant gradient
length scale $L$. The equations of motion in the $(x,y)$-plane perpendicular
to the magnetic field are

$\left.\begin{array}[]{rc l}\ddot{x}&=&\Omega_{0}\,(1-y/L)\,\dot{y}\\\
\ddot{y}&=&-\,\Omega_{0}\,(1-y/L)\,\dot{x}\end{array}\right\\},$ (1)

where $\Omega_{0}=qB_{0}/M$ denotes the constant gyrofrequency for a particle
of charge $q$ and mass $M$. By using the canonical momentum conservation law

$P_{x}=M\,\left[\dot{x}\;-\frac{}{}\Omega_{0}\left(y-y^{2}/2L\right)\right]\equiv
M\Omega_{0}L\,\left(u-\frac{1}{2}\right),$ (2)

we can now express the $y$-dynamics in terms of the energy conservation law
$E=M\dot{y}^{2}/2+V(y,u)$, where the potential energy $V(y,u)\equiv
M\dot{x}^{2}/2$ represents a quartic potential in $y$, with minima
$(\dot{x}=0)$ at $y/L=1\pm\sqrt{2u}$ and a local maximum at $y=L$, where
$V(L,u)=\frac{1}{2}\,M\,(u\,\Omega_{0}L)^{2}$.

The study of charged particle motion in the presence of a magnetic neutral
sheet in a magnetized space plasma has been carried out elsewhere Speiser_1965
; Rothwell_Yates_1984 ; Parks_2003 . In the present paper, we expand the
recent work of Kabin Kabin_2021  on explicit solutions of this system and
calculate, for each orbit type, the action-angle coordinates and the
generating function for the canonical transformation to these coordinates.
Because the associated 2D orbits are labeled by two exact invariants (energy
and canonical momentum along the $x$-axis), the motion is exactly integrable
and the construction of the magnetic moment is not required.

The remainder of the paper is organized as follows. In Sec. II, the equations
of motion (1) in the $(x,y)$ plane are transformed into the problem of a
particle moving in a symmetric quartic potential, whose solutions in terms of
the Jacobi elliptic functions and integrals are well known Reichl_1984 ;
Brizard_Westland_2017 . In Sec. III, these orbital solutions are used to
calculate an action integral $J\equiv(1/2\pi)\oint P_{y}\,dy$ for each
particle orbit in terms of elliptic integrals. In addition, orbital periods
are recovered from these action integrals in accordance to standard
definitions Goldstein_2002 . In Sec. IV, an angle coordinate $\theta$ is
introduced for each orbit and these angles are shown to be canonical to the
action coordinates introduced in Sec. III. In Sec. V, the function $S$ that
generates the canonical transformation $(P_{y},y)\rightarrow(J,\theta)$ is
derived from the relation Arnold_1989  $P_{y}\,dy=J\,d\theta+dS$. In Sec. VI,
we investigate the problem of a particle moving in an asymmetric quartic
potential in the presence of a constant electric field. Finally, our work is
summarized in Sec. VII, while useful definitions of the Jacobi epsilon and
zeta functions are presented in App. A.

##  II Particle orbits in a straight magnetic field with constant gradient

In this Section, we derive solutions to all possible orbits associated with
the equations of motion (1). We begin with their dimensionless form:

$\displaystyle x^{\prime\prime}$ $\displaystyle=$
$\displaystyle(1-y)\,y^{\prime},$ (3) $\displaystyle y^{\prime\prime}$
$\displaystyle=$ $\displaystyle-\,(1-y)\,x^{\prime},$ (4)

where the coordinates $(x,y)$ are normalized to the length scale $L$ and time
derivatives (denoted by primes) have been normalized to the dimensionless time
$\tau=\Omega_{0}t$. Equations (3)-(4) have two conservation laws: the total
(kinetic) energy

$\left(x^{\prime}\right)^{2}\;+\;\left(y^{\prime}\right)^{2}\;=\;v^{2},$ (5)

and the canonical momentum along the $x$-axis

$x^{\prime}\;=\;u\;-\;\frac{1}{2}\,\left(1-y\right)^{2},$ (6)

where the normalized perpendicular kinetic energy $v^{2}$ and the normalized
shifted canonical $x$-momentum $u>0$ are constants. Kabin Kabin_2021
considered the action integrals for periodic motion along the $y$-axis
representing a particle moving in the quartic potential

$V(y,u)\;\equiv\;\frac{1}{2}\,(x^{\prime})^{2}\;=\;\frac{1}{8}\left[2u\;-\frac{}{}(1-y)^{2}\right]^{2},$
(7)

with total energy $E=v^{2}/2$. This quartic potential, as shown in Fig. 1, is
composed of two wells with minima at $y=1\pm\sqrt{2u}$, which are separated by
a barrier (local maximum) at $y=1$ with height $u^{2}/2$. By inserting the
momentum conservation law (6) into the equation of motion (4), we obtain the
second-order ordinary differential equation

$y^{\prime\prime}\;=\;-\,u\,(1-y)\;+\;\frac{1}{2}(1-y)^{3},$

which can be transformed into the standard Duffing equation, whose solution in
terms of the Jacobi elliptic functions is well known Nayfeh_1973 ; Lawden_1989
. In the remainder of this section, we will use these elliptic solutions to
solve the problem of finding all possible orbits associated with the coupled
equations of motion (1).

Figure 1: Plot of the quartic potential (7) in the range $-2\leq y\leq 3$ for
$u=1/2$ (solid) and $u=1$ (dashed). The two wells (with minima at at
$y=1\pm\sqrt{2u}$) are separated by a barrier (local maximum) at $y=1$ with
height $u^{2}/2$.

First, using the transformation $s\equiv 1-y$ and the definition $v\equiv
u\,{\sf e}$, Eq. (4) is transformed into the energy conservation law

$\left(s^{\prime}\right)^{2}\;=\;\left(u\,(1+{\sf
e})\;-\;\frac{1}{2}\,s^{2}\right)\left(u\,({\sf
e}-1)\;+\;\frac{1}{2}\,s^{2}\right).$ (8)

The solution for $s(\tau)$ is periodic, with two turning point $\pm s_{1}$,
where $s_{1}=\sqrt{2u\,(1+{\sf e})}$ for all ${\sf e}$. If ${\sf e}<1$,
however, two additional turning points appear at $\pm s_{2}$, where
$s_{2}=\sqrt{2u\,(1-{\sf e})}$. In what follows, the initial conditions used
to integrate Eqs. (3)-(4) are $(x_{0},x_{0}^{\prime})=(0,-\,u\,{\sf e})$ and
$(y_{0},y_{0}^{\prime})=(1-s_{1},0)$, which are consistent with the
conservation laws (5)-(6).

###  II.1 Type I Orbit

When ${\sf e}<1$, the particle is either trapped in the left well
$1-s_{1}<y<1-s_{2}$ or the right well $1+s_{2}<y<1+s_{1}$, i.e., $y$ does not
go through $1$ (Kabin Kabin_2021  calls these type I orbits). We now consider
periodic solutions for $s(\tau)=s_{1}\,P(\nu\tau)$, expressed in terms of
solutions to the Jacobi elliptic differential equation

$\left(P^{\prime}\right)^{2}\;=\;\left(1-P^{2}\right)\left(P^{2}\;-\;m^{\prime}\right),$
(9)

where $0<m^{\prime}=(1-{\sf e})/(1+{\sf e})=s_{2}^{2}/s_{1}^{2}<1$. For this
case Lawden_1989 ; NIST_Jacobi , the solution of Eq. (9) is $P(\zeta)={\rm
dn}(\zeta|m)$, where $\zeta=\nu\tau$ and $\nu=\sqrt{u\,(1+{\sf e})/2}$, with
$m=1-m^{\prime}=2\,{\sf e}/(1+{\sf e})<1$. Hence, the solution for a type-I
orbit is

$y(\tau;u,{\sf e})\;=\;1\;-\;s_{1}\;{\rm dn}(\zeta\,|\,m),$ (10)

where the orbit is trapped in the left well (see Fig. 1): $1-s_{1}\leq y\leq
1-s_{2}$. The particle velocity in this case is

$y^{\prime}(\tau;u,{\sf e})\;=\;2\,u\,{\sf e}\;{\rm cn}(\zeta\,|\,m)\,{\rm
sn}(\zeta\,|\,m),$ (11)

and the period for this motion (normalized to $\Omega_{0}$) is $T=2\,{\sf
K}(m)/\nu$, where

${\sf
K}(m)\;\equiv\;\int_{0}^{\pi/2}\frac{d\varphi}{\sqrt{1-m\,\sin^{2}\varphi}}$

denotes the complete elliptic integral of the first kind footnote . We can now
calculate the orbit-averaged position $\langle
y\rangle\equiv\int_{0}^{T}y(\tau)\,d\tau/T$ using Eq. (10), which can be
expressed as

$\displaystyle\langle y\rangle(u,{\sf e})$ $\displaystyle=$ $\displaystyle
1\;-\;\frac{\nu}{{\sf K}(m)}\int_{0}^{2{\sf K}(m)}{\rm dn}(\zeta|m)\,d\zeta$
(12) $\displaystyle=$ $\displaystyle 1\;-\;\frac{\pi\nu}{{\sf K}(m)}.$

Figure 2 (top) shows the orbit average $\langle y\rangle(u,{\sf e})$ in the
range $0\leq{\sf e}\leq 1$ for $u=0.5$ (A) and 1.0 (B). For all $u>0$, the
orbit average is $\langle y\rangle(u,1)=1$. The dashed curve in the top frame
of Fig. 2 shows the guiding-center approximation for the polarization shift
calculated in Ref. Brizard_2017 . We note that the comparison of predictions
from guiding-center theory Cary_Brizard_2009  with the exact solution is
possible because the particle orbit never crosses the neutral sheet (where the
magnetic field vanishes) at $y=1$.

From the solution (10) for $y(\tau;u,{\sf e})$, we can also obtain an
expression for $x^{\prime}(\tau;u,{\sf e})$ from Eq. (6):

$x^{\prime}(\tau;u,{\sf e})\;=\;u\left[1\;-\frac{}{}(1+{\sf e})\,{\rm
dn}^{2}(\nu\tau|m)\right].$ (13)

First, when averaging the right side of Eq. (13) over an entire orbital period
$T=2{\sf K}/\nu$, we find the orbit-averaged drift velocity

$\displaystyle\left\langle x^{\prime}\right\rangle(u,{\sf e})$
$\displaystyle=$ $\displaystyle u\left[1\;-\;\frac{1+{\sf e}}{2\,{\sf
K}}\int_{0}^{2{\sf K}}{\rm dn}^{2}(\zeta|m)\,d\zeta\right]$ (14)
$\displaystyle=$ $\displaystyle u\left[1\;-\;(1+{\sf e})\;\frac{{\sf
E}(m)}{{\sf K}(m)}\right],$

where ${\sf E}(m)=\int_{0}^{\pi/2}\sqrt{1-m\,\sin^{2}\varphi}\,d\varphi$
denotes the complete elliptic integral of the second kind. The orbit-averaged
drift velocity is shown in Fig. 2 (bottom) in addition to its guiding-center
approximation Brizard_2017  shown as a dashed curve. Next, the explicit
integration of Eq. (13), using the initial condition $x(0)=0$, yields the
expression

$\displaystyle x(\tau;u,{\sf e})$ $\displaystyle=$ $\displaystyle
u\,\tau\;-\;2\nu\int_{0}^{\nu\tau}{\rm dn}^{2}(\zeta|m)\,d\zeta$ (15)
$\displaystyle=$ $\displaystyle u\,\tau\;-\;2\nu\;{\mathcal{E}}(\nu\tau|m)$
$\displaystyle\equiv$ $\displaystyle\left\langle
x^{\prime}\right\rangle\,\tau\;-\;2\nu\;{\mathcal{Z}}(\nu\tau,m),$

where the $2{\sf K}$-periodic Jacobi zeta function ${\mathcal{Z}}(\nu\tau,m)$
is expressed in terms of the quasi-periodic Jacobi epsilon function
${\mathcal{E}}(\nu\tau|m)$ according to the expression NIST_Jacobi :
${\mathcal{Z}}(\zeta,m)={\mathcal{E}}(\zeta|m)-{\sf E}(m)\zeta/{\sf K}(m)$;
see App. A for details about the Jacobi epsilon and zeta functions. A plot of
$x(\tau;u,{\sf e})$ is shown in the top frame of Fig. 3 for $u=1/2$ and ${\sf
e}=0.9$ in the range $0\leq\tau\leq 2{\sf K}/\nu$, with the drift motion
$\langle x^{\prime}\rangle\,\tau$ shown as a dashed line.

Figure 2: (Top) Plot of the orbit-averaged position $\langle y\rangle(u,{\sf
e})$ in the range $0\leq{\sf e}\leq 1$ for $u=0.5$ (A) and 1.0 (B). At the
separatrix $({\sf e}=1)$, the orbit-averaged position is $\langle
y\rangle(u,1)=1$ for all $u>0$. (Bottom) Plot of the normalized orbit-averaged
drift velocity $\langle x^{\prime}\rangle/u$ in the range $0\leq{\sf e}\leq
1$. In both plots, the dashed curves represent the guiding-center
approximations discussed in Ref. Brizard_2017 . Figure 3: (Top) Plot of
$x(\tau;u,{\sf e})$ for $u=1/2$ and ${\sf e}=0.9$ (solid curve), with the
drift motion $\langle x^{\prime}\rangle\,\tau$ shown as a dashed line.
(Bottom) Plot of $\overline{x}(\tau;u,{\sf e})$ for $u=1/2$ and ${\sf e}=1.8$
(solid curve), with the drift motion
$\langle\overline{x}^{\prime}\rangle\,\tau$ shown as a dashed line.

###  II.2 Type II Orbit

When ${\sf e}>1$, the particle is trapped in the well $1-s_{1}<y<1+s_{1}$
above the potential barrier (see Fig. 1) at $y=1$ (Kabin Kabin_2021  calls
these type II orbits). Since $1<m=2{\sf e}/(1+{\sf e})<2$, the orbit solution
can be expressed using the identity $(m>1)$: ${\rm dn}(\zeta|m)={\rm
cn}(\sqrt{m}\,\zeta|m^{-1})\equiv{\rm cn}(\xi|\overline{m})$, so that Eq. (10)
becomes Reichl_1984

$\overline{y}(\tau;u,{\sf e})\;=\;1\;-\;s_{1}\;{\rm
cn}\left(\xi\,|\,\overline{m}\right),$ (16)

where $\xi=\overline{\nu}\tau$, with $\overline{\nu}=\sqrt{m}\,\nu=\sqrt{u{\sf
e}}$, while the particle velocity in this case is

$\overline{y}^{\prime}(\tau;u,{\sf e})\;=\;u\,(1+{\sf e})\;\sqrt{m}\;{\rm
sn}\left(\xi\,|\,\overline{m}\right)\,{\rm
dn}\left(\xi\,|\,\overline{m}\right),$ (17)

and the period for this motion (normalized to $\Omega_{0}$) is
$\overline{T}=4\,{\sf K}(\overline{m})/\overline{\nu}\equiv 4\overline{\sf
K}/\overline{\nu}$. We can now calculate the orbit-averaged position
$\langle\overline{y}\rangle\equiv\int_{0}^{\overline{T}}\overline{y}(\tau)\,d\tau/\overline{T}$
using Eq. (16), which can be expressed as

$\langle\overline{y}\rangle(u,{\sf
e})\;=\;1\;-\;\frac{\sqrt{\overline{m}}}{4\overline{\sf
K}}\int_{0}^{4\overline{\sf K}}{\rm cn}(\xi|\overline{m})\,d\xi\;=\;1.$ (18)

Hence, no polarization shift is observed for a type II orbit, which is
consistent with the breakdown of the guiding-center approximation since the
particle orbit crosses the neutral sheet where the magnetic field vanishes.

From the solution (16) for $\overline{y}(\tau;u,{\sf e})$, we can also obtain
an expression for $\overline{x}^{\prime}(\tau;u,{\sf e})$ from Eq. (6):

$\overline{x}^{\prime}(\tau;u,{\sf e})\;=\;u\left[1\;-\frac{}{}(1+{\sf
e})\,{\rm cn}^{2}(\overline{\nu}\,\tau|\overline{m})\right].$ (19)

When averaging the right side of Eq. (19) over an entire orbital period
$\overline{T}=4\overline{\sf K}/\overline{\nu}$, we find the orbit-averaged
drift velocity

$\displaystyle\left\langle\overline{x}^{\prime}\right\rangle(u,{\sf e})$
$\displaystyle=$ $\displaystyle u\left[1\;-\;\frac{1+{\sf e}}{4\,\overline{\sf
K}}\int_{0}^{4\overline{\sf K}}{\rm cn}^{2}(\xi)\,d\xi\right]$ (20)
$\displaystyle=$ $\displaystyle u\,{\sf e}\left[1\;-\;2\,\frac{{\sf
E}(\overline{m})}{{\sf K}(\overline{m})}\right],$

which is shown in Fig. 4, where a drift reversal is observed when ${\sf
E}(\overline{m})/{\sf K}(\overline{m})>\frac{1}{2}$. A similar reversal is
observed in the bounce-averaged toroidal drift precession frequency for the
banana orbit of a trapped particle in an axisymmetric tokamak magnetic field
Brizard_2011 . Next, the explicit integration of Eq. (19) yields the
expression

$\displaystyle\overline{x}(\tau;u,{\sf e})$ $\displaystyle=$ $\displaystyle
u\,\tau\;-\;\frac{u}{\overline{\nu}}\,(1+{\sf
e})\int_{0}^{\overline{\nu}\tau}{\rm cn}^{2}(\xi|\overline{m})\,d\xi$ (21)
$\displaystyle=$ $\displaystyle
u\,\tau\;-\;2\,\overline{\nu}\left[{\mathcal{E}}(\overline{\nu}\tau|\overline{m})\;-\frac{}{}(1-\overline{m})\,\overline{\nu}\,\tau\right]$
$\displaystyle\equiv$
$\displaystyle\left\langle\overline{x}^{\prime}\right\rangle\,\tau\;-\;2\,\overline{\nu}\;{\mathcal{Z}}(\overline{\nu}\tau,\overline{m}).$

A plot of $\overline{x}(\tau;u,{\sf e})$ is shown in the bottom frame of Fig.
3 for $u=1/2$ and ${\sf e}=1.8$ in the range $0\leq\tau\leq 2\overline{\sf
K}/\overline{\nu}$, with the reversed drift motion
$\langle\overline{x}^{\prime}\rangle\,\tau$ shown as a dashed line.

Figure 4: Plot of the drift velocity $\langle\overline{x}^{\prime}\rangle/u$
in the range $1\leq{\sf e}\leq 2$. Note that the drift velocity starts at
$\langle\overline{x}^{\prime}\rangle(u,1)=u$ and experiences drift reversal
when ${\sf E}(\overline{m})/{\sf K}(\overline{m})>1/2$. Figure 5: Plots of
particle orbits in the $(x,y)$-plane: (I) $y(\tau;u,{\sf e})$ versus
$x(\tau;u,{\sf e})$ for $u=1/2$ and ${\sf e}=0.8$; and (II)
$\overline{y}(\tau;u,{\sf e})$ versus $\overline{x}(\tau;u,{\sf e})$ for
$u=1/2$ and ${\sf e}=1.8$. The drift motion along the $x$-axis for each orbit
is shown by the arrow.

###  II.3 Separatrix Solution

At the interface between these two orbit types is the separatrix solution for
${\sf e}=1$ (i.e., $m=1$):

$y(\tau;u,1)\;=\;1\;-\;2\,\sqrt{u}\;{\rm sech}\left(\sqrt{u}\,\tau\right),$
(22)

which is recovered from both solutions (10) and (16). The $x$-velocity, on the
other hand, is

$x^{\prime}(\tau;u,1)\;=\;u\left[1\;-\frac{}{}2\;{\rm
sech}^{2}(\sqrt{u}\,\tau)\right],$ (23)

which yields the $x$-position

$x(\tau;u,1)\;=\;u\,\tau\;-\;2\sqrt{u}\;\tanh(\sqrt{u}\,\tau).$

Figure 6 shows the phase-space plots $(y,y^{\prime})$ for type I and II for
$u=1/2$, with the separatrix solution $({\sf e}=1)$ shown as a dashed curve.

Figure 6: Phase-space plots $(y,y^{\prime})$ for type I $(0<{\sf e}<1)$ and
type II $({\sf e}>1)$. In both cases, the parameter $u=1/2$ is used and the
separatrix $({\sf e}=1)$ is shown as a dashed curve.

##  III Action Integrals and Orbital Periods

The area enclosed by each phase-space curve in Fig. 6 is related to the action
integral $J(E)=(1/2\pi)\oint y^{\prime}\,dy$, from which the orbital period
can be derived as $T(E)\equiv 2\pi\,dJ(E)/dE$. In our model, the energy is
expressed as

$E\;=\;u^{2}{\sf e}^{2}/2,$ (24)

so that the orbital period $T(u,{\sf e})$ for each orbit type can be
calculated as

$T(u,{\sf e})\;\equiv\;\frac{2\pi}{u^{2}{\sf e}}\;\frac{\partial J(u,{\sf
e})}{\partial{\sf e}}\;\equiv\;\frac{2\pi}{\omega(u,{\sf e})},$ (25)

where $\omega(u,{\sf e})$ denotes the angular frequency for the orbit.

###  III.1 Action Integrals

The action integral for a type I orbit $({\sf e}<1)$ is calculated from the
solution (10) as

$\displaystyle J$ $\displaystyle\equiv$ $\displaystyle\frac{1}{2\pi}\oint
y^{\prime}\,dy\;=\;\frac{1}{2\pi}\int_{0}^{T}y^{\prime 2}\;d\tau$ (26)
$\displaystyle=$ $\displaystyle\frac{2\,\nu^{3}}{\pi}\;m^{2}\int_{0}^{2\,{\sf
K}}{\rm cn}^{2}\zeta\;{\rm sn}^{2}\zeta\;d\zeta,$

where $u\,(1+{\sf e})=2\,\nu^{2}$ and $m=2{\sf e}/(1+{\sf e})<1$. Using
standard properties of the the Jacobi elliptic functions [see Eq. (83) in App.
A], we find

$J(u,{\sf e})\;=\;\frac{4\,\nu^{3}}{3\pi}\left[(2-m)\,{\sf
E}(m)\;-\frac{}{}2\,(1-m)\,{\sf K}(m)\right].$ (27)

As ${\sf e}\rightarrow 1$ $(m\rightarrow 1)$, we find
$J(u,1)=4u^{\frac{3}{2}}/(3\pi)$, while the action $J(u,{\sf
e})\sim(u/2)^{\frac{3}{2}}{\sf e}^{2}$ vanishes as ${\sf e}\rightarrow 0$.

The action integral for a type II orbit (${\sf e}>1$), on the other hand, is
calculated from the solution (16) as

$\displaystyle\overline{J}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2\pi}\oint\overline{y}^{\prime}\,dy\;=\;\frac{1}{2\pi}\int_{0}^{\overline{T}}\overline{y}^{\prime
2}\;d\tau$ (28) $\displaystyle=$
$\displaystyle\frac{2\,\nu^{3}}{\pi}\;\sqrt{m}\int_{0}^{4\,\overline{\sf
K}}{\rm sn}^{2}\xi\;{\rm dn}^{2}\xi\;d\xi.$

Using standard properties of the the Jacobi elliptic functions [see Eq. (85)
in App. A], we find

$\overline{J}(u,{\sf
e})\;=\;\frac{8\,\nu^{3}}{3\pi}\;\sqrt{m}\left[(2-m)\,\overline{\sf
E}\;-\frac{}{}(1-m)\;\overline{\sf K}\right].$ (29)

As ${\sf e}\rightarrow 1$ $(m\rightarrow 1)$, we find
$\overline{J}(u,1)=8u^{\frac{3}{2}}/(3\pi)$, which is double the limit for
$J(u,1)$ since only one separatrix well is considered in the type I case (see
Fig. 6). On the other hand, as ${\sf e}\rightarrow\infty$ $(m\rightarrow 2)$,
we find $\overline{J}(u,{\sf e})\rightarrow 4\,(u\,{\sf e})^{\frac{3}{2}}{\sf
K}(\frac{1}{2})/(3\pi)$.

###  III.2 Orbital Periods

For an orbit of type I, using Eq. (27), we find

$\displaystyle T(u,{\sf e})$ $\displaystyle=$
$\displaystyle\left(\frac{2\sqrt{2}\,u^{\frac{3}{2}}}{3\,u^{2}{\sf
e}}\right)\left[\frac{3\,{\sf e}}{\sqrt{1+{\sf e}}}\;{\sf K}(m)\right]$ (30)
$\displaystyle=$ $\displaystyle\frac{2\,{\sf K}(m)}{\sqrt{u\,(1+{\sf
e})/2}}\;\equiv\;\frac{2\,{\sf K}}{\nu}\;=\;\frac{2\pi}{\omega},$

which is the orbital period observed from the type-I solution (10). For an
orbit of type II, using Eq. (29), we find

$\displaystyle\overline{T}(u,{\sf e})$ $\displaystyle=$
$\displaystyle\left(\frac{4\sqrt{2}\,u^{\frac{3}{2}}}{3\,u^{2}{\sf
e}}\right)\left[3\,\sqrt{\frac{{\sf e}}{2}}\;{\sf K}(\overline{m})\right]$
(31) $\displaystyle=$ $\displaystyle\frac{4\,{\sf
K}(\overline{m})}{\sqrt{u\,{\sf e}}}\;\equiv\;\frac{4\,\overline{\sf
K}}{\overline{\nu}}\;=\;\frac{2\pi}{\overline{\omega}},$

which is the orbital period observed from the type-II solution (16).

##  IV Angle Coordinates and Canonical Condition

The canonical transformation from the canonical coordinates $(y,P_{y})$ to the
action-angle coordinates $(\theta,J)$ satisfies the canonical condition

$\frac{\partial y}{\partial\theta}\,\frac{\partial P_{y}}{\partial
J}\;-\;\frac{\partial y}{\partial J}\,\frac{\partial
P_{y}}{\partial\theta}\;=\;1,$ (32)

where $y(\tau;u,{\sf e})\equiv y[\theta(\tau;u,{\sf e}),J(u,{\sf e})]$ and
$P_{y}(\tau;u,{\sf e})\equiv\partial y(\tau;u,{\sf
e})/\partial\tau=[\partial\theta(\tau;u,{\sf e})/\partial\tau]\,\partial
y(\theta,J)/\partial\theta$. By definition, using the Hamiltonian $H=u^{2}{\sf
e}^{2}/2$, the angle coordinate evolves according to the canonical Hamilton
equation

$\frac{\partial\theta}{\partial\tau}\;\equiv\;\frac{\partial H}{\partial
J}\;=\;u^{2}{\sf e}\;\frac{\partial\sf e}{\partial J}\;=\;\omega(J),$ (33)

which was shown in Eq. (25) to be correct for both orbital types. Since the
Hamiltonian is independent of the canonical angle $\theta$, we obtain
$\partial J/\partial\tau=-\,\partial H/\partial\theta\equiv 0$.

We note that the canonical transformation from the canonical coordinates
$(y,P_{y})$ to the action-angle coordinates $(\theta,J)$ can be inverted,
yielding the canonical condition

$\frac{\partial\theta}{\partial y}\,\frac{\partial J}{\partial
P_{y}}\;-\;\frac{\partial\theta}{\partial P_{y}}\,\frac{\partial J}{\partial
y}\;=\;1,$ (34)

for the transformation from the action-angle coordinates $(\theta,J)$ to the
canonical coordinates $(y,P_{y})$.

###  IV.1 Angle coordinates

For type I orbits, we use the definition $\zeta(\theta,J)=\nu\tau\;=\;2\,{\sf
K}(m)\,\theta/(2\pi)$, which yields the angle coordinate

$\theta\;=\;\pi\zeta/{\sf K}(m)\;\equiv\;\omega(u,{\sf e})\,\tau,$ (35)

with $\partial\theta/\partial\tau=\pi\nu/{\sf K}(m)\equiv\omega(u,{\sf e})$.
Hence,

$\left.\begin{array}[]{rcl}y(\theta,J)&=&1\;-\;2\nu\;{\rm
dn}\,\zeta(\theta,J)\\\ &&\\\ P_{y}(\theta,J)&=&2u\,{\sf e}\,{\rm
cn}\,\zeta(\theta,J)\;{\rm sn}\,\zeta(\theta,J)\end{array}\right\\},$ (36)

where the action dependence appears through $J(u,{\sf e})$.

For type II orbits, on the other hand, we use the definition
$\xi(\overline{\theta},\overline{J})=\overline{\nu}\tau\;=\;4\,\overline{\sf
K}\,\overline{\theta}/(2\pi)$, which yields the angle coordinate

$\overline{\theta}\;=\;\pi\xi/(2\,\overline{\sf
K})\;\equiv\;\overline{\omega}(u,{\sf e})\,\tau,$ (37)

with $\partial\overline{\theta}/\partial\tau=\pi\overline{\nu}/(2\overline{\sf
K})\equiv\overline{\omega}(u,{\sf e})$. Hence,

$\left.\begin{array}[]{rcl}\overline{y}(\overline{\theta},\overline{J})&=&1\;-\;2\nu\;{\rm
cn}\,\xi(\overline{\theta},\overline{J})\\\ &&\\\
\overline{P_{y}}(\overline{\theta},\overline{J})&=&2u\,{\sf
e}\;\sqrt{\overline{m}}\,{\rm sn}\,\xi(\overline{\theta},\overline{J})\;{\rm
dn}\,\xi(\overline{\theta},\overline{J})\end{array}\right\\},$ (38)

where the action dependence appears through $\overline{J}(u,{\sf e})$.

###  IV.2 Canonical condition

####  IV.2.1 Type I orbit

We now verify the canonical condition (32) for orbit type I. First, we note
that the derivatives with respect to the action coordinate must be calculated
more carefully according to the expression

$\displaystyle\left.\frac{\partial}{\partial J}\right|_{\theta}$
$\displaystyle=$ $\displaystyle\;\frac{\partial\sf e}{\partial
J}\left(\left.\frac{\partial}{\partial\sf
e}\right|_{\theta}\;+\;\frac{1}{\Omega_{\sf
e}}\;\frac{\partial}{\partial\tau}\right)$ (39) $\displaystyle=$
$\displaystyle\frac{\omega(u,{\sf e})}{u^{2}\,{\sf
e}}\left(\left.\frac{\partial}{\partial\sf
e}\right|_{\theta}\;+\;\frac{1}{\Omega_{\sf
e}}\;\frac{\partial}{\partial\tau}\right),$

where $\Omega_{\sf e}^{-1}\equiv\partial\tau/\partial{\sf e}|_{\theta}$, with

$d\theta\;=\;(\pi\nu/{\sf K})\,d\tau\;+\;\pi\tau\,d{\sf e}\;\partial(\nu/{\sf
K})/\partial{\sf e}\;=\;0,$

so that we find

$\frac{1}{\Omega_{\sf e}}\;=\;-\;\frac{\tau\,{\sf
K}}{\nu}\frac{\partial}{\partial\sf e}\left(\frac{\nu}{\sf
K}\right)\;=\;\frac{\tau}{2{\sf e}\,(1-{\sf e})}\left(\frac{\sf E}{\sf
K}-(1-{\sf e})\right).$

Second, the angular derivatives can be evaluated as
$\partial/\partial\theta=\omega^{-1}\,\partial/\partial\tau$, so that Eq. (32)
becomes

$\displaystyle\frac{\partial y}{\partial\theta}\,\frac{\partial
P_{y}}{\partial J}-\frac{\partial y}{\partial J}\,\frac{\partial
P_{y}}{\partial\theta}$ $\displaystyle=$
$\displaystyle\omega^{-1}\frac{\partial
y}{\partial\tau}\frac{\omega}{u^{2}{\sf e}}\left(\frac{\partial
P_{y}}{\partial\sf e}+\frac{1}{\Omega_{\sf e}}\frac{\partial
P_{y}}{\partial\tau}\right)$ (40) $\displaystyle-\frac{\omega}{u^{2}{\sf
e}}\left(\frac{\partial y}{\partial\sf e}+\frac{1}{\Omega_{\sf
e}}\frac{\partial y}{\partial\tau}\right)\omega^{-1}\frac{\partial
P_{y}}{\partial\tau}$ $\displaystyle=$ $\displaystyle\frac{\partial
y}{\partial\tau}\frac{\partial P_{y}}{\partial E}-\frac{\partial y}{\partial
E}\frac{\partial P_{y}}{\partial\tau},$

where $\partial/\partial E=(u^{2}{\sf e})^{-1}\partial/\partial{\sf e}$. Here,
using Eq. (36), we find

$\displaystyle\frac{\partial y}{\partial\tau}\;\frac{\partial P_{y}}{\partial
E}$ $\displaystyle=$ $\displaystyle 4\,{\rm sn}^{2}\zeta\;{\rm
cn}^{2}\zeta+m\,(2-m)\;\frac{\partial}{\partial m}\left({\rm
sn}^{2}\zeta\frac{}{}{\rm cn}^{2}\zeta\right),$ $\displaystyle\frac{\partial
y}{\partial E}\;\frac{\partial P_{y}}{\partial\tau}$ $\displaystyle=$
$\displaystyle\left({\rm cn}^{2}\zeta\;-\frac{}{}{\rm sn}^{2}\zeta\right){\rm
dn}^{2}\zeta$ $\displaystyle+(2-m)\left({\rm cn}^{2}\zeta\;-\frac{}{}{\rm
sn}^{2}\zeta\right)\;\frac{\partial\,{\rm dn}^{2}\zeta}{\partial m},$

where we used ${\sf e}=m/(2-m)$ and the $m$-derivatives of ${\rm
pq}\,\zeta\equiv{\rm pq}(\zeta|m)={\rm pq}({\sf K}(m)\theta/\pi|m)$ are
evaluated at constant $\theta$. Next, using the identities ${\rm
cn}^{2}\zeta+{\rm sn}^{2}\zeta=1$ and ${\rm dn}^{2}\zeta+m\,{\rm
sn}^{2}\zeta=1$, the canonical condition (40) becomes

$\frac{\partial y}{\partial\theta}\,\frac{\partial P_{y}}{\partial
J}-\frac{\partial y}{\partial J}\,\frac{\partial
P_{y}}{\partial\theta}\;=\;1,$ (41)

where the final result is obtained by cancellations without evaluating the
$m$-derivatives.

####  IV.2.2 Type II orbit

Next, we verify the canonical condition (32) for orbit type II. Here, we
replace Eq. (39) with the expression

$\displaystyle\left.\frac{\partial}{\partial\overline{J}}\right|_{\theta}$
$\displaystyle=$ $\displaystyle\;\frac{\partial\sf
e}{\partial\overline{J}}\left(\left.\frac{\partial}{\partial\sf
e}\right|_{\overline{\theta}}\;+\;\frac{1}{\overline{\Omega}_{\sf
e}}\;\frac{\partial}{\partial\tau}\right)$ (42) $\displaystyle=$
$\displaystyle\frac{\overline{\omega}(u,{\sf e})}{u^{2}\,{\sf
e}}\left(\left.\frac{\partial}{\partial\sf
e}\right|_{\overline{\theta}}\;+\;\frac{1}{\overline{\Omega}_{\sf
e}}\;\frac{\partial}{\partial\tau}\right),$

where $\overline{\Omega}_{\sf e}^{-1}\equiv\partial\tau/\partial{\sf
e}|_{\overline{\theta}}$, with

$d\overline{\theta}\;=\;(\pi\overline{\nu}/2\overline{\sf
K})\,d\tau\;+\;(\pi\tau/2)\,d{\sf e}\;\partial(\overline{\nu}/\overline{\sf
K})/\partial{\sf e}\;=\;0,$

so that we find

$\frac{1}{\overline{\Omega}_{\sf e}}\;=\;-\;\frac{\tau\,\overline{\sf
K}}{\overline{\nu}}\frac{\partial}{\partial\sf
e}\left(\frac{\overline{\nu}}{\overline{\sf
K}}\right)\;=\;\frac{\tau}{2(1-{\sf e}^{2})}\left(\frac{\overline{\sf
E}}{\overline{\sf K}}-(1-{\sf e})\right).$

Second, the angular derivatives can be evaluated as
$\partial/\partial\overline{\theta}=\overline{\omega}^{-1}\,\partial/\partial\tau$,
so that Eq. (32) becomes

$\displaystyle\frac{\partial\overline{y}}{\partial\overline{\theta}}\,\frac{\partial\overline{P}_{y}}{\partial\overline{J}}-\frac{\partial\overline{y}}{\partial\overline{J}}\,\frac{\partial\overline{P}_{y}}{\partial\overline{\theta}}$
$\displaystyle=$
$\displaystyle\overline{\omega}^{-1}\frac{\partial\overline{y}}{\partial\tau}\frac{\overline{\omega}}{u^{2}{\sf
e}}\left(\frac{\partial\overline{P}_{y}}{\partial\sf
e}+\frac{1}{\overline{\Omega}_{\sf
e}}\frac{\partial\overline{P}_{y}}{\partial\tau}\right)$ (43)
$\displaystyle-\;\frac{\overline{\omega}}{u^{2}{\sf
e}}\left(\frac{\partial\overline{y}}{\partial\sf
e}+\frac{1}{\overline{\Omega}_{\sf
e}}\frac{\partial\overline{y}}{\partial\tau}\right)\;\overline{\omega}^{-1}\frac{\partial\overline{P}_{y}}{\partial\tau}$
$\displaystyle=$
$\displaystyle\frac{\partial\overline{y}}{\partial\tau}\frac{\partial\overline{P}_{y}}{\partial
E}-\frac{\partial\overline{y}}{\partial
E}\frac{\partial\overline{P}_{y}}{\partial\tau}.$

Here, using Eq. (36), we find

$\displaystyle\frac{\partial\overline{y}}{\partial\tau}\;\frac{\partial\overline{P}_{y}}{\partial
E}$ $\displaystyle=$ $\displaystyle(2\overline{m}+1)\,\overline{\rm
sn}^{2}\xi\,\overline{\rm dn}^{2}\xi$
$\displaystyle-\;\overline{m}\,(2\overline{m}-1)\;\frac{\partial}{\partial\overline{m}}\left(\overline{\rm
sn}^{2}\xi\frac{}{}\overline{\rm dn}^{2}\xi\right),$
$\displaystyle\frac{\partial\overline{y}}{\partial
E}\;\frac{\partial\overline{P}_{y}}{\partial\tau}$ $\displaystyle=$
$\displaystyle\overline{\rm cn}^{2}\xi\;\left(\overline{\rm
dn}^{2}\xi\;-\frac{}{}\overline{m}\;\overline{\rm sn}^{2}\xi\right)$
$\displaystyle+\;\overline{m}\,(2\overline{m}-1)\;\left(\overline{\rm
dn}^{2}\xi\;-\frac{}{}\overline{m}\;\overline{\rm
sn}^{2}\xi\right)\;\frac{\partial\,\overline{\rm
cn}^{2}\xi}{\partial\overline{m}},$

where we used ${\sf e}=1/(2\overline{m}-1)$ and the $\overline{m}$-derivatives
of $\overline{\rm pq}\,\xi\equiv{\rm pq}(\xi|\overline{m})={\rm pq}(2{\sf
K}(\overline{m})\overline{\theta}/\pi|\overline{m})$ are evaluated at constant
$\overline{\theta}$. Next, using the identities $\overline{\rm
cn}^{2}\xi+\overline{\rm sn}^{2}\xi=1$ and $\overline{\rm
dn}^{2}\xi+\overline{m}\,\overline{\rm sn}^{2}\xi=1$, the canonical condition
(43) becomes

$\frac{\partial\overline{y}}{\partial\overline{\theta}}\,\frac{\partial\overline{P}_{y}}{\partial\overline{J}}-\frac{\partial\overline{y}}{\partial\overline{J}}\,\frac{\partial\overline{P}_{y}}{\partial\overline{\theta}}\;=\;1,$
(44)

where the final result is obtained by cancellations without evaluating the
$\overline{m}$-derivatives.

##  V Canonical Transformation

The canonical condition (32) proves the existence of a generating function $S$
for the canonical transformation Goldstein_2002
$(y,P_{y})\rightarrow(\theta,J)$. Here, the periodic function $S$ is defined
from the expression $P_{y}\;dy=J\;d\theta+dS$, which leads to two
complementary equations

$\left.\begin{array}[]{rcl}J\;+\;\partial S/\partial\theta&=&P_{y}\;\partial
y/\partial\theta\\\ &&\\\ \partial S/\partial J&=&P_{y}\;\partial y/\partial
J\end{array}\right\\},$ (45)

which are consistent with the canonical condition (32). We note that the
expression $P_{y}\,dy-dS=J\,d\theta$ is also consistent with the inverse
canonical condition (34).

In Eq. (45), the partial derivative $\partial/\partial J$ is understood in the
sense of Eq. (39). The first equation in Eq. (45) yields the action-integral
definition

$J\;=\;\frac{1}{2\pi}\;\int_{0}^{2\pi}P_{y}\,\frac{\partial
y}{\partial\theta}\;d\theta\;=\;\frac{1}{2\pi}\int_{0}^{T}\left(y^{\prime}\right)^{2}d\tau.$

This equation also yields a defining expression for the generating function

$\displaystyle S(\theta,J)$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\theta/\omega}P_{y}\,\frac{\partial
y}{\partial\tau}\;d\tau\;-\;J\,\theta$ (46) $\displaystyle=$
$\displaystyle\int_{0}^{\theta/\omega}\left(\frac{\partial
y}{\partial\tau}\right)^{2}\;d\tau\;-\;J\,\theta,$

where $S(0,J)\equiv 0$ and we will show that $S(\theta,J)$ is explicitly
dependent on the angle coordinate $\theta$, i.e., the orbit average $\langle
S\rangle=0$ vanishes for both orbit types.

###  V.1 Type I Orbit

For an orbit of type I, we use the orbit velocity (11) and Eq. (46) yields

$S(\theta,J)=4\nu^{3}\,m^{2}\int_{0}^{\zeta}{\rm sn}^{2}(u|m)\;{\rm
cn}^{2}(u|m)\,du\;-\;J\,\theta.$ (47)

Next, using the Jacobi zeta function (87), the integral can be evaluated as

$\displaystyle 4\nu^{3}\,m^{2}\int_{0}^{\zeta}{\rm sn}^{2}(u|m)\;{\rm
cn}^{2}(u|m)\,du$ $\displaystyle=$
$\displaystyle\frac{4}{3}\,\nu^{3}\left[(2-m)\,{\mathcal{Z}}(\zeta,m)\;+\;\frac{1}{2}\frac{\partial}{\partial\zeta}{\rm
dn}^{2}(\zeta|m)\right]\;+\;J\,\theta,$

and the canonical generating function for orbit type I is

$\displaystyle S(\theta,J)$ $\displaystyle=$
$\displaystyle\frac{4}{3}\,\nu^{3}\left[(2-m)\,{\mathcal{Z}}(\zeta,m)\;+\;\frac{1}{2}\frac{\partial}{\partial\zeta}{\rm
dn}^{2}(\zeta|m)\right]$ (48) $\displaystyle\equiv$
$\displaystyle\frac{2}{3}\,u\,(1+{\sf e})\,\omega(u,{\sf e})\;\frac{\partial
R(\theta,J)}{\partial\theta},$

where the periodic function

$\displaystyle R(\theta,J)$ $\displaystyle\equiv$
$\displaystyle(2-m)\,\ln\vartheta_{4}(\theta/2,{\sf q}(m))$ (49)
$\displaystyle+\;\frac{1}{2}\,(1-m)^{\frac{1}{2}}\;\frac{\vartheta_{3}^{2}(\theta/2,{\sf
q}(m))}{\vartheta_{4}^{2}(\theta/2,{\sf q}(m))}$

is expressed in terms of the theta functions Lawden_1989 ; NIST_Jacobi

$\left.\begin{array}[]{rcl}\vartheta_{3}(u,{\sf
q})&=&1+2\,\sum_{n=1}^{\infty}{\sf q}^{n^{2}}\,\cos(2nu)\\\ &&\\\
\vartheta_{4}(u,{\sf q})&=&1+2\,\sum_{n=1}^{\infty}(-1)^{n}{\sf
q}^{n^{2}}\,\cos(2nu)\end{array}\right\\},$ (50)

which are defined in terms of the nome

${\sf q}(m)\;=\;\exp[-\pi{\sf K}(1-m)/{\sf K}(m)].$

We note that Eq. (48) automatically yields the vanishing orbit average
$\langle S\rangle=0$.

###  V.2 Type II Orbit

For an orbit of type II, we use the orbit velocity (17) and Eq. (46) yields

$\overline{S}(\overline{\theta},\overline{J})=4\nu^{3}\,\sqrt{m}\int_{0}^{\xi}{\rm
sn}^{2}(u|\overline{m})\;{\rm
dn}^{2}(u|\overline{m})\,du-\overline{J}\,\overline{\theta}.$ (51)

Next, using the Jacobi zeta function (88), the integral can be evaluated as

$\displaystyle 4\nu^{3}\,\sqrt{m}\int_{0}^{\xi}{\rm
sn}^{2}(u|\overline{m})\;{\rm dn}^{2}(u|\overline{m})\,du$ $\displaystyle=$
$\displaystyle\frac{4}{3}\,\nu^{3}\,\sqrt{m}\left[(2-m)\,{\mathcal{Z}}(\xi,\overline{m})\;+\;\frac{1}{2}\frac{\partial}{\partial\xi}{\rm
cn}^{2}(\xi|\overline{m})\right]\;+\;\overline{J}\,\overline{\theta},$

and the canonical generating function for orbit type II is

$\displaystyle\overline{S}(\overline{\theta},\overline{J})$ $\displaystyle=$
$\displaystyle\frac{4}{3}\,\nu^{3}\,\sqrt{m}\left[(2-m)\,{\mathcal{Z}}(\xi,\overline{m})\;+\;\frac{1}{2}\frac{\partial}{\partial\xi}{\rm
cn}^{2}(\xi|\overline{m})\right]$ (52) $\displaystyle\equiv$
$\displaystyle\frac{1}{3}\,u\,(1+{\sf e})\,\overline{\omega}(u,{\sf
e})\;\frac{\partial\overline{R}(\overline{\theta},\overline{J})}{\partial\overline{\theta}},$

where the periodic function

$\displaystyle\overline{R}(\overline{\theta},\overline{J})$
$\displaystyle\equiv$
$\displaystyle(2-m)\,\ln\vartheta_{4}(\overline{\theta},{\sf
q}(\overline{m}))$ (53)
$\displaystyle+\;\frac{1}{2}(m-1)^{\frac{1}{2}}\;\frac{\vartheta_{2}^{2}(\overline{\theta},{\sf
q}(\overline{m}))}{\vartheta_{4}^{2}(\overline{\theta},{\sf
q}(\overline{m}))}$

is expressed in terms of the theta function
$\vartheta_{4}(\overline{\theta},{\sf q}(\overline{m}))$ and the theta
function

$\vartheta_{2}(\overline{\theta},\overline{\sf
q})\;=\;2\,\sum_{n=0}^{\infty}\overline{\sf
q}^{(n+\frac{1}{2})^{2}}\,\cos[(2n+1)\,\overline{\theta}],$ (54)

which are defined in terms of the nome

$\overline{\sf q}\;=\;{\sf q}(\overline{m})\;=\;\exp[-\pi{\sf
K}(1-\overline{m})/{\sf K}(\overline{m})].$

Once again, we note that Eq. (52) automatically yields the vanishing orbit
average $\langle S\rangle=0$.

##  VI Motion in the presence of an electric field

In this Section, we consider the impact of a constant electric field on the
dynamics of a charged particle in a straight magnetic field with constant
perpendicular gradient. The case of a constant electric field ${\bf
E}=E_{0}\,\widehat{\sf x}$ in the ignorable $x$-direction has been studied in
Ref. Parks_2003  and will not be considered here.

Instead, we consider the constant electric field ${\bf E}=E_{0}\,\widehat{\sf
y}$, so that the equations of motion (1) become

$\left.\begin{array}[]{rc l}\ddot{x}&=&\Omega_{0}\,(1-y/L)\,\dot{y}\\\
\ddot{y}&=&qE_{0}/M\;-\;\Omega_{0}\,(1-y/L)\,\dot{x}\end{array}\right\\}.$
(55)

These equations still conserve the canonical momentum (2) while the energy
conservation law

$E\;=\;M(\dot{x}^{2}+\dot{y}^{2})/2-qE_{0}\,(y-1)$

is now expressed in terms of an electrostatic potential that is assumed to
vanish at $y=1$. The normalized equation of motion (4) is now replaced with

$y^{\prime\prime}\;=\;w\;-\;(1-y)\,x^{\prime}\;=\;w\;-\;(1-y)\,u\;+\;\frac{1}{2}\,(1-y)^{3},$
(56)

where $w\equiv qE_{0}/(M\Omega_{0}^{2}L)$ is the normalized $E\times B$
velocity.

###  VI.1 Asymmetric quartic potential

We now set $s=1-y$, so that Eq. (56) becomes

$s^{\prime\prime}\;=\;-\,w\;+\;u\,s\;-\;\frac{1}{2}\,s^{3}\;\equiv\;-\,V^{\prime}(s),$
(57)

where the asymmetric quartic potential

$\displaystyle V(s)$ $\displaystyle=$ $\displaystyle
w\,s\;+\;\frac{1}{2}\left(u-\frac{1}{2}\,s^{2}\right)^{2}$ (58)
$\displaystyle=$
$\displaystyle\frac{1}{2}\,u^{2}\;+\;w\,s\;-\;\frac{1}{2}\,u\,s^{2}\;+\;\frac{1}{8}\,s^{4}$

has extrema points $(s_{1},s_{2},s_{3})$ at

$\left.\begin{array}[]{rcl}s_{1}&=&2\,\sigma\;\cos(\pi/6+\psi/3)\\\
s_{2}&=&2\,\sigma\;\sin(\psi/3)\\\
s_{3}&=&-\,2\,\sigma\;\cos(\pi/6-\psi/3)\end{array}\right\\},$ (59)

where $\sigma^{2}=2\,u/3$ and $\psi=\sin^{-1}(w/\sigma^{3})$. These points are
real when $-\sigma^{3}\leq w\leq\sigma^{3}$, so that the quartic potential
(58) has two minima $(s_{1},s_{3})$ and a local maximum $(s_{2})$, where
$s_{3}<s_{2}<s_{1}$. As shown in Fig. 7 (where $-\sigma^{3}<w<0$), a particle
can now be trapped either in a shallow well (left) or a deep well (right) if
the energy level falls below the barrier height $V_{2}=V(s_{2})$.

Figure 7: Plot of the asymmetric quartic potential (58) in the range $-2\leq
s\leq 2$ for $u=1/2$ and $w=-0.08$ (solid). The symmetric quartic potential,
with $w=0$, is shown as a dashed curve.

The problem of a particle moving in the asymmetric quartic potential (58) was
solved exactly in terms of Jacobi and Weierstrass elliptic functions in Ref.
Brizard_Westland_2017 . There, it was shown that a particle trapped in the
shallow well at a given energy level has exactly the same orbital period as a
particle trapped in the deep well at the same energy level. While it should be
possible to construct action-angle coordinates for the trapped motion in the
shallow well and the deep well separately, as well as for the trapped motion
above the potential barrier, we will not pursue this avenue in the present
paper.

###  VI.2 Electric field perturbation

In previous work, Reichl and Zheng Reichl_1984  considered a symmetric double-
well potential perturbed by a small time-dependent (monochromatic) electric
field represented in dimensionless form as $w=\epsilon\,\cos(\Omega\,\tau)$,
which allows a fraction of particles trapped in one well to migrate over the
top of the potential barrier. The analysis of this perturbed Duffing equation
can also be found in Ref. Nayfeh_1973 .

In this Section, we consider the case of a weak time-independent electric
field $w=\epsilon\,\sigma^{3}$, with $\epsilon\ll 1$, so that the asymmetry in
the quartic potential is weak. This perturbed dynamical problem can be
represented as a Hamiltonian perturbation problem, where, using the action-
angle coordinates $(J,\theta)$, the Hamiltonian is expressed as

$H(J,\theta)\;=\;H_{0}(J)\;+\;\epsilon\,\sigma^{3}\,s(J,\theta),$ (60)

where $H_{0}(J)=\frac{1}{2}\,u^{2}\,{\sf e}^{2}(J)$ and $s(J,\theta)$
describes a solution of the unperturbed (symmetric) quartic potential problem.
Because of the electric-field perturbation, the unperturbed action coordinate
$J$ is no longer an invariant

$\frac{\partial J}{\partial\tau}\;=\;-\;\frac{\partial
H}{\partial\theta}(J,\theta)\;=\;-\;\epsilon\,\sigma^{3}\,\frac{\partial
s}{\partial\theta}(J,\theta),$ (61)

while the angle coordinate evolves according the Hamilton equation

$\frac{\partial\theta}{\partial\tau}\;=\;\frac{\partial H}{\partial
J}(J,\theta)\;=\;\omega_{0}(J)\;+\;\epsilon\,\sigma^{3}\,\frac{\partial
s}{\partial J}(J,\theta),$ (62)

where $\omega_{0}(J)=u^{2}{\sf e}\,\partial{\sf e}/\partial J$ is obtained
from Eq. (25).

According to standard perturbation methods Nayfeh_1973 , we now seek a
canonical transformation $(J,\theta)\rightarrow({\mathcal{J}},\Theta)$, where
the new action-angle coordinates $({\mathcal{J}},\Theta)$ are constructed by
Lie-transform methods (see App. C of Ref. Brizard_2015 )

$\displaystyle{\mathcal{J}}$ $\displaystyle=$ $\displaystyle
J\;+\;\epsilon\,\\{S_{1},J\\}\;+\;\epsilon^{2}\\{S_{2},J\\}$ (63)
$\displaystyle+\;\frac{\epsilon^{2}}{2}\left\\{S_{1},\frac{}{}\\{S_{1},J\\}\right\\}+\cdots,$
$\displaystyle\Theta$ $\displaystyle=$
$\displaystyle\theta\;+\;\epsilon\,\\{S_{1},\theta\\}\;+\;\epsilon^{2}\\{S_{2},\theta\\}$
(64)
$\displaystyle+\;\frac{\epsilon^{2}}{2}\left\\{S_{1},\frac{}{}\\{S_{1},\theta\\}\right\\}+\cdots.$

Here, the generating functions $S_{1}(J,\theta)$ and $S_{2}(J,\theta)$ are
used to remove the angular dependence from the new Hamiltonian up to second
order

$\displaystyle{\mathcal{H}}(J)$ $\displaystyle\equiv$ $\displaystyle
H(J,\theta)\;-\;\epsilon\,\\{S_{1},H\\}\;-\;\epsilon^{2}\,\\{S_{2},H\\}$ (65)
$\displaystyle+\;\frac{\epsilon^{2}}{2}\;\left\\{S_{1},\frac{}{}\\{S_{1},H\\}\right\\}\;+\;\cdots$
$\displaystyle=$
$\displaystyle{\mathcal{H}}_{0}(J)\;+\;\epsilon\,{\mathcal{H}}_{1}(J)\;+\;\epsilon^{2}{\mathcal{H}}_{2}(J)+\cdots,$

and the action-angle Poisson bracket $\\{\;,\;\\}$ is defined as

$\\{f,\;g\\}=\frac{\partial f}{\partial\theta}\;\frac{\partial g}{\partial
J}-\frac{\partial f}{\partial J}\;\frac{\partial
g}{\partial\theta}\equiv\frac{\partial}{\partial J}\left(g\;\frac{\partial
f}{\partial\theta}\right)-\frac{\partial}{\partial\theta}\left(g\;\frac{\partial
f}{\partial J}\right).$

By definition, the zeroth-order Hamiltonian is ${\mathcal{H}}_{0}(J)\equiv
H_{0}(J)$, and the lowest-order (unperturbed) angular velocity is
$\omega_{0}(J)\equiv{\mathcal{H}}_{0}^{\prime}(J)$.

From Eq. (65), the new first-order Hamiltonian is defined by the relation

${\mathcal{H}}_{1}(J)\;=\;\sigma^{3}\,s(J,\theta)\;-\;\omega_{0}(J)\;\frac{\partial
S_{1}}{\partial\theta},$ (66)

and, since the left side of Eq. (66) is independent of the angle $\theta$, we
readily find the solution for the new first-order Hamiltonian

${\mathcal{H}}_{1}(J)\;\equiv\;\sigma^{3}\,\langle
s\rangle\;=\;\frac{\sigma^{3}}{2\pi}\int_{0}^{2\pi}s(J,\theta)\,d\theta,$ (67)

while the first-order generating function is defined as the non-secular
solution $S_{1}(J,\theta)$ obtained from the angle-dependent part of Eq. (66):

$\frac{\partial
S_{1}}{\partial\theta}\;=\;\frac{\sigma^{3}}{\omega_{0}}\left(s\;-\frac{}{}\langle
s\rangle\right)\;\equiv\;\frac{\sigma^{3}}{\omega_{0}}\,\widetilde{s}.$ (68)

At second order, the new Hamiltonian is defined from the expression

${\mathcal{H}}_{2}(J)=-\;\omega_{0}\;\frac{\partial
S_{2}}{\partial\theta}-\omega_{1}\;\frac{\partial
S_{1}}{\partial\theta}-\frac{\sigma^{3}}{2}\left\\{S_{1},\;\widetilde{s}(J,\theta)\right\\},$
(69)

where $\omega_{1}(J)\equiv{\mathcal{H}}_{1}^{\prime}(J)$. Since the left side
of Eq. (69) is $\theta$-independent, the new second-order Hamiltonian is then
defined as

$\displaystyle{\mathcal{H}}_{2}(J)$ $\displaystyle=$
$\displaystyle-\;\frac{\sigma^{3}}{2}\left\langle\left\\{S_{1},\frac{}{}s(J,\theta)\right\\}\right\rangle\equiv-\;\frac{\sigma^{3}}{2}\frac{\partial}{\partial
J}\left\langle\widetilde{s}\;\frac{\partial
S_{1}}{\partial\theta}\right\rangle$ (70) $\displaystyle=$
$\displaystyle-\,\frac{\sigma^{6}}{2}\,\frac{\partial}{\partial
J}\left[\frac{1}{\omega_{0}}\left(\langle s^{2}\rangle\;-\frac{}{}\langle
s\rangle^{2}\right)\right].$

The second-order generating function is defined as the non-secular solution
$S_{2}(J,\theta)$ obtained from the angle-dependent part of Eq. (69).

Once the new Hamiltonian (65) is constructed up to any desired order in
$\epsilon$ (say $\epsilon^{n}$), the new action ${\mathcal{J}}$ is conserved
up the $n$th-order, i.e.,
$d{\mathcal{J}}/d\tau={\mathcal{O}}(\epsilon^{n+1})$, and the angular motion
$d\Theta/d\tau={\mathcal{H}}^{\prime}({\mathcal{J}})\equiv\omega({\mathcal{J}})$
is exactly solved as $\Theta(\tau,{\mathcal{J}})=\omega({\mathcal{J}})\,\tau$.

####  VI.2.1 Type I orbit

We begin our perturbation analysis of the particle motion in an asymmetric
quartic potential with type-I particle orbit, where $s(J,\theta)=2\nu\,{\rm
dn}(\theta\,{\sf K}/\pi|m)$. Using the definition for the quasi-periodic
Jacobi amplitude function NIST_Jacobi

$\int_{0}^{\zeta}{\rm
dn}(u|m)\,du\;=\;\phi(\zeta|m)\;\equiv\;\sin^{-1}\left({\rm
sn}(\zeta|m)\right),$

with $\phi(\zeta+2{\sf K}|m)=\phi(\zeta|m)+\pi$, we find the orbit average

$\langle s\rangle\;=\;\frac{\nu}{\sf K}\int_{0}^{2{\sf K}}{\rm
dn}\,\zeta\;d\zeta\;=\;\frac{\pi\nu}{\sf K}\;\equiv\;\omega_{0}(J),$

so that the new first-order Hamiltonian (67) becomes

${\mathcal{H}}_{1}(J)\;=\;\sigma^{3}\,\omega_{0}(J).$ (71)

The generating function $S_{1}$ can be obtained from Eq. (66) by direct
integration

$\displaystyle S_{1}(J,\theta)$ $\displaystyle=$
$\displaystyle\frac{2\nu\,\sigma^{3}}{\omega_{0}}\left(\int{\rm
dn}(\theta\,{\sf K}/\pi|m)\,d\theta\;-\;\frac{\pi\theta}{2{\sf K}}\right)$
(72) $\displaystyle=$ $\displaystyle\sigma^{3}\left[2\,\phi(\theta\,{\sf
K}/\pi|m)\;-\frac{}{}\theta\right],$

which can be used in the first-order canonical transformation in Eqs.
(63)-(64).

Next, we calculate the new second-order Hamiltonian from Eq. (70), where
$\langle s\rangle=\omega_{0}$ and

$\langle s^{2}\rangle\;=\;\frac{2\nu^{2}}{\sf K}\int_{0}^{2{\sf K}}{\rm
dn}^{2}\zeta\,d\zeta\;=\;4\nu^{2}\,{\sf E}/{\sf K}=(4\nu/\pi)\,\omega_{0}{\sf
E},$

so that Eq. (70) yields

${\mathcal{H}}_{2}(J)\;=\;-\,\frac{2\,\sigma^{6}}{\pi}\,\frac{\partial}{\partial
J}\left(\nu\,{\sf E}\;-\;\frac{\pi}{4}\;\omega_{0}\right),$ (73)

which yields the new Hamiltonian (65) up to second order.

####  VI.2.2 Type II orbit

Next, we consider the perturbation analysis of the particle motion in an
asymmetric quartic potential with type-II particle orbits, where
$s(\overline{J},\overline{\theta})=2\nu\,{\rm
cn}(2\overline{\theta}\,\overline{\sf K}/\pi|\overline{m})$. Using the
integral

$\int_{0}^{\xi}{\rm
cn}(u|\overline{m})\,du\;=\;\sqrt{m}\;\tan^{-1}\left(\sqrt{\overline{m}}\frac{}{}{\rm
sd}(\xi|\overline{m})\right),$

we find $\langle s\rangle=0$ and the new first-order Hamiltonian is
${\mathcal{H}}_{1}(\overline{J})=0$. The first-order generating function is
thus defined by the relation

$\frac{\partial
S_{1}}{\partial\overline{\theta}}\;=\;\frac{2\nu\sigma^{3}}{\overline{\omega}_{0}}\;{\rm
cn}(2\overline{\theta}\,\overline{\sf K}/\pi|\overline{m}),$ (74)

which yields the solution

$\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle\frac{2\nu\sigma^{3}}{\overline{\omega}_{0}}\int_{0}^{\overline{\theta}}{\rm
cn}(2\theta\,\overline{\sf
K}/\pi|\overline{m})\;d\theta=\frac{\pi\nu\sigma^{3}}{\overline{\sf
K}\overline{\omega}_{0}}\int_{0}^{\xi}{\rm cn}(u|\overline{m})\,du$ (75)
$\displaystyle=$ $\displaystyle
2\,\sigma^{3}\tan^{-1}\left(\sqrt{\overline{m}}\frac{}{}{\rm
sd}(\xi|\overline{m})\right),$

which can be used in the first-order canonical transformation in Eqs.
(63)-(64).

Next, we calculate the new second-order Hamiltonian from Eq. (70), where

$\langle s^{2}\rangle\;=\;\frac{\nu^{2}}{\overline{\sf
K}}\int_{0}^{4\overline{\sf K}}{\rm
cn}^{2}\xi\,d\xi\;=\;\frac{4m\nu^{2}}{\overline{\sf K}}\,\left[\overline{\sf
E}\;-\frac{}{}(1-\overline{m})\,\overline{\sf K}\right],$

so that Eq. (70) yields

${\mathcal{H}}_{2}(\overline{J})\;=\;-\,\frac{4\,\sigma^{6}}{\pi}\,\frac{\partial}{\partial\overline{J}}\left[\overline{\nu}\,\left(\overline{\sf
E}\;-\frac{}{}(1-\overline{m})\,\overline{\sf K}\right)\right],$ (76)

which yields the new Hamiltonian (65) up to second order.

####  VI.2.3 Approximate particle orbit

The next step in the perturbative solution of the equations of motion (55) is
to invert the action-angle transformations (63)-(64)

$\displaystyle J(\tau,{\mathcal{J}})$ $\displaystyle=$
$\displaystyle{\mathcal{J}}\;-\;\epsilon\,\\{S_{1},{\mathcal{J}}\\}\;-\;\epsilon^{2}\\{S_{2},{\mathcal{J}}\\}$
(77)
$\displaystyle+\;\frac{\epsilon^{2}}{2}\left\\{S_{1},\frac{}{}\\{S_{1},{\mathcal{J}}\\}\right\\}+\cdots,$
$\displaystyle\theta(\tau,{\mathcal{J}})$ $\displaystyle=$
$\displaystyle\Theta\;-\;\epsilon\,\\{S_{1},\Theta\\}\;-\;\epsilon^{2}\\{S_{2},\Theta\\}$
(78)
$\displaystyle+\;\frac{\epsilon^{2}}{2}\left\\{S_{1},\frac{}{}\\{S_{1},\Theta\\}\right\\}+\cdots,$

where $\Theta(\tau,{\mathcal{J}})=\omega({\mathcal{J}})\,\tau$, with
$\omega({\mathcal{J}})\equiv{\mathcal{H}}^{\prime}({\mathcal{J}})$, is used
after derivatives with respect to $\Theta$ have been taken. Lastly, we insert
Eqs. (77)-(78) into either the solution $y(J,\theta)$ for type I or II orbits,
given by Eq. (36) or (38), respectively. This perturbative solution was shown
in Ref. Brizard_2015  to be valid up to time scales of order $\epsilon^{-1}$
for the case of the perturbed harmonic oscillator.

##  VII Summary

In this paper, we presented in Sec. II all orbital solutions for the equations
of motion (1) describing the motion of a charged particle in a straight
magnetic field with constant gradient. In addition, we derived expressions for
the action-angle coordinates for each orbit type in Secs. III-IV and obtained
an explicit expression for the generating function for each canonical
transformation in Sec. V. This work follows in the steps of previous work on
the trapped and passing particle orbits in simple tokamak geometry
Brizard_Duthoit_2014 , based on the canonical transformation for the pendulum
problem Brizard_2013 .

Next, we considered the case of motion in an asymmetric quartic potential due
to the presence of a constant electric field in Sec. VI, which is solved
perturbatively for the case of a weak electric field. In future work, we will
construct the canonical transformation to action-angle coordinates based on
the exact solution in terms of elliptic functions presented in Ref.
Brizard_Westland_2017 .

##  Appendix A Integrals of Jacobi Elliptic Functions

In this Appendix, we present a brief review of the Jacobi epsilon and zeta
functions Lawden_1989 ; NIST_Jacobi .

###  A.1 Jacobi Epsilon Function

In calculating action integrals in Sec. III, we need to compute integrals of
even powers of ${\rm sn}(u|m)$:

$I_{2n}(z|m)\;\equiv\;\int_{0}^{z}{\rm sn}^{2n}(u|m)\,du$

for $n\geq 0$. First, $I_{0}(z|m)=z$ and

$I_{2}(z|m)\;=\;\int_{0}^{z}{\rm
sn}^{2}(u|m)du\;=\;m^{-1}\left[z\;-\;{\mathcal{E}}(z|m)\right],$ (79)

where ${\mathcal{E}}(z|m)$ is the Jacobi epsilon function, which is quasi-
periodic: ${\mathcal{E}}(z+2{\sf K}|m)={\mathcal{E}}(z|m)+2\,{\sf E}(m)$.
Next, we use the recurrence relation $(n\geq 1)$ Lawden_1989

$\displaystyle m\,(1+2n)\,I_{2n+2}(z)$ $\displaystyle=$
$\displaystyle(1+m)\,2n\;I_{2n}(z)$ $\displaystyle-\;(2n-1)\,I_{2n-2}(z)$
$\displaystyle+\;{\rm sn}^{2n-1}(z)\,{\rm cn}(z)\,{\rm dn}(z),$

which, for $n=1$, yields

$\displaystyle 3\,m\,I_{4}(z|m)$ $\displaystyle=$ $\displaystyle
2(1+m)\,I_{2}(z|m)\;-\;I_{0}(z|m)$ $\displaystyle+\;{\rm sn}(z|m)\,{\rm
cn}(z|m)\,{\rm dn}(z|m),$

or

$\displaystyle m\int_{0}^{z}{\rm sn}^{4}(u|m)du$ $\displaystyle=$
$\displaystyle\frac{1}{3}\,\left({\rm sn}(z|m)\,{\rm cn}(z|m)\,{\rm
dn}(z|m)\;-\frac{}{}z\right)$ (81)
$\displaystyle+\;\frac{2}{3m}(1+m)\left(z\;-\frac{}{}{\mathcal{E}}(z|m)\right).$

First, using the identity ${\rm cn}^{2}(z|m)=1-{\rm sn}^{2}(z|m)$, we find

$\displaystyle m^{2}\int_{0}^{\zeta}{\rm sn}^{2}u\;{\rm cn}^{2}u\,du$
$\displaystyle=$ $\displaystyle m^{2}\int_{0}^{\zeta}\left({\rm
sn}^{2}u\;-\frac{}{}{\rm sn}^{4}u\right)du$ (82) $\displaystyle=$
$\displaystyle m^{2}\left[I_{2}(\zeta|m)\;-\frac{}{}I_{4}(\zeta|m)\right]$
$\displaystyle=$ $\displaystyle
m\left(\zeta\;-\frac{}{}{\mathcal{E}}(\zeta)\right)$
$\displaystyle-\;\frac{m}{3}\,\left({\rm sn}\,\zeta\;{\rm cn}\,\zeta\;{\rm
dn}\,\zeta\;-\frac{}{}\zeta\right)$
$\displaystyle-\;\frac{2}{3}\,(1+m)\left(\zeta\;-\frac{}{}{\mathcal{E}}(\zeta)\right)$
$\displaystyle=$
$\displaystyle\frac{1}{3}\left[(2-m)\,{\mathcal{E}}(\zeta)\;-\frac{}{}2\,(1-m)\,\zeta\right]$
$\displaystyle-\;\frac{m}{3}\;{\rm sn}\,\zeta\;{\rm cn}\,\zeta\;{\rm
dn}\,\zeta.$

When evaluated at $\zeta=2\,{\sf K}(m)$, this expression yields

$\displaystyle m^{2}\int_{0}^{2{\sf K}}{\rm sn}^{2}u\,{\rm cn}^{2}u\,du$
$\displaystyle=$ $\displaystyle\frac{2}{3}(2-m)\,{\sf E}(m)$ (83)
$\displaystyle-\;\frac{4}{3}\,(1-m)\,{\sf K}(m),$

where we used ${\mathcal{E}}(2{\sf K})=2\,{\sf E}$.

Second, using the identity ${\rm
dn}^{2}(\xi|\overline{m})=1-\overline{m}\,{\rm sn}^{2}(\xi|\overline{m})$, we
find

$\displaystyle\int_{0}^{\xi}\overline{\rm sn}^{2}u\;\overline{\rm
dn}^{2}u\,du$ $\displaystyle=$ $\displaystyle\int_{0}^{\xi}\left(\overline{\rm
sn}^{2}u\;-\;\overline{m}\int_{0}^{\xi}\overline{\rm sn}^{4}u\right)du$ (84)
$\displaystyle=$ $\displaystyle
I_{2}(\xi|\overline{m})\;-\;\overline{m}\;I_{4}(\xi|\overline{m})$
$\displaystyle=$ $\displaystyle
m\left(\xi\;-\frac{}{}\overline{\mathcal{E}}(\xi)\right)$
$\displaystyle-\frac{1}{3}\left(\overline{\rm sn}\,\xi\;\overline{\rm
cn}\,\xi\;\overline{\rm dn}\,\xi\;-\frac{}{}\xi\right)$
$\displaystyle-\frac{2}{3}(1+m)\left(\xi\;-\frac{}{}\overline{\mathcal{E}}(\xi)\right)$
$\displaystyle=$
$\displaystyle\frac{1}{3}\left[(2-m)\,\overline{\mathcal{E}}(\xi)\;-\frac{}{}(1-m)\,\xi\right]$
$\displaystyle-\;\frac{1}{3}\;\overline{\rm sn}\,\xi\;\overline{\rm
cn}\,\xi\;\overline{\rm dn}\,\xi,$

where $\overline{m}\equiv m^{-1}$. When evaluated at $\xi=4\,\overline{\sf
K}$, this expression yields

$\int_{0}^{4\overline{\sf K}}\overline{\rm sn}^{2}u\;\overline{\rm
dn}^{2}u\,du\;=\;\frac{4}{3}\left[(2-m)\,\overline{\sf
E}\;-\frac{}{}(1-m)\,\overline{\sf K}\right],$ (85)

where we used $\overline{\mathcal{E}}(4\overline{\sf K})=4\,\overline{\sf E}$.

###  A.2 Jacobi Zeta and Theta Functions

We now introduce the Jacobi zeta function Lawden_1989 ; NIST_Jacobi

$\displaystyle{\mathcal{Z}}(z,m)$ $\displaystyle\equiv$
$\displaystyle{\mathcal{Z}}\left(\phi(z|m)\;|\;m\right)$ (86) $\displaystyle=$
$\displaystyle{\mathcal{E}}(z|m)\;-\;\frac{{\sf E}(m)}{{\sf K}(m)}\;z,$

where the Jacobi amplitude function $\phi(z|m)$ is defined by the relation
$\sin\phi(z|m)={\rm sn}(z|m)$. If we substitute the definition (86) into Eq.
(82), we obtain

$\displaystyle 3\,m^{2}\int_{0}^{\zeta}{\rm sn}^{2}u\;{\rm cn}^{2}u\,du$
$\displaystyle=$
$\displaystyle\frac{3\,J\,\theta}{4\,\nu^{3}}+(2-m)\,{\mathcal{Z}}(\zeta,m)$
(87) $\displaystyle+\;\frac{1}{2}\,\frac{\partial}{\partial\zeta}{\rm
dn}^{2}(\zeta|m),$

where we used the definition (27) and

$-\;m\;{\rm sn}\,\zeta\;{\rm cn}\,\zeta\;{\rm
dn}\,\zeta\;=\;\frac{1}{2}\frac{\partial}{\partial\zeta}{\rm dn}^{2}\zeta.$

If we substitute the definition (86) into Eq. (84), on the other hand, we
obtain

$\displaystyle\int_{0}^{\xi}\overline{\rm sn}^{2}u\;\overline{\rm
dn}^{2}u\,du$ $\displaystyle=$
$\displaystyle\frac{\overline{J}\,\overline{\theta}}{4\,\nu^{3}\sqrt{m}}+\frac{1}{3}(2-m)\,{\mathcal{Z}}(\xi,\overline{m})$
(88) $\displaystyle+\;\frac{1}{6}\,\frac{\partial}{\partial\xi}\overline{\rm
cn}^{2}(\xi|\overline{m}),$

where we used the definition (29) and we used

$-\;\overline{\rm sn}\,\xi\;\overline{\rm cn}\,\xi\;\overline{\rm
dn}\,\xi\;=\;\frac{1}{2}\frac{\partial}{\partial\xi}\overline{\rm cn}^{2}\xi.$

We note that Eqs. (87) and (88) may be expressed as partial derivatives of
periodic theta functions Lawden_1989  if we express the Jacobi zeta function
as

${\mathcal{Z}}(z,m)\;=\;\frac{\partial}{\partial
z}\ln\vartheta_{4}\left(\frac{\pi z}{2{\sf K}(m)},{\sf q}(m)\right),$ (89)

and the Jacobi elliptic functions

$\displaystyle{\rm dn}^{2}(z|m)$ $\displaystyle=$
$\displaystyle\sqrt{1-m}\;\frac{\vartheta_{3}^{2}(\pi z/2{\sf K},{\sf
q})}{\vartheta_{4}^{2}(\pi z/2{\sf K},{\sf q})},$ (90) $\displaystyle{\rm
cn}^{2}(z|\overline{m})$ $\displaystyle=$
$\displaystyle\sqrt{m-1}\;\frac{\vartheta_{2}^{2}(\pi z/2\overline{\sf
K},\overline{\sf q})}{\vartheta_{4}^{2}(\pi z/2\overline{\sf K},\overline{\sf
q})},$ (91)

where $\vartheta_{2}(u,\overline{\sf q})=2\sum_{n=0}^{\infty}\overline{\sf
q}^{(n+\frac{1}{2})^{2}}\,\cos[(2n+1)u]$ denotes another Jacobi theta
function.

## References

  * (1) A.J. Brizard, Phys. Plasmas 24, 042115 (2017).
  * (2) J.R. Cary and A.J. Brizard, Rev. Mod. Phys. 81, 693 (2009).
  * (3) In the present work, we primarily use the notation (also used in the Mathematica software) found in L.M. Milne-Thomson, Jacobian Elliptic Functions and Theta Functions (Chap. 16), in Handbook of Mathematical Functions, M. Abramowitz and I.A. Stegun, eds., (Dover, 1972).
  * (4) T.W. Speiser, J. Geophys. Res. 70, 4219 (1965).
  * (5) P.L. Rothwell and G.K. Yates, Global single ion effects within the Earth’s plasma sheet, in Magnetic Reconnection in Space and Laboratory Plasmas, Proc. Chapman Conf., Geophysical Monograph Series 30, 51 (American Geophysical Union, 1984). doi: 10.1029/GM030p0051.
  * (6) G.K. Parks, Physics of Space Plasmas, 2nd ed. (Westview, 2003), Chap. 4.
  * (7) K. Kabin, Phys. Plasmas 28, 122101 (2021).
  * (8) L.E. Reichl and W.M. Zheng, Phys. Rev. A 29, 2186 (1984); 30, 1068 (1984).
  * (9) A.J. Brizard and M.C. Westland, Commun. Nonlinear Sci. Numer. Simulat. 43, 351 (2017).
  * (10) H. Goldstein, C. Poole, and J. Safko, Classical Mechanics, 3rd ed. (Addison Wesley, 2002).
  * (11) V.I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed. (Springer-Verlag, 1989).
  * (12) A.H. Nayfeh, Perturbation Methods (Wiley, 1973).
  * (13) D.F. Lawden, Elliptic Functions and Applications (Springer-Verlag, 1989).
  * (14) W. P. Reinhardt and P. L. Walker, Jacobian Elliptic Functions (Chap. 22), in NIST Handbook of Mathematical Functions (Cambridge University Press, 2010).
  * (15) A.J. Brizard, Phys. Plasmas 18, 022508 (2011).
  * (16) A.J. Brizard, An Introduction to Lagrangian Mechanics, 2nd ed. (World Scientific, 2015).
  * (17) A.J. Brizard and F.-X. Duthoit, Phys. Plasmas 21, 052509 (2014).
  * (18) A.J. Brizard, Comm. Nonlinear Sci. Numer. Simulat. 18, 511 (2013).


# Generation of Radio Frequency Radiation by Femtosecond Filaments

Travis Garrett  Jennifer Elle  Michael White  Remington Reid  Air Force
Research Laboratory, Directed Energy Directorate, Albuquerque, NM 87123, USA
Alexander Englesbe  Naval Research Laboratory, Plasma Physics Division,
Washington, DC 20375, USA  Ryan Phillips  Peter Mardahl  Erin Thornton  James
Wymer  Air Force Research Laboratory, Directed Energy Directorate,
Albuquerque, NM 87123, USA  Anna Janicek  Oliver Sale  Leidos Innovations
Center, Albuquerque, NM 87106, USA  Andreas Schmitt-Sody  Air Force Research
Laboratory, Directed Energy Directorate, Albuquerque, NM 87123, USA

###### Abstract

Recent experiments have shown that femtosecond filamentation plasmas generate
ultra-broadband radio frequency radiation (RF). We show that a combination of
plasma dynamics is responsible for the RF: a plasma wake field develops behind
the laser pulse, and this wake excites (and copropagates with) a surface wave
on the plasma column. The surface wave proceeds to detach from the end of the
plasma and propagates forward as the RF pulse. We have developed a four stage
model of these plasma wake surface waves and find that it accurately predicts
the RF from a wide range of experiments, including both 800 nm and 3.9 $\mu$m
laser systems.

The development of chirped pulse amplification Strickland and Mourou (1985)
has enabled an array of exotic physics Faure  _et al._ (2004); Liu  _et al._
(1997); Di Piazza  _et al._ (2012), including atmospheric filamentation Braun
_et al._ (1995). Filamentation in turn produces interesting and useful
effects, from supercontinuum generation Brodeur and Chin (1999); Kandidov  _et
al._ (2003), broadband THz radiation Amico  _et al._ (2008); Andreeva  _et
al._ (2016), and RF pulses, which have recently been explored in great detail
Forestier  _et al._ (2010); Englesbe  _et al._ (2018); Janicek  _et al._
(2020); Mitrofanov  _et al._ (2021a); Englesbe  _et al._ (2021); Mitrofanov
_et al._ (2021b). We have discovered that a novel combination of a plasma wake
field and a surface wave is responsible for this filament RF. Femtosecond
lasers are becoming important spectroscopic tools (both in the lab Kampfrath
_et al._ (2013) and for remote sensing Labutin  _et al._ (2016)), and the
added production of pulsed RF promises to expand the utility of these systems
Park and Field (2016); Doll  _et al._ (2013); Capmany and Novak (2007), and to
provide insight into the physics of filamentation.

During typical $\lambda=$ 800 nm filamentation Kerr self-focusing and plasma
defocusing are in rough balance, resulting a plasma column with a density of
$n_{e}\sim 10^{22}$ m-3 Couairon and Mysyrowicz (2007); Chin (2010). This
electron density corresponds to a plasma frequency of $\omega_{pl}\simeq
6\times 10^{12}$ rad/s, which is closely linked to the generation of THz
radiation. In the single color THz theory of Amico  _et al._ (2008) the
freshly ionized electrons receive a ponderomotive push from the laser pulse,
thus exciting a coherent longitudinal current $I_{z}$ in the plasma (which is
quickly damped by the high collision frequency). This short current pulse
translates at $c$ behind the laser pulse, thereby producing a conical shell of
radially polarized THz radiation.

The RF resembles the THz radiation in some ways, and differs critically in
others. It has a similar conical shell spatial profile and radial
polarization, which indicates that a GHz scale longitudinal current pulse is
also generated within the plasma. However it also has a well defined broadband
peak in the 5-15 GHz range (and is not the tail of the THz radiation Englesbe
_et al._ (2018, 2021); Mitrofanov  _et al._ (2021a)) and the RF amplitude
grows strongly with decreasing pressure Englesbe  _et al._ (2018); Mitrofanov
_et al._ (2021a), in contrast to the THz which falls off sharply Rodriguez and
Dakovski (2010). A distinct physical mechanism thus must be responsible for
the RF.

Exploratory Particle In Cell (PIC) simulations Birdsall and Langdon (2004);
Peterkin and Luginsland (2002); Taflove and Hagness (2005); Yee (1966); Boris
(1970); Villasenor and Buneman (1992); Berenger (1994) (see section $1$ in sup
for more details) provided the key insights that explain the RF generation.
They revealed that a hot shell of electrons expands off the plasma surface
into the surrounding atmosphere over roughly 50 ps (as anticipated by Zhou
_et al._ (2011)), and that the corresponding current density $J_{r}$ depends
strongly on the electron-neutral collision frequency. Subsequent simulations
demonstrated that this Plasma Wake field (PW) excites a Surface Wave (SW) on
the outer boundary of the plasma (namely, a long wavelength Surface Plasmon
Polariton (SPP)Pitarke  _et al._ (2006); Maier  _et al._ (2006), with a speed
approaching $c$). Finally the SW was found to effectively detach from the end
of the plasma, with most of the energy being then converted into a broadband
RF pulse.

The PWSW numerical model is split into 4 stages, as length scales spanning 7
orders of magnitude need to be resolved. First the initial electron
distribution function $\mathcal{N}_{e}$ is determined. In a typical
$\lambda=800$ nm experiment a linearly polarized, single color, 40 mJ pulse is
focused into a quartz cell with a ${f}/60$ mirror, which sets the beam waist
$\omega_{0}$ to 30 $\mu$m (10 mJ and ${f}/15$ focusing are common for
$\lambda=$ 3.9 $\mu$m). The resulting Keldysh parameter
$\gamma=\omega_{laser}\sqrt{2m_{e}I_{pot}}/(q_{e}E)$ is close to 1, so we use
the general ionization rate $\scriptstyle\mathcal{W}$ from Popruzhenko  _et
al._ (2008), which is supported by experimental data Sharma  _et al._ (2018)
near atmospheric pressure (other ionization mechanisms will likely be needed
for 10 $\mu$m lasers Tochitsky  _et al._ ). A laser intensity of $I\sim
5\times 10^{17}$ W/m2 generates enough plasma at 1 atm for intensity clamping
Couairon and Mysyrowicz (2007) to occur (with $\scriptstyle\mathcal{W}$ $\sim$
$10^{12}$ s-1), which leads to an estimated plasma length $L_{pl}$ of $25$ cm,
and a radius $r_{pl}$ of $\sim$ $1$ mm (see section 2 of sup  for more
details). At 1 Torr $I$ increases to $\sim 2\times 10^{18}$ W/m2, which
ionizes much of the O2, and the plasma is reduced to $L_{pl}\simeq 13$ cm and
$r_{pl}\simeq 0.5$ mm. In turn, the lower energy 3.9 $\mu$m laser plasmas are
shorter, with $L_{pl}$ ranging from about $1$ to $2$ cm.

These estimates match well (within 25%) with photographic measurements of the
plasma taken at different pressures (available in section 2 of sup ).
Simulations of filamentation with these parameters have also been performed
using a Unidirectional Pulse Propagation Equations (UPPE) solver Kolesik and
Moloney (2004). These agree with the overall dimensions of the plasma column,
but also indicate that there is rich fine scale structure throughout the
plasma (as expected as we are in the multifilamentary regime Bergé  _et al._
(2004); Ettoumi  _et al._ (2015a, b)). For this letter we simplify and do not
include fine scale plasma perturbations. The plasma column is approximated as
having a constant density central core due to intensity clamping, and then a
100 $\mu$m radial ramp down to zero density at the outer radius $r_{pl}$.

The velocity distribution is determined by our PIC code: given $E(\vec{x},t)$
the air is ionized at the $\scriptstyle\mathcal{W}$ rate Popruzhenko  _et al._
(2008) and the new electrons are accelerated through the remainder of the
pulse Boris (1970) (code to perform these calculations is included in sup ).
In general the velocity distribution is influenced by both the details of
strong field ionization near $\gamma=1$ (e.g. Corkum (1993)), and by the
deformation of the laser pulse during filamentation. For this work we further
simplify and assume that the electrons are ionized with zero initial velocity
and are then accelerated by the remainder of a Gaussian pulse (with $\hat{x}$
polarization and propagating in $+z$). As a whole the initial
$\mathcal{N}_{e}$ is highly non-Maxwellian, and at 100 Torr has a peak kinetic
energy $K_{tail}\simeq$ 5 eV, and an average kinetic energy $K_{avg}\simeq
0.6$ eV, while these increase at 1 Torr to $K_{tail}\simeq 16$ eV and
$K_{avg}\simeq 2$ eV. For 3.9 $\mu$m lasers the kinetic energy is $\sim 25$
times larger, as the laser intensities are comparable and the energy scales as
$\lambda^{2}$.

We next consider the evolution of the plasma column. We construct thin
transverse slices of plasma given $\mathcal{N}_{e}$, with periodic boundary
conditions used in the longitudinal direction $\hat{z}$ (which is valid to
leading order as the electron velocities are a small fraction of $c$) and
simulate the radial evolution with our PIC code. The Debye length is quite
small: $\lambda_{Debye}\simeq 10$ nm, so we use an energy conserving method
Pointon (2008) to calculate the Lorentz force. The electron-neutral elastic
collision frequency $\nu_{eN}$ depends on the cross sections for O2 and N2,
which for our energies is roughly 10 Å2 Itikawa (2006). In turn the electron-
ion momentum-transfer collision frequency is given by
$\nu_{ei}=7.7{\times}10^{-12}n_{e}\ln(\Lambda_{C})/K_{eV}^{3/2}$, where
$\Lambda_{C}=6\pi n_{e}\lambda_{Debye}^{3}$ Thorne and Blandford (2017). The
resulting radial current density $J_{r}$ and electron density $n_{e}$ are then
recorded as functions of radius and time (more detail can be found in section
3 of sup ). These are well resolved with a mesh resolution of $\Delta x=\Delta
y=$ 2 $\mu$m, and a macro-particle weight of ${\sim}10$ on the outer edge of
the plasma.

Figure 1: Log plots of the electron density $n_{e}(r,t)$ showing the plasma
wake variation at pressures of 100 (a), 10 (b), and 1 Torr (c) for a
$\lambda=800$ nm laser pulse. The plasma outer radius has been set to
$r_{pl}=0.5$ mm for all 3 pressures for comparison purposes. Figure 2: Log
plots of the radial current density $J_{r}(r,t)$ at 1 Torr (800 nm source), as
it evolves in a transverse slice of plasma. The initial density at the edge of
the plasma has been given a step function profile, with a value of
$n_{e}=10^{19}$ m-3 in the top plot and $n_{e}=10^{20}$ m-3 on the bottom. The
profiles of the plasma wakes outside the column ($r>r_{pl}=0.5$ mm) are quite
similar, while at the surface SPPs are excited at roughly ${{\sim}}20$ and
${{\sim}}60$ GHz respectively, along with higher frequency volume plasmons in
the interior.

The electron number density for $\lambda=800$ nm simulations of the PW at 100,
10, and 1 Torr can be seen in Fig. 1. The outer edge of the plasma at $t=0$
has a simplified step function profile, with $n_{e}=10^{20}$ m-3 at radius
$r_{pl}=0.5$ mm. Thus, in addition to the PW being launched off of the edge of
this plasma, a coherent radial plasma frequency oscillation of ${\sim}90$ GHz
is stimulated in the interior Dawson (1959) and a SPP at ${\sim}63$ GHz is
excited on the surface Pitarke  _et al._ (2006); Maier  _et al._ (2006); Sakai
_et al._ (2018). The PW that expands out into the neutral atmosphere
($r>r_{pl}$) is insensitive to the density of the outer edge of the plasma
(see Fig. 2), as opposed to the excited surface and volume plasmons.

Analytic approximations provide a useful complement to the simulations.
Consider the late time electron density at a small distance $r_{\Delta}$ off
the surface of the plasma at radius $r_{pl}$. Using Gauss’s law and
approximating $\ln((r_{pl}+r_{\Delta})/r_{pl})$ as $r_{\Delta}/r_{pl}$ for
$r_{\Delta}\ll r_{pl}$ we find the that the number of electrons per unit
length $n_{L}$ that escape to $r_{\Delta}$ scales as:

$n_{L}(r_{\Delta})\simeq\frac{2\pi\epsilon_{0}K_{eV}}{q_{e}}\frac{r_{pl}}{r_{\Delta}},$
(1)

where $K_{eV}\simeq K_{tail}$ is the kinetic energy of the hot electrons that
lead the escape. The radial electrostatic field $E_{r,stat}$ that corresponds
to $n_{L}$ has a simple form: $E_{r,stat}(r_{\Delta})=K_{eV}/r_{\Delta}$. At
$\lambda=800$ nm and 1 Torr $K_{tail}\simeq 16$ eV and $r_{\Delta}\simeq$ 60
$\mu$m, giving $E_{r,stat}\simeq 2.7\times 10^{5}$ V/m. This allows for an
estimate of the PW evolution timescale: an electron with energy $K_{tail}=16$
eV in an electric field of this magnitude follows a parabolic trajectory over
100 ps. This compares well with the spread of electrons seen in Fig. 1. The
magnitude of the radial current density can likewise be approximated as:

$J_{r}\simeq v_{eff}(t)\frac{q_{e}n_{L}}{2\pi
r_{pl}r_{\Delta}};v_{eff}(t)=\min\begin{cases}v_{tail}\\\
v_{diff}(t)\end{cases}$ (2)

where the effective velocity $v_{eff}(t)$ is the minimum of the tail velocity
and the Fickian diffusion speed $v_{diff}(t)=(K_{tail}/(m_{e}\nu t))^{1/2}$.
At 1 Torr this gives $J_{r}\simeq 1\times 10^{5}$ A/m2. Finally we note that
caution is needed for small $r_{\Delta}$ values: Eq.s (1) and (2) only hold up
to densities that are comparable to the original plasma edge density.

We next switch from these transverse PIC simulations to a continuum FDTD-Drude
model in a 2D axisymmetric coordinate system Taflove and Hagness (2005) that
spans the length of the plasma. The Drude model:
$\partial_{t}\vec{J}_{pl}+\nu_{pl}\vec{J}_{pl}=\epsilon_{0}\omega_{pl}^{2}\vec{E}$
is integrated into the code via an auxiliary differential equation Okoniewski
_et al._ (1997). A mixed resolution of $\Delta r=$ 2 $\mu$m and $\Delta z=$ 50
$\mu$m suffices to resolve the SWs. The PW current profile $J_{r}$ is driven
across the surface of the plasma column in the longitudinal direction at the
speed of light. In this work only the external radial currents $J_{r>r_{pl}}$
are used to drive the SW and subsequent RF (see section 4 of sup  for more
detail). The higher frequency volume plasmons can also excite waves in
neighboring lower density plasma, but we suspect that including fine scale
plasma structure will act as surface roughness Kretschmann (1972) and cause
these to continuously detach.

As expected the PW current excites a broadband surface wave with similar
frequency content. The frequency $\omega$ of the wave is considerably lower
than the limit SPP frequency and it thus travels at velocity
$c((2\omega^{2}-\omega_{pl}^{2})/(\omega^{2}-\omega_{pl}^{2}))^{1/2}$ which
approaches the speed of light. The resulting phase-matched copropagation leads
to a steady growth in the surface wave intensity (see also Smith and Purcell
(1953); Dahan  _et al._ (2020)).

The SW is well approximated by the Sommerfeld-Goubau Sommerfeld (1899); Goubau
(1950) solution for SPPs on a cylinder of finite conductivity. The external
radial component of the electric field $E_{r}$ has the form:

$E_{r}(r,z,t)=-\frac{\pi r_{pl}E_{0}}{2r_{outer}}e^{i(\omega
t-hz)}H^{(1)}_{1}(r\sqrt{k^{2}-h^{2}}),$ (3)

for waves of frequency $\omega$ propagating in the $+z$ direction, with an
amplitude of $E_{0}$ at the surface. The complex $h$ wavenumber describes the
SW wavelength and attenuation length scale, $k$ is the free space wavenumber
$k=\omega/c$, $H^{(1)}_{1}$ is a Hankel function of the first kind with order
1, and $r_{outer}=1/|\sqrt{k^{2}-h^{2}}|$. For radii $r_{pl}<r<r_{outer}$ the
Hankel function with complex argument is approximately
$H^{(1)}_{1}\simeq-2r_{outer}/(\pi r)$, and for larger radii it falls off
exponentially (as is typical for SPPs).

Recent work Jaisson (2014); Mendonça (2019) has simplified the calculation of
$h$ through the use of the Lambert $W$ function. Given a typical plasma
conductivity $\sigma$ a complex $\xi$ can be defined:

$\xi=-\frac{ke^{-2\gamma_{euler}}}{2}\sqrt{\frac{\epsilon_{0}\omega}{2\sigma}}r_{pl}(1+i),$
(4)

(where $\gamma_{euler}$ is the Euler-Mascheroni constant), and $h$ is then
given by:

$h=\sqrt{k^{2}+\frac{4\xi e^{-2\gamma_{euler}}}{r_{pl}^{2}W_{1}(\xi)}},$ (5)

where the first branch of $W$ is used. At 1 Torr the plasma has a radius
$r_{pl}\simeq 0.5$ mm and conductivity $\sigma\simeq 2\times 10^{3}$ S/m, and
thus for a 10 GHz SW we find $h\simeq 210+3.2i$. The dispersive SW thus
travels at about $0.988c$, attenuates by a factor of $1/e$ every $30$ cm (if
not being actively excited), and has a $1/r$ transverse envelope out to
$r_{outer}\simeq 2.4$ cm, all of which compares well with simulations.

Figure 3: Log plot of electric field strength shortly after the detachment of
a $\lambda=3.9$ $\mu$m, 10 Torr surface wave from a 2 cm plasma column, as
simulated by the axisymmetric FDTD-Drude model.

We next consider the excitation of the surface wave by the plasma wake
currents. An isolated Hertzian dipole with uniform current $I$ along an
antenna of length $L$ will radiate electromagnetic waves at a wavelength
$\lambda_{Dipole}$ with power:
$P_{Dipole}=(\pi\eta_{0}/3)I^{2}L^{2}/\lambda^{2}_{Dipole}$. In our case we
have a radial current distribution of approximate length $r_{\Delta}$ that is
normal to a conducting cylinder of radius $r_{pl}$. When there is a separation
of length scales $r_{\Delta}\ll r_{pl}\ll\lambda_{SW}$ and an overall radial
current of $I_{r}$ then the power that is primarily radiated into SWs is:

$P_{SW}\simeq\frac{\pi\eta_{0}}{60}I_{r}^{2}\frac{r_{\Delta}^{2}}{r_{pl}^{5/3}\lambda^{1/3}_{SW}},$
(6)

(where we have simplified with
$\lambda_{SW}=2\pi/\operatorname{\mathbb{R}e}(h)$). We note that the
dimensionless term $r_{\Delta}^{2}/(r_{pl}^{5/3}\lambda^{1/3}_{SW})$ (compare
to $L^{2}/\lambda^{2}_{Dipole}$ for the Hertzian dipole), has different
asymptotic scalings in other length scale regimes. Equation (6) applies both
to a stationary radial antenna and the $\pm z$ wavetrains it emits, and to our
case, where the antenna translates at $c$ with its comoving SW pulse.

The total effective current $I_{r}(z)$ is found by integrating the current
density $J_{r}$ from 0 to $z$, and absent dissipation the amplitude $E_{0}$ of
the SW in (3) would grow linearly (for constant $J_{r}$). However there is
dissipation due to the finite conductivity of the plasma, and the attenuation
coefficient $\alpha=\operatorname{\mathbb{I}m}(h)$ is similar in magnitude to
$1/L_{pl}$. The SW thus grows linearly initially and then saturates after
several $1/\alpha$, and the effective total radial current after propagating a
distance $z$ is:

$I_{r}(z)=2\pi r_{pl}J_{r}\frac{1-e^{-\alpha z}}{\alpha}.$ (7)

We express the amplitude of the final SW in terms of its longitudinal current,
as $I_{z}=2\pi r_{pl}E_{0}/\eta_{0}$. Integrating (3) to find the SW peak
power and setting this equal to (6) we find:

$I_{z}(z)=\frac{2\pi^{2}r_{\Delta}J_{r}(1-e^{-\alpha
z})}{\alpha(30\ln(r_{outer}/r_{pl}))^{1/2}}\left(\frac{r_{pl}}{\lambda_{SW}}\right)^{1/6}.$
(8)

For the $\lambda=800$ nm laser this gives $I_{z}\simeq 0.6$ A at the end of
the plasma at 1 Torr, a current of $I_{z}\simeq 0.2$ A at 10 Torr, and
$I_{z}\simeq 0.05$ A at 100 Torr, which are close to the measured currents in
the corresponding simulations.

When the surface wave reaches the end of plasma the majority of the energy
detaches and propagates forward as the RF pulse, with the reflected wave
decaying as it travels back down the plasma. Fig. 3 shows an example of this,
where the SW has just detached from a 2 cm plasma column produced by a
$\lambda=$ 3.9 $\mu$m laser source at 10 Torr. In this way the system
resembles a surface wave end-fire antenna Andersen (1967); Sporer  _et al._
(2017). The second half of the primarily one cycle RF pulse forms as the SW
separates from the plasma, with some variation in the profile dependent on the
abruptness of the detachment. The recorded $I_{z}(z,t)$ currents from the SW
simulation are used in a ${\sim}1$ m2 laboratory scale axisymmetric FDTD
model, with resolution $\Delta r=\Delta z=0.5$ mm. Plots of the resulting
Fourier transformed RF field strengths at 10 Torr for both $\lambda=800$ nm
and $\lambda=3.9$ $\mu$m simulations can be seen in Fig. 4 (see section 5 of
sup  for more detail).

Figure 4: Plots of simulated electric field strength as a function of
frequency and angle for RF pulses from a 40 mJ, f${}_{N}=60$, 800 nm laser
source (left, Englesbe  _et al._ (2018)) and a 12 mJ, f${}_{N}=15$, 3.9 $\mu$m
source (right, Mitrofanov  _et al._ (2021a)), with both cases run at 10 Torr.
The numerical profiles and amplitudes closely match the experimental results.
Figure 5: Polar plots of the detector peak to peak Voltage $V_{pp}$ for the
$\lambda=800$ nm laser source as a function of angle and pressure as measured
in the lab (left), and as calculated by both the simulation framework (right,
triangles), and by the associated analytic approximations (right, $\times$
ticks) that culminate in (9).

To directly compare with lab results the simulated electric fields are
converted to peak to peak voltages $V_{pp}$ via a frequency dependent
calibration factor $C_{E/V}(f_{GHz})\simeq 50+10f_{GHz}$ (which captures the
horn antenna response Englesbe  _et al._ (2018)). We adapt an equation for
energy spectral density radiated per unit solid angle derived for THz
radiation Amico  _et al._ (2008):

$\frac{d^{2}U}{d\omega
d\Omega}=\frac{|I_{z}(\omega)|^{2}\sin^{2}\theta}{4\pi\epsilon_{0}c(1-\cos\theta)^{2}}\sin^{2}\left(\frac{L_{pl}\omega}{2c}(1-\cos\theta)\right),$
(9)

to construct voltage profiles for our analytic approximations. The average
amplitude of $I_{z}$ (8) over the length of the plasma is applied to a
${\sim}10$ GHz pulse, which is Fourier transformed for use with (9) (along
with the calibration factor), and then converted back to an effective $V_{pp}$
(code that implements all of these analytic approximations is available in sup
).

In general simulations with varying pressures, focal geometries, and laser
pulse energies and wavelengths all produce RF that nicely matches recent
experimental results Englesbe  _et al._ (2018); Janicek  _et al._ (2020);
Mitrofanov  _et al._ (2021a); Englesbe  _et al._ (2021). The magnitude of the
simulated RF is a bit lower than the lab data for the runs shown in Fig. 5,
but is within the errors for both the simulations and in the experimental
data. Subtle changes in Ti:Sapphire laser alignment generates $\pm 50\%$
amplitude variations at a given pressure (with the spatial and frequency
content being much more stable). It is suspected that the slight laser pulse
variations are amplified during filamentation, leading to fine scale changes
in the plasma radius, density, and electron velocities. Plausible variations
in the same parameters lead to comparable changes in the simulated RF
amplitude as well.

We have demonstrated that a plasma wake field that excites a surface wave is
responsible for the generation of filament RF. The rich internal dynamics of
the plasma merits further research: we expect that stochastic higher frequency
RF may inform on the fine scale structure of the plasma and the nonlinear
optics that generate it. The ionization model needs to be expanded to include
more physical effects (e.g. Tochitsky  _et al._ ; Corkum (1993); Amini  _et
al._ (2019)), and strong field ionization can be further explored with noble
gas filament RF. The plasma can also be driven into the strongly coupled
regime with harder focusing, and we intend to model resulting RF with the use
of Baalrud and Daligault (2013); Donkó (2009). Finally we note that it may be
possible to recreate the core mechanism – an antenna that copropagates with
the wave that it excites – in other settings and thereby drive strong coherent
radiation.

###### Acknowledgements.

The authors thank Edward Ruden, Serge Kalmykov and Charles Armstrong for
useful discussions, and the Air Force Office of Scientific Research (AFOSR)
for support via Laboratory Tasks No. FA9550-18RDCOR0009 and No.
FA9550-19RDCOR027. This work was supported in part by high-performance
computer time and resources from the DoD High Performance Computing
Modernization Program. Approved for public release; distribution is unlimited.
Public Affairs release approval No. AFRL-2020-0489.

## References

  * Strickland and Mourou (1985) Donna Strickland and Gerard Mourou, “Compression of amplified chirped optical pulses,” Optics communications 56, 219–221 (1985).
  * Faure  _et al._ (2004) Jérôme Faure, Yannick Glinec, A Pukhov, S Kiselev, S Gordienko, E Lefebvre, J-P Rousseau, F Burgy, and Victor Malka, “A laser–plasma accelerator producing monoenergetic electron beams,” Nature 431, 541–544 (2004).
  * Liu  _et al._ (1997) X Liu, D Du, and G Mourou, “Laser ablation and micromachining with ultrashort laser pulses,” IEEE journal of quantum electronics 33, 1706–1716 (1997).
  * Di Piazza  _et al._ (2012) A Di Piazza, C Müller, KZ Hatsagortsyan, and Ch H Keitel, “Extremely high-intensity laser interactions with fundamental quantum systems,” Reviews of Modern Physics 84, 1177 (2012).
  * Braun  _et al._ (1995) A Braun, G Korn, X Liu, D Du, J Squier, and G Mourou, “Self-channeling of high-peak-power femtosecond laser pulses in air,” Optics letters 20, 73–75 (1995).
  * Brodeur and Chin (1999) A Brodeur and SL Chin, “Ultrafast white-light continuum generation and self-focusing in transparent condensed media,” JOSA B 16, 637–650 (1999).
  * Kandidov  _et al._ (2003) VP Kandidov, OG Kosareva, IS Golubtsov, W Liu, A Becker, N Akozbek, CM Bowden, and SL Chin, “Self-transformation of a powerful femtosecond laser pulse into a white-light laser pulse in bulk optical media (or supercontinuum generation),” Applied Physics B 77, 149–165 (2003).
  * Amico  _et al._ (2008) CD Amico, Aurélien Houard, Selcuk Akturk, Yi Liu, J Le Bloas, Michel Franco, Bernard Prade, Arnaud Couairon, VT Tikhonchuk, and André Mysyrowicz, “Forward thz radiation emission by femtosecond filamentation in gases: theory and experiment,” New Journal of Physics 10, 013015 (2008).
  * Andreeva  _et al._ (2016) VA Andreeva, OG Kosareva, NA Panov, DE Shipilo, PM Solyankin, MN Esaulkov, P González de Alaiza Martínez, AP Shkurinov, VA Makarov, L Bergé, _et al._ , “Ultrabroad terahertz spectrum generation from an air-based filament plasma,” Physical review letters 116, 063902 (2016).
  * Forestier  _et al._ (2010) Benjamin Forestier, Aurélien Houard, Magali Durand, Yves-Bernard André, Bernard Prade, J-Y Dauvignac, Franck Perret, Ch Pichot, Michel Pellet, and André Mysyrowicz, “Radiofrequency conical emission from femtosecond filaments in air,” Applied Physics Letters 96, 141111 (2010).
  * Englesbe  _et al._ (2018) Alexander Englesbe, Jennifer Elle, Remington Reid, Adrian Lucero, Hugh Pohle, Matthew Domonkos, Serge Kalmykov, Karl Krushelnick, and Andreas Schmitt-Sody, “Gas pressure dependence of microwave pulses generated by laser-produced filament plasmas,” Optics Letters 43, 4953–4956 (2018).
  * Janicek  _et al._ (2020) A Janicek, E Thornton, T Garrett, A Englesbe, J Elle, and A Schmitt-Sody, “Length dependence on broadband microwave emission from laser-generated plasmas,” IEEE Transactions on Plasma Science 48, 1979–1983 (2020).
  * Mitrofanov  _et al._ (2021a) AV Mitrofanov, DA Sidorov-Biryukov, MM Nazarov, AA Voronin, MV Rozhko, AB Fedotov, and AM Zheltikov, “Coherently enhanced microwave pulses from midinfrared-driven laser plasmas,” Optics Letters 46, 1081–1084 (2021a).
  * Englesbe  _et al._ (2021) Alexander Englesbe, Jennifer Elle, Robert Schwartz, Travis Garrett, Daniel Woodbury, Dogeun Jang, Ki-Yong Kim, Howard Milchberg, Remington Reid, Adrian Lucero, _et al._ , “Ultrabroadband microwave radiation from near-and mid-infrared laser-produced plasmas in air,” Physical Review A 104, 013107 (2021).
  * Mitrofanov  _et al._ (2021b) AV Mitrofanov, AA Voronin, MV Rozhko, DA Sidorov-Biryukov, MM Nazarov, AB Fedotov, and AM Zheltikov, “Polarization and spatial mode structure of mid-infrared-driven terahertz-to-microwave radiation,” ACS Photonics 8, 1988–1996 (2021b).
  * Kampfrath  _et al._ (2013) Tobias Kampfrath, Koichiro Tanaka, and Keith A Nelson, “Resonant and nonresonant control over matter and light by intense terahertz transients,” Nature Photonics 7, 680–690 (2013).
  * Labutin  _et al._ (2016) Timur A Labutin, Vasily N Lednev, Alexey A Ilyin, and Andrey M Popov, “Femtosecond laser-induced breakdown spectroscopy,” Journal of Analytical Atomic Spectrometry 31, 90–118 (2016).
  * Park and Field (2016) G Barratt Park and Robert W Field, “Perspective: The first ten years of broadband chirped pulse fourier transform microwave spectroscopy,” The Journal of chemical physics 144, 200901 (2016).
  * Doll  _et al._ (2013) Andrin Doll, Stephan Pribitzer, René Tschaggelar, and Gunnar Jeschke, “Adiabatic and fast passage ultra-wideband inversion in pulsed epr,” Journal of magnetic resonance 230, 27–39 (2013).
  * Capmany and Novak (2007) José Capmany and Dalma Novak, “Microwave photonics combines two worlds,” Nature photonics 1, 319–330 (2007).
  * Couairon and Mysyrowicz (2007) Arnaud Couairon and André Mysyrowicz, “Femtosecond filamentation in transparent media,” Physics reports 441, 47–189 (2007).
  * Chin (2010) See Leang Chin,  _Femtosecond laser filamentation_ , Vol. 55 (Springer, 2010).
  * Rodriguez and Dakovski (2010) George Rodriguez and Georgi L Dakovski, “Scaling behavior of ultrafast two-color terahertz generation in plasma gas targets: energy and pressure dependence,” Optics express 18, 15130–15143 (2010).
  * Birdsall and Langdon (2004) Charles K Birdsall and A Bruce Langdon,  _Plasma physics via computer simulation_ (CRC press, 2004).
  * Peterkin and Luginsland (2002) Robert E Peterkin and John W Luginsland, “A virtual prototyping environment for directed-energy concepts,” Computing in Science & Engineering 4, 42–49 (2002).
  * Taflove and Hagness (2005) Allen Taflove and Susan C Hagness,  _Computational electrodynamics: the finite-difference time-domain method_ (Artech house, 2005).
  * Yee (1966) Kane Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic mdeia,” IEEE Transactions on antennas and propagation 14, 302–307 (1966).
  * Boris (1970) Jay P Boris, “Relativistic plasma simulation-optimization of a hybrid code,” in  _Proc. Fourth Conf. Num. Sim. Plasmas_ (1970) pp. 3–67.
  * Villasenor and Buneman (1992) John Villasenor and Oscar Buneman, “Rigorous charge conservation for local electromagnetic field solvers,” Computer Physics Communications 69, 306–316 (1992).
  * Berenger (1994) Jean-Pierre Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of computational physics 114, 185–200 (1994).
  * (31) See the Supplemental Material at http:// link.aps.org/supplemental/10.1103/PhysRevE.104.L063201 for details on the numerical methods used in this Letter, along with code that implements the analytic approximations, and simulation movies.
  * Zhou  _et al._ (2011) Bing Zhou, Aurélien Houard, Yi Liu, Bernard Prade, André Mysyrowicz, Arnaud Couairon, Patrick Mora, Christopher Smeenk, Ladan Arissian, and Paul Corkum, “Measurement and control of plasma oscillations in femtosecond filaments,” Physical review letters 106, 255002 (2011).
  * Pitarke  _et al._ (2006) JM Pitarke, VM Silkin, EV Chulkov, and PM Echenique, “Theory of surface plasmons and surface-plasmon polaritons,” Reports on progress in physics 70, 1 (2006).
  * Maier  _et al._ (2006) Stefan A Maier, Steve R Andrews, L Martin-Moreno, and FJ Garcia-Vidal, “Terahertz surface plasmon-polariton propagation and focusing on periodically corrugated metal wires,” Physical Review Letters 97, 176805 (2006).
  * Popruzhenko  _et al._ (2008) SV Popruzhenko, VD Mur, VS Popov, and D Bauer, “Strong field ionization rate for arbitrary laser frequencies,” Physical review letters 101, 193003 (2008).
  * Sharma  _et al._ (2018) Animesh Sharma, Mikhail N Slipchenko, Mikhail N Shneider, Xingxing Wang, Kazi Arafat Rahman, and Alexey Shashurin, “Measurements of electron numbers in femtosecond laser induced plasmas using rayleigh microwave scattering,” in  _2018 AIAA Aerospace Sciences Meeting_ (2018) p. 0177.
  * (37) Sergei Tochitsky, Eric Welch, Mikhail Polyanskiy, Igor Pogorelsky, Paris Panagiotopoulos, Miroslav Kolesik, Ewan M Wright, Stephan W Koch, Jerome V Moloney, Jeremy Pigeon, _et al._ , “Megafilament in air formed by self-guided terawatt long-wavelength infrared laser,” Nature Photonics 13.
  * Kolesik and Moloney (2004) Miroslav Kolesik and Jerome V Moloney, “Nonlinear optical pulse propagation simulation: From maxwell’s to unidirectional equations,” Physical Review E 70, 036604 (2004).
  * Bergé  _et al._ (2004) L Bergé, S Skupin, F Lederer, G Méjean, J Yu, J Kasparian, E Salmon, JP Wolf, M Rodriguez, L Wöste, _et al._ , “Multiple filamentation of terawatt laser pulses in air,” Physical Review Letters 92, 225002 (2004).
  * Ettoumi  _et al._ (2015a) Wahb Ettoumi, Jérôme Kasparian, and J-P Wolf, “Laser filamentation as a new phase transition universality class,” Physical Review Letters 114, 063903 (2015a).
  * Ettoumi  _et al._ (2015b) Wahb Ettoumi, Jérôme Kasparian, and J-P Wolf, “Spin-glass model governs laser multiple filamentation,” Physical review letters 115, 033902 (2015b).
  * Corkum (1993) Paul B Corkum, “Plasma perspective on strong field multiphoton ionization,” Physical review letters 71, 1994 (1993).
  * Pointon (2008) Timothy D Pointon, “Second-order, exact charge conservation for electromagnetic particle-in-cell simulation in complex geometry,” Computer Physics Communications 179, 535–544 (2008).
  * Itikawa (2006) Yukikazu Itikawa, “Cross sections for electron collisions with nitrogen molecules,” Journal of physical and chemical reference data 35, 31–53 (2006).
  * Thorne and Blandford (2017) Kip S Thorne and Roger D Blandford,  _Modern Classical Physics: Optics, Fluids, Plasmas, Elasticity, Relativity, and Statistical Physics_ (Princeton University Press, 2017).
  * Dawson (1959) John M Dawson, “Nonlinear electron oscillations in a cold plasma,” Physical review 113, 383 (1959).
  * Sakai  _et al._ (2018) Osamu Sakai, Akinori Iwai, Yoshiharu Omura, Satoshi Iio, and Teruki Naito, “Wave propagation in and around negative-dielectric-constant discharge plasma,” Physics of Plasmas 25, 031901 (2018).
  * Okoniewski  _et al._ (1997) M Okoniewski, M Mrozowski, and MA Stuchly, “Simple treatment of multi-term dispersion in fdtd,” IEEE Microwave and Guided Wave Letters 7, 121–123 (1997).
  * Kretschmann (1972) E Kretschmann, “The angular dependence and the polarisation of light emitted by surface plasmons on metals due to roughness,” Optics Communications 5, 331–336 (1972).
  * Smith and Purcell (1953) Stephen Judson Smith and EM Purcell, “Visible light from localized surface charges moving across a grating,” Physical Review 92, 1069 (1953).
  * Dahan  _et al._ (2020) Raphael Dahan, Saar Nehemia, Michael Shentcis, Ori Reinhardt, Yuval Adiv, Xihang Shi, Orr Be’er, Morgan H Lynch, Yaniv Kurman, Kangpeng Wang, _et al._ , “Resonant phase-matching between a light wave and a free-electron wavefunction,” Nature Physics 16, 1123–1131 (2020).
  * Sommerfeld (1899) Arnold Sommerfeld, “Ueber die fortpflanzung elektrodynamischer wellen längs eines drahtes,” Annalen der Physik 303, 233–290 (1899).
  * Goubau (1950) Georg Goubau, “Surface waves and their application to transmission lines,” Journal of Applied Physics 21, 1119–1128 (1950).
  * Jaisson (2014) Denis Jaisson, “Simple formula for the wave number of the goubau line,” Electromagnetics 34, 85–91 (2014).
  * Mendonça (2019) J Ricardo G Mendonça, “Electromagnetic surface wave propagation in a metallic wire and the lambert w function,” American Journal of Physics 87, 476–484 (2019).
  * Andersen (1967) JB Andersen, “Radiation from surface-wave antennas,” Electronics Letters 3, 251–252 (1967).
  * Sporer  _et al._ (2017) Michael Sporer, Robert Weigel, and Alexander Koelpin, “A 24 ghz dual-polarized and robust dielectric rod antenna,” IEEE Transactions on Antennas and Propagation 65, 6952–6959 (2017).
  * Amini  _et al._ (2019) Kasra Amini, Jens Biegert, Francesca Calegari, Alexis Chacón, Marcelo F Ciappina, Alexandre Dauphin, Dmitry K Efimov, Carla Figueira de Morisson Faria, Krzysztof Giergiel, Piotr Gniewek, _et al._ , “Symphony on strong field approximation,” Reports on Progress in Physics 82, 116001 (2019).
  * Baalrud and Daligault (2013) Scott D Baalrud and Jérôme Daligault, “Effective potential theory for transport coefficients across coupling regimes,” Physical review letters 110, 235001 (2013).
  * Donkó (2009) Zoltán Donkó, “Molecular dynamics simulations of strongly coupled plasmas,” Journal of Physics A: Mathematical and Theoretical 42, 214029 (2009).