# On the validity of the guiding-center approximation in a magnetic dipole
field

Alain J. Brizard1,a and Danielle G. Markowski1,2 1Department of Physics, Saint
Michael’s College, Colchester, VT 05439, USA
2Department of Physics and Astronomy, Clemson University, Clemson, SC
29634-0978, USA
aAuthor to whom correspondence should be addressed<EMAIL_ADDRESS>
###### Abstract

The charged-particle dynamics in an axisymmetric magnetic-dipole geometry is
used to assess the validity of Hamiltonian guiding-center theory, which
includes higher-order corrections associated with guiding-center polarization
induced by magnetic-field nonuniformity.

##  I Introduction

Higher-order guiding-center corrections Littlejohn_1983 ; Cary_Brizard_2009 ;
Burby_2013 ; Tronko_Brizard_2015  play an important role in the faithful
representation of charged-particle dynamics expressed in terms of a reduced
dynamical description which decouples the fast gyromotion about a local
magnetic field line and the slow bounce and drift motions along and across
magnetic-field lines. The validity of the guiding-center approximation was
recently explored in a straight magnetic field with constant gradients
Brizard_2017 . There, it was shown that the guiding-center approximation
remains valid even in the presence of strong gradients.

Here, we investigate how a truncated guiding-center expression for an exact
invariant for particle dynamics can accurately represent this particle
invariant. The higher-order guiding-center corrections have been shown to be
relevant in realistic magnetic-field geometries Belova_2003  for magnetically-
confined energetic charged particles, where the faithful guiding-center
representation of the toroidal angular momentum for these energetic particles
required higher-order corrections (up to first order in magnetic-field
nonuniformity).

The purpose of our present work is to explore how higher-order guiding-center
dynamics is faithful to charged single-particle dynamics in an axisymmetric
dipole magnetic field, in which the azimuthal canonical angular momentum is an
exact invariant. In particular, we compare the faithfulness of the truncated
expressions of the guiding-center azimuthal angular momentum derived in the
original work of Littlejohn Littlejohn_1983  and the recent work of Tronko and
Brizard Tronko_Brizard_2015 , which includes effects due to guiding-center
polarization Kaufman_1986 ; Brizard_2013

The remainder of the paper is organized as follows. In Sec. II, we present a
brief description of the axisymmetric magnetic-dipole geometry considered in
this work, while in Sec. III, we present the numerical solutions of the
equations of motion for a charged particle (with mass $m$ and charge $q$)
moving in an axisymmetric dipole magnetic field. These solutions satisfy the
conservation laws of energy and azimuthal angular canonical momentum exactly.
They also exhibit the expected three orbital time scales involving the fast
gyromotion, the intermediate bounce motion, and the slow azimuthal precession
motion. In Sec. IV, we present the guiding-center equations of motion
expressed as Euler-Lagrange equations derived from a Lagrangian that contains
higher-order guiding-center corrections. These equations conserve guiding-
center energy and guiding-center azimuthal angular canonical momentum exactly.
In Sec. V, the validity of the guiding-center approximation is explored in
magnetic-field dipole geometry. In particular, we compare the standard
guiding-center correction $-\,\mu\,({\bf R}+\frac{1}{2}\,\tau\,\widehat{{\sf
b}})$, first derived by Littlejohn Littlejohn_1983 , where ${\bf R}$ denotes
the gyrogauge vector field and $\tau\equiv\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}$ denotes the magnetic twist, with the guiding-center correction
$-\,\mu\,({\bf R}+\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}})$, which was derived Tronko_Brizard_2015  by requiring that the guiding-
center transformation yield the correct expression for the guiding-center
polarization Kaufman_1986 . It will be shown that, in magnetic-dipole geometry
(where the magnetic twist vanishes $\tau\equiv 0$), the guiding-center
correction $-\,\mu\,({\bf
R}+\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}})$ yields a
guiding-center description that is faithful to the exact particle dynamics.

##  II Magnetic-Dipole Geometry

We begin by considering charged-particle motion in magnetic-dipole geometry,
where the magnetic-dipole field is represented as

${\bf B}\;=\;\frac{B_{\rm e}\,r_{\rm
e}^{3}}{r^{3}}\left(2\,\sin\lambda\;\widehat{r}\;-\;\cos\lambda\,\widehat{\lambda}\right),$
(1)

where modified spherical coordinates $(r,\varphi,\lambda)$ are used (with the
latitude angle $-\pi/2\leq\lambda\leq\pi/2$ used instead of the polar angle
$\theta=\pi/2-\lambda$). The spatial Jacobian is ${\cal J}=(\nabla
r\,\mbox{\boldmath$\times$}\,\nabla\varphi\,\mbox{\boldmath$\cdot$}\,\nabla\lambda)^{-1}=r^{2}\,\cos\lambda$,
the unit vectors
$(\widehat{r},\widehat{\lambda}=\partial\widehat{r}/\partial\lambda)$ on the
poloidal plane are

$\left.\begin{array}[]{rcl}\widehat{r}&=&\cos\lambda\left(\cos\varphi\;\widehat{\sf
x}+\sin\varphi\;\widehat{\sf y}\right)\;+\;\sin\lambda\;\widehat{\sf z}\\\
\widehat{\lambda}&=&-\;\sin\lambda\left(\cos\varphi\;\widehat{\sf
x}+\sin\varphi\;\widehat{\sf y}\right)\;+\;\cos\lambda\;\widehat{\sf
z}\end{array}\right\\},$ (2)

and $B_{\rm e}=B_{\rm E}\,(r_{\rm E}/r_{\rm e})^{3}$ denotes the strength of
the equatorial magnetic field (with $B_{\rm E}\simeq 0.312$ G and $r_{\rm
E}=6378$ km) at $(r,\lambda)=(r_{\rm e},0)$. The magnitude of the magnetic
field is expressed as

$B(r,\lambda)\;=\;\frac{B_{\rm e}\,r_{\rm
e}^{3}}{r^{3}}\;\sqrt{1+3\,\sin^{2}\lambda}\;\equiv\;B_{\rm e}\,\overline{B},$
(3)

while the unit vector along the magnetic field is

$\widehat{{\sf
b}}\;=\;\frac{2\,\sin\lambda\;\widehat{r}\;-\;\cos\lambda\,\widehat{\lambda}}{\sqrt{1+3\,\sin^{2}\lambda}}.$
(4)

In the present work, dimensionless variables and fields are normalized to
their equatorial values and are represented with an overbar (e.g.,
$\overline{r}=r/r_{\rm e}$).

The dipole magnetic field (1) may also be written in terms of the magnetic
vector potential ${\bf A}=\psi\,\nabla\varphi$ as

${\bf B}\;=\;\nabla\,\mbox{\boldmath$\times$}\,{\bf
A}\;=\;\nabla\psi\,\mbox{\boldmath$\times$}\,\nabla\varphi,$ (5)

which guarantees that the dipole magnetic field is divergenceless:
$\nabla\,\mbox{\boldmath$\cdot$}\,{\bf B}=0$, where the dipole magnetic flux
is

$\psi(r,\lambda)\;=\;\psi_{\rm e}\,(r_{\rm
e}/r)\;\cos^{2}\lambda\;\equiv\;\psi_{\rm e}\,\overline{\psi},$ (6)

with $\psi_{\rm e}=B_{\rm e}\,r_{\rm e}^{2}$.

Since the dipole magnetic field satisfies the condition ${\bf
B}\,\mbox{\boldmath$\cdot$}\,\nabla\psi=0$, we conclude that the dipole
magnetic-field lines (1) lie on a constant-$\psi$ surface $\psi=\psi_{\rm e}$,
i.e., $r=r_{\rm e}\,\cos^{2}\lambda$. On this constant-$\psi$ surface, the
magnitude (3) of the dipole magnetic field is a function of latitude $\lambda$
alone: $B(\lambda)=B_{\rm
e}\,(1+3\,\sin^{2}\lambda)^{\frac{1}{2}}/\cos^{6}\lambda$, which becomes
infinite as we approach the polar regions (i.e.,
$\lambda\rightarrow\pm\,\pi/2$). We should note, however, that the region near
$\lambda=\pm\,\pi/2$ on a single magnetic-field line is actually located well
within Earth’s surface.

Next, we note that, in the absence of current sources, the curl of the
magnetic field (1) vanishes: $\nabla\,\mbox{\boldmath$\times$}\,{\bf B}=0$,
which implies that

$\widehat{{\sf b}}\,\mbox{\boldmath$\times$}\,\nabla\ln
B\;=\;\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}\;=\;\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,(\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}})\;=\;K(\lambda)\;\nabla\varphi,$ (7)

where the magnetic curvature $\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}=K\,\nabla\varphi\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}$ is
expressed in terms of the dimensionless latitude-dependent magnetic-curvature
function

$K(\lambda)\;=\;\frac{3\;(1-\sin^{4}\lambda)}{(1+3\,\sin^{2}\lambda)^{\frac{3}{2}}}.$
(8)

This magnetic-curvature function will play a key role in our discussion of the
faithfulness of the guiding-center approximation.

##  III Phase-space Lagrangian Dynamics

In order to derive the equations of motion, we begin with the noncanonical
phase-space Lagrangian

$\displaystyle L$ $\displaystyle=$ $\displaystyle\left(\frac{q}{c}\,{\bf
A}\;+\;{\bf p}\right)\,\mbox{\boldmath$\cdot$}\,\dot{\bf
x}\;-\;\left(\frac{|{\bf p}|^{2}}{2m}\right)\;\equiv\;{\bf
P}\,\mbox{\boldmath$\cdot$}\,\dot{\bf x}\;-\;H$ (9) $\displaystyle=$
$\displaystyle\frac{q}{c}\,\psi\,\dot{\varphi}\;+\;\frac{m}{2}\left(\dot{r}^{2}+r^{2}\,\dot{\lambda}^{2}+r^{2}\,\cos^{2}\lambda\;\dot{\varphi}^{2}\right),$

where the canonical momentum ${\bf P}=q{\bf A}/c+{\bf p}$ is the sum of the
kinetic momentum ${\bf p}=m\,{\bf v}$ and the time-independent magnetic vector
potential ${\bf A}=\psi\,\nabla\varphi$.

The Euler-Lagrange equations are expressed in terms of the spatial coordinates
$(r,\lambda,\varphi)$ in dimensionless form:

$\displaystyle\overline{r}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\overline{r}\,\lambda^{\prime
2}\;+\;\varphi^{\prime}\cos^{2}\lambda\left(\overline{r}\,\varphi^{\prime}-\epsilon^{-1}/\overline{r}^{2}\right),$
(10) $\displaystyle\lambda^{\prime\prime}$ $\displaystyle=$
$\displaystyle-\;\varphi^{\prime}\cos\lambda\sin\lambda\left(\varphi^{\prime}+\epsilon^{-1}/\overline{r}^{3}\right)\;-\;2\,\overline{r}^{\prime}\lambda^{\prime}/\overline{r},$
(11) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(\overline{\psi}\;+\frac{}{}\epsilon\,\varphi^{\prime}\;\overline{r}^{2}\,\cos^{2}\lambda\right)^{\prime},$
(12)

where all time derivatives are normalized to $\epsilon\,\Omega_{\rm
e}\ll\Omega_{\rm e}=qB_{\rm e}/mc$. The fast gyromotion time scale at $r_{\rm
e}=2\,r_{\rm E}$ is $2\pi/\Omega_{\rm e}\simeq 9\,\mu{\rm sec}$ for electrons
and $2\pi/\Omega_{\rm e}\simeq 17$ msec for protons. In addition, at $r_{\rm
e}=2\,r_{\rm E}$, the particle velocity is normalized to $\epsilon\,r_{\rm
e}\Omega_{\rm e}\simeq 16\,\epsilon\,c$ (for protons), which implies that
$\epsilon\ll 1/16$ in order to satisfy the non-relativistic assumption.

We note that, according to Eq. (12), the covariant azimuthal component of the
canonical momentum ${\bf P}$:

$\displaystyle P_{\varphi}$ $\displaystyle=$ $\displaystyle{\bf
P}\,\mbox{\boldmath$\cdot$}\,\partial{\bf
x}/\partial\varphi\;=\;(q/c)\,\psi\;+\;m\,\dot{\varphi}\,|\partial{\bf
x}/\partial\varphi|^{2}$ $\displaystyle=$ $\displaystyle(q/c)\psi_{\rm
e}\left(\overline{r}^{-3}+\epsilon\,\varphi^{\prime}\right)\overline{r}^{2}\,\cos^{2}\lambda\;\equiv\;(q/c)\psi_{\rm
e}\,\overline{P}_{\varphi}$

is an exact dynamical invariant, which follows from the azimuthal symmetry of
the dipole magnetic field.

Figure 1: Plots of (left) the normalized radial position $r(t)/r_{\rm
e}=\overline{r}(t)$, (center) the latitude $\lambda(t)$, and (right) the
normalized magnetic flux $\psi(t)/\psi_{\rm
e}=\cos^{2}\lambda(t)/\overline{r}(t)$ over a full bounce period for
$\epsilon=1/50$. The fast oscillations are due to the gyromotion of charged
particles about a single magnetic-field line while the slow bounce motion is
seen in the left and center frames. The normalized gyroradius is
$|\mbox{\boldmath$\rho$}_{0}|/r_{\rm e}=[2\mu B_{\rm e}/(mr_{\rm
e}^{2}\Omega_{\rm
e}^{2})]^{\frac{1}{2}}=\epsilon\,\overline{\mu}^{\frac{1}{2}}\ll 1$. Figure 2:
Plot of the azimuthal angle $\varphi(t)$ over a full bounce period, with the
fast oscillations due to gyromotion.

Figures 1-2 show the numerical solutions of the equations of motion (10)-(12)
over one bounce-motion (normalized) period $T_{\rm b}\simeq 2.5$, with initial
conditions $(\overline{r}_{0},\lambda_{0},\varphi_{0})=(1,0,0)$ and
$(\overline{r}^{\prime},\lambda_{0}^{\prime},\varphi_{0}^{\prime})=(0,1,1)$,
so that a particle orbit begins on the equatorial plane. These figures show
the fast gyromotion oscillations and the slow bounce-motion modulation, as a
particle oscillates between latitude values $-\lambda_{\rm
b}\leq\lambda\leq\lambda_{\rm b}$, as can be seen in the center frame of Fig.
1. We note that because of finite-gyroradius effects, a charged particle does
not move on a constant-$\psi$ surface (right frame of Fig. 1) and its
departure is described by the gyroradius vector

$\displaystyle\mbox{\boldmath$\rho$}_{0}$ $\displaystyle=$
$\displaystyle\frac{\widehat{{\sf b}}}{\Omega}\,\mbox{\boldmath$\times$}\,{\bf
v}\;=\;\frac{1}{B\Omega}\left(\dot{\psi}\;\nabla\varphi\;-\;\dot{\varphi}\;\nabla\psi\right)$
(14) $\displaystyle=$
$\displaystyle\frac{1}{\Omega}\left(\frac{\dot{\psi}}{|\nabla\psi|}\;\widehat{\varphi}\;-\;\frac{\dot{\varphi}}{|\nabla\varphi|}\;\widehat{\psi}\right)\;\equiv\;\epsilon\,r_{\rm
e}\,\overline{\mbox{\boldmath$\rho$}}_{0},$

where
$B=|\nabla\psi\,\mbox{\boldmath$\times$}\,\nabla\varphi|=|\nabla\psi|\,|\nabla\varphi|$
and
$(\widehat{\varphi},\widehat{\psi})=(\nabla\varphi/|\nabla\varphi|,\nabla\psi/|\nabla\psi|)$
are covariant unit vectors. Here, using Eq. (6), we find

$\displaystyle\nabla\psi$ $\displaystyle=$ $\displaystyle-\,\psi_{\rm
e}\,(r_{\rm e}/r^{2})\;\cos\lambda\;\left(\cos\lambda\;\nabla
r+2\,r\,\sin\lambda\frac{}{}\nabla\lambda\right),$ $\displaystyle\dot{\psi}$
$\displaystyle=$ $\displaystyle-\,\psi_{\rm e}\,(r_{\rm
e}/r^{2})\;\cos\lambda\;\left(\cos\lambda\;\dot{r}+2\,r\,\sin\lambda\frac{}{}\dot{\lambda}\right),$

so that

$\frac{\dot{\psi}}{|\nabla\psi|}\;=\;-\,\frac{\cos\lambda\;\dot{r}+2\,r\,\sin\lambda\frac{}{}\dot{\lambda}}{\sqrt{1+3\,\sin^{2}\lambda}}$
(15)

and
$\widehat{\psi}=-\,(\cos\lambda\,\widehat{r}+2\,\sin\lambda\,\widehat{\lambda})/\sqrt{1+3\,\sin^{2}\lambda}$.

Figure 3: Plot of the normalized parallel particle velocity
$\overline{v}_{\|}=v_{\|}/(\epsilon\,r_{\rm e}\Omega_{\rm e})$, with $v_{\|}$
given by Eq. (III), over a bounce period. The particle reaches a bounce
turning-point $\pm\lambda_{\rm b}$ as $v_{\|}$ vanishes and reverses sign.

The kinetic energy ${\cal E}=m|{\bf v}|^{2}/2$ is also an exact dynamical
invariant since $\dot{\cal E}=m\dot{\bf v}\,\mbox{\boldmath$\cdot$}\,{\bf
v}=(q/c)\,({\bf v}\,\mbox{\boldmath$\times$}\,{\bf
B})\,\mbox{\boldmath$\cdot$}\,{\bf v}=0$. Here, it is possible to separate the
parallel and perpendicular components of the particle velocity as follows.
First, we define the parallel particle velocity

$\displaystyle v_{\|}$ $\displaystyle=$ $\displaystyle{\bf
v}\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\;=\;\frac{2\,\dot{r}\,\sin\lambda\;-\;r\,\dot{\lambda}\,\cos\lambda}{\sqrt{1+3\,\sin^{2}\lambda}}$
$\displaystyle=$ $\displaystyle\epsilon\,r_{\rm e}\Omega_{\rm
e}\left(\frac{2\,\overline{r}^{\prime}\,\sin\lambda-\overline{r}\,\lambda^{\prime}\,\cos\lambda}{\sqrt{1+3\,\sin^{2}\lambda}}\right)\equiv\epsilon\,r_{\rm
e}\Omega_{\rm e}\,\overline{v}_{\|},$

while the perpendicular velocity can be expressed as

$\displaystyle{\bf v}_{\bot}$ $\displaystyle=$ $\displaystyle(\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,{\bf
v})\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}\;=\;\frac{\dot{\psi}}{|\nabla\psi|}\;\widehat{\psi}\;+\;\frac{\dot{\varphi}}{|\nabla\varphi|}\;\widehat{\varphi}$
(17) $\displaystyle\equiv$ $\displaystyle\epsilon\,r_{\rm e}\Omega_{\rm
e}\,\overline{\mbox{\boldmath$\rho$}}_{0}\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}.$

The parallel velocity is shown in Fig. 3 over one bounce period, where large-
amplitude gyro-oscillations are observed as the particle moves near the
equatorial plane. The 3D particle orbit is shown in Fig. 4, where we observe
the fast gyromotion (spiraling orbit about a single magnetic-field line),
intermediate bounce motion (projected on the $z$-axis) between turning points
where the parallel velocity (III) vanishes.

Figure 4: Normalized particle orbit (over several bounce periods) showing the
fast gyromotion, the intermediate bounce motion (predominantly along the
$z$-axis), and the slow azimuthal precession motion (projected onto the
$(x,y)$ equatorial plane).

##  IV Guiding-center Phase-space Lagrangian Dynamics

In order to remove the fast gyromotion time scale $(\ll 1\,{\rm sec})$ from
the equations of motion (10)-(12), and extract slower orbital bounce
$(>1\,{\rm min})$ and azimuthal precession $(>1\,{\rm hour})$ time scales, we
proceed with the guiding-center representation for charged-particle motion in
a nonuniform magnetic field (see Refs. Littlejohn_1983 ; Cary_Brizard_2009 ;
Tronko_Brizard_2015 ).

Using the guiding-center coordinates $({\bf X},P_{\|},\mu,\zeta)$, the phase-
space Lagrangian (9) is transformed into the guiding-center phase-space
Lagrangian Tronko_Brizard_2015

$\displaystyle L_{\rm gc}$ $\displaystyle=$
$\displaystyle\left(\frac{q}{c}\,{\bf A}\;+\;P_{\|}\;\widehat{{\sf
b}}\right)\,\mbox{\boldmath$\cdot$}\,\dot{{\bf
X}}\;+\;J\,\left(\dot{\zeta}\;-\;{\bf
R}^{*}\,\mbox{\boldmath$\cdot$}\,\dot{{\bf X}}\right)$ (18)
$\displaystyle-\;\left(\frac{P_{\|}^{2}}{2m}\;+\;\mu\,B\right),$

where ${\bf X}=(R,\Phi,\Lambda)$ denotes the guiding-center position, $P_{\|}$
denotes the parallel guiding-center momentum, and the guiding-center
gyroaction $J=\mu\,B_{\rm e}/\Omega_{\rm e}=\mu\,(mc/q)$ is canonically
conjugate to the guiding-center gyroangle $\zeta$. The guiding-center position
${\bf X}({\bf x},{\bf v})={\bf x}-\epsilon\,\mbox{\boldmath$\rho$}_{0}-\cdots$
is defined with the help of the (lowest-order) gyroradius vector (14), so that
the guiding-center spherical coordinates are

$\left.\begin{array}[]{rcl}R&=&r\;-\;\epsilon\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla
r\;+\;\cdots\\\
\Lambda&=&\lambda\;-\;\epsilon\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\lambda\;+\;\cdots\\\
\Phi&=&\varphi\;-\;\epsilon\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\varphi\;+\;\cdots\end{array}\right\\},$
(19)

where higher-order corrections depend on the magnetic-field nonuniformity. The
guiding-center transformation from the local momentum coordinates
$(p_{\|}=\widehat{{\sf b}}\,\mbox{\boldmath$\cdot$}\,{\bf p},J_{0}=m|{\bf
v}\,\mbox{\boldmath$\times$}\,\widehat{{\sf b}}|^{2}/2\Omega,\zeta_{0})$ to
the guiding-center momentum coordinates $(P_{\|},J\equiv\mu\,B_{\rm
e}/\Omega_{\rm e},\zeta)$, on the other hand, depends on magnetic-field
nonuniformity Cary_Brizard_2009 . In the present work, the first-order
correction $\mu_{1}$, which is derived in the Appendix, will be needed in
order to verify the validity of the guiding-center approximation.

The vector functions ${\bf A}=\Psi\,\nabla\Phi$ and ${\bf B}=B\,\widehat{{\sf
b}}$ in Eq. (18) are evaluated at the guiding-center position ${\bf X}$. We
note that, as a result of the guiding-center transformation, the magnetic
moment $\mu$ is an invariant of the guiding-center dynamics, while it is not
an invariant of the particle dynamics. In Eq. (18), the guiding-center
correction is expressed in terms of the guiding-center vector field

${\bf R}^{*}\;=\;{\bf
R}\;+\;\left\\{\begin{array}[]{lr}\frac{1}{2}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}\;\equiv\;\frac{1}{2}\,K\,\nabla\Phi&({\sf A})\\\ &\\\
\frac{1}{2}\,(\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}})\;\widehat{{\sf b}}\;\equiv\;0&({\sf B})\end{array}\right.$ (20)

defined in terms of the gyrogauge vector field Littlejohn_1983  ${\bf
R}=\nabla\widehat{1}\,\mbox{\boldmath$\cdot$}\,\widehat{2}$, which is
expressed in terms of the orthogonal unit vector fields
$(\widehat{1},\widehat{2},\widehat{{\sf
b}}=\widehat{1}\,\mbox{\boldmath$\times$}\,\widehat{2})$, and a higher-order
correction that either involves (A) magnetic curvature Tronko_Brizard_2015  or
(B) a correction that vanishes in magnetic-dipole geometry Littlejohn_1983 .
In magnetic-dipole geometry, where we can use $\widehat{1}=\widehat{\Phi}$ and
$\widehat{2}=\widehat{{\sf b}}\,\mbox{\boldmath$\times$}\,\widehat{\Phi}$, we
find ${\bf
R}=\nabla\Phi\,(\partial\widehat{\Phi}/\partial\Phi)\,\mbox{\boldmath$\cdot$}\,(\widehat{{\sf
b}}\,\mbox{\boldmath$\times$}\,\widehat{\Phi})=b_{z}(\Lambda)\,\nabla\Phi$,
where $\partial\widehat{\Phi}/\partial\Phi=\widehat{\sf
z}\,\mbox{\boldmath$\times$}\,\widehat{\Phi}$ and $b_{z}\equiv\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\widehat{\sf z}$, while the dimensionless
magnetic curvature $K$ is defined in Eq. (8). Hence, the vector field ${\bf
R}^{*}=b_{z}^{*}\;\nabla\Phi$ is expressed in terms of the curvature-modified
function

$\displaystyle b_{z}^{*}$ $\displaystyle\equiv$ $\displaystyle
b_{z}+\frac{1}{2}\,K\;=\;\frac{9\,\sin^{4}\Lambda-1}{(1+3\,\sin^{2}\Lambda)^{\frac{3}{2}}}\;+\;\frac{3}{2}\frac{(1-\sin^{4}\Lambda)}{(1+3\,\sin^{2}\Lambda)^{\frac{3}{2}}}$
(21) $\displaystyle=$
$\displaystyle\frac{(15\,\sin^{4}\Lambda+1)}{2\,(1+3\,\sin^{2}\Lambda)^{\frac{3}{2}}}.$

Figure 5 shows that, while $b_{z}(\Lambda)$ remains negative until
$\Lambda>\arcsin(1/\sqrt{3})\simeq\pi/5\,(\simeq 35^{o})$, the curvature-
modified function $b_{z}^{*}(\Lambda)$ remains positive throughout the range
$0\leq\Lambda\leq\pi/2$.

###  IV.1 Guiding-center Euler-Lagrange equations

The guiding-center Euler-Lagrange equations associated with variations of the
guiding-center Lagrangian (18) in $({\bf X},P_{\|},\mu,\zeta)$ yield,
respectively,

$\displaystyle\dot{P}_{\|}\,\widehat{{\sf b}}-\frac{q}{c}\,\dot{\bf
X}\,\mbox{\boldmath$\times$}\,{\bf B}^{*}-\dot{J}\,{\bf R}^{*}$
$\displaystyle=$ $\displaystyle-\,\nabla H_{\rm gc}=-\,\mu\nabla B,$ (22)
$\displaystyle\widehat{{\sf b}}\,\mbox{\boldmath$\cdot$}\,\dot{\bf X}$
$\displaystyle=$ $\displaystyle\partial H_{\rm gc}/\partial P_{\|}=P_{\|}/m,$
(23) $\displaystyle\dot{\zeta}\;-\;\dot{\bf X}\,\mbox{\boldmath$\cdot$}\,{\bf
R}^{*}$ $\displaystyle=$ $\displaystyle\partial H_{\rm gc}/\partial J=\Omega,$
(24) $\displaystyle\dot{J}$ $\displaystyle=$ $\displaystyle-\;\partial H_{\rm
gc}/\partial\zeta=0,$ (25)

where the curvature-modified magnetic field

$\displaystyle{\bf B}^{*}$ $\displaystyle=$
$\displaystyle\nabla\,\mbox{\boldmath$\times$}\,\left[{\bf
A}\;+\;(c/q)\left(P_{\|}\;\widehat{{\sf b}}\;-\frac{}{}\epsilon\,J\,{\bf
R}^{*}\right)\right]$ (26) $\displaystyle=$
$\displaystyle\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi\;+\;(c/q)\;P_{\|}K(\Lambda)\,\nabla\Phi$

is defined in terms of the gyroaction-corrected magnetic flux

$\Psi^{*}\;=\;\Psi(R,\Lambda)\;-\;(c/q)\,J\,b_{z}^{*}(\Lambda).$ (27)

Although Eq. (25) implies that the guiding-center gyroaction $J$ is a guiding-
center invariant (since $H_{\rm gc}$ is independent of the gyroangle $\zeta$
up to second order in magnetic-field nonuniformity), it is not an exact
particle invariant.

Making use of Eq. (25), Eqs. (22)-(23) yield the reduced guiding-center
equations of motion

$\displaystyle\dot{\bf X}$ $\displaystyle=$
$\displaystyle\frac{P_{\|}}{m}\;\frac{{\bf
B}^{*}}{B_{\|}^{*}}\;+\;\frac{c\widehat{{\sf
b}}}{qB_{\|}^{*}}\,\mbox{\boldmath$\times$}\,\mu\,\nabla B$ $\displaystyle=$
$\displaystyle\frac{c}{qB_{\|}^{*}}\left(\mu\,B+\frac{P_{\|}^{2}}{m}\right)K\,\nabla\Phi\;+\;\frac{P_{\|}}{mB_{\|}^{*}}\,\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi,$
$\displaystyle\dot{P}_{\|}$ $\displaystyle=$ $\displaystyle-\;\frac{{\bf
B}^{*}}{B_{\|}^{*}}\,\mbox{\boldmath$\cdot$}\,\mu\nabla
B\;=\;-\;\frac{1}{B_{\|}^{*}}\;\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi\,\mbox{\boldmath$\cdot$}\,\mu\nabla
B,$ (29)

where we used Eq. (7) and

$\displaystyle B_{\|}^{*}$ $\displaystyle=$ $\displaystyle\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf B}^{*}\;=\;\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi$
$\displaystyle=$ $\displaystyle\frac{-\,1}{{\cal
J}\sqrt{1+3\,\sin^{2}\Lambda}}\left(2\,\sin\Lambda\;\frac{\partial\Psi^{*}}{\partial\Lambda}+R\,\cos\Lambda\;\frac{\partial\Psi^{*}}{\partial
R}\right),$

where ${\cal J}=(\nabla
R\,\mbox{\boldmath$\times$}\,\nabla\Phi\,\mbox{\boldmath$\cdot$}\,\nabla\Lambda)^{-1}=R^{2}\,\cos\Lambda$
is the spherical Jacobian. In Eq. (29), we used the fact that the magnitude of
the dipole magnetic field $B$ is independent of the guiding-center azimuthal
angle $\Phi$ (i.e., $\nabla\Phi\,\mbox{\boldmath$\cdot$}\,\nabla B=0$). We
note that $P_{\|}=m\dot{\bf X}\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf b}}$ is
indeed the guiding-center parallel momentum and the fast gyromotion, now
described by Eq. (24), is decoupled from the slow guiding-center dynamics
(IV.1)-(29).

Figure 5: Plots of $b_{z}^{*}$ (solid) and $b_{z}$ (dashed) as functions of
the guiding-center latitude angle $0\leq\Lambda\leq\pi/2$. The impact of the
dimensionless magnetic curvature $b_{z}^{*}-b_{z}=\frac{1}{2}\,K$ is clearly
seen with $b_{z}^{*}>0$ in the range $0\leq\Lambda\leq\pi/2$.

In addition, we note that, while the guiding-center energy ${\cal
E}=P_{\|}^{2}/2m+\mu\,B$ is an exact guiding-center invariant:

$\dot{\cal E}\;=\;P_{\|}\,\dot{P}_{\|}/m\;+\;\mu\,\dot{\bf
X}\,\mbox{\boldmath$\cdot$}\,\nabla B\;=\;0,$

the guiding-center azimuthal canonical angular momentum

$P_{{\rm gc}\Phi}\;=\;(q/c)\,\Psi^{*}\;=\;(q/c)\,\Psi\;-\;J\;b_{z}^{*}$ (31)

is also an exact guiding-center invariant, since

$\dot{\Psi}^{*}\;=\;\dot{\bf
X}\,\mbox{\boldmath$\cdot$}\,\nabla\Psi^{*}\;=\;\dot{R}\;\frac{\partial\Psi^{*}}{\partial
R}\;+\;\dot{\Lambda}\;\frac{\partial\Psi^{*}}{\partial\Lambda}\;=\;0,$ (32)

where, from Eq. (IV.1), we used

$\displaystyle\dot{R}$ $\displaystyle=$ $\displaystyle\dot{\bf
X}\,\mbox{\boldmath$\cdot$}\,\nabla
R\;=\;\frac{P_{\|}}{mB_{\|}^{*}}\,\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi\,\mbox{\boldmath$\cdot$}\,\nabla
R$ (33) $\displaystyle=$ $\displaystyle-\;\frac{P_{\|}}{m{\cal
J}B_{\|}^{*}}\,\frac{\partial\Psi^{*}}{\partial\Lambda},$
$\displaystyle\dot{\Lambda}$ $\displaystyle=$ $\displaystyle\dot{\bf
X}\,\mbox{\boldmath$\cdot$}\,\nabla\Lambda\;=\;\frac{P_{\|}}{mB_{\|}^{*}}\,\nabla\Psi^{*}\,\mbox{\boldmath$\times$}\,\nabla\Phi\,\mbox{\boldmath$\cdot$}\,\nabla\Lambda$
(34) $\displaystyle=$ $\displaystyle\frac{P_{\|}}{m{\cal
J}B_{\|}^{*}}\,\frac{\partial\Psi^{*}}{\partial R}.$

Lastly, the guiding-center azimuthal angular velocity

$\dot{\Phi}\;=\;\dot{\bf
X}\,\mbox{\boldmath$\cdot$}\,\nabla\Phi\;=\;\frac{c}{qB_{\|}^{*}}\left(\mu\,B+\frac{P_{\|}^{2}}{m}\right)K\,|\nabla\Phi|^{2}$
(35)

is entirely driven by magnetic curvature, while the guiding-center parallel
force equation (29) is driven by the magnetic-field gradient.

###  IV.2 Normalized guiding-center equations

In order to extract explicit orbital time scales from the guiding-center
equations of motion (29) and (33)-(35), we now introduce the following
normalized definitions. First, we introduce $B=B_{\rm e}\,\overline{B}$,
$R=r_{\rm e}\,\rho$, and $\Psi=B_{\rm e}r_{\rm e}^{2}\,\overline{\Psi}$, so
that

$\Psi^{*}\;=\;B_{\rm e}r_{\rm e}^{2}\left(\overline{\Psi}\;-\;\frac{\nu_{\rm
e}^{2}}{2\;\Omega_{\rm e}^{2}}\;b_{z}^{*}(\Lambda)\right)\;\equiv\;B_{\rm
e}r_{\rm e}^{2}\,\overline{\Psi}^{*},$ (36)

where $\overline{B}=\rho^{-3}(1+3\sin^{2}\Lambda)^{\frac{1}{2}}$ and
$\overline{\Psi}=\rho^{-1}\cos^{2}\Lambda$, and we introduced the energy-
dependent bounce frequency

$\nu_{\rm e}\;\equiv\;\sqrt{\frac{2\,\mu B_{\rm e}}{mr_{\rm
e}^{2}}}\;\equiv\;\epsilon\,\Omega_{\rm e}\,\overline{\mu}^{\frac{1}{2}},$
(37)

and the equatorial radius $r_{\rm e}>r_{\rm E}$ is above the surface of the
Earth (note that the dipole-field approximation for Earth’s magnetic field,
however, is valid for $r_{\rm e}<3\,r_{\rm E}$).

Next, the normalized guiding-center motion in the poloidal
$(\rho,\Lambda)$-plane

$\left(\dot{\rho},\;\dot{\Lambda}\right)\;=\;\frac{\nu_{\rm
e}\,\overline{P}_{\|}}{\overline{\cal
J}\overline{B}_{\|}^{*}}\,\left(-\;\frac{\partial\overline{\Psi}^{*}}{\partial\Lambda},\frac{\partial\overline{\Psi}^{*}}{\partial\rho}\right),$
(38)

clearly takes places on a time scale associated with the frequency $\nu_{\rm
e}$, where we introduced the normalized guiding-center parallel momentum
$\overline{P}_{\|}=P_{\|}/(mr_{\rm e}\nu_{\rm e})$, which satisfies the
normalized guiding-center parallel-force equation

$\dot{\overline{P}}_{\|}\;=\;-\;\frac{\nu_{\rm e}}{2\overline{\cal
J}\overline{B}_{\|}^{*}}\left(\frac{\partial\overline{\Psi}^{*}}{\partial\rho}\,\frac{\partial\overline{B}}{\partial\Lambda}\;-\;\frac{\partial\overline{\Psi}^{*}}{\partial\Lambda}\,\frac{\partial\overline{B}}{\partial\rho}\right),$
(39)

so that the total guiding-center kinetic energy becomes

${\cal E}\;=\;\mu\,B_{\rm
e}\left(\overline{B}\;+\;\overline{P}_{\|}^{2}\right)\;=\;\mu\,B_{\rm
t}\;\equiv\;\mu B_{\rm e}\;\overline{B}_{\rm t},$ (40)

where $B_{\rm t}$ denotes the magnetic field at which the parallel guiding-
center momentum vanishes. Here, we easily see that
$|\overline{P}_{\|}|=\sqrt{\overline{B}_{t}-\overline{B}}$ (which has a
maximum at the equator $|\overline{P}_{\|{\rm e}}|=\sqrt{\overline{B}_{\rm
t}-1}$), where $1\leq\overline{B}\leq\overline{B}_{\rm t}={\cal E}/(\mu B_{\rm
e})$, and the poloidal guiding-center motion (38) also comes to a stop at a
turning point. For this reason, the equations (38) represent the periodic
bounce motion between turning points located symmetrically in the northern and
southern hemispheres. Lastly, the normalized guiding-center azimuthal motion
is expressed as

$\dot{\Phi}\;=\;\frac{\nu_{\rm e}^{2}}{2\,\Omega_{\rm
e}\,\overline{B}_{\|}^{*}}\left(\overline{B}\;+\;2\,\overline{P}_{\|}^{2}\right)\,K\,|\overline{\nabla}\Phi|^{2},$
(41)

which does not stop at a turning point since
$\overline{B}+2\,\overline{P}_{\|}^{2}=2\,\overline{B}_{\rm
t}-\overline{B}\geq\overline{B}_{\rm t}$ does not vanish.

We conclude this Section by noting that Eqs. (38) and (41) imply that time
derivatives can be normalized to the bounce time scale $\tau\equiv\nu_{\rm
e}t$, so that we obtain the normalized guiding-center equations of motion

$\displaystyle\left(\rho^{\prime}(\tau),\;\Lambda^{\prime}(\tau)\right)$
$\displaystyle=$ $\displaystyle\frac{\overline{P}_{\|}}{\overline{\cal
J}\overline{B}_{\|}^{*}}\,\left(-\;\frac{\partial\overline{\Psi}^{*}}{\partial\Lambda},\frac{\partial\overline{\Psi}^{*}}{\partial\rho}\right),$
(42) $\displaystyle\Phi^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle\frac{\nu_{\rm e}}{2\,\Omega_{\rm
e}\,\overline{B}_{\|}^{*}}\left(\overline{B}\;+\;2\,\overline{P}_{\|}^{2}\right)\,K\,|\overline{\nabla}\Phi|^{2},$
(43) $\displaystyle\overline{P}_{\|}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle\frac{-1}{2\overline{\cal
J}\overline{B}_{\|}^{*}}\left(\frac{\partial\overline{\Psi}^{*}}{\partial\rho}\,\frac{\partial\overline{B}}{\partial\Lambda}-\frac{\partial\overline{\Psi}^{*}}{\partial\Lambda}\,\frac{\partial\overline{B}}{\partial\rho}\right).$
(44)

For electrons (with $\mu B_{\rm e}=100$ keV) or protons (with $\mu B_{\rm
e}=100$ MeV) at $r_{\rm e}=2\,r_{\rm E}$, for example, we find $\nu_{\rm
e}=10.4$ Hz (electrons) or 7.67 Hz (protons) and $\nu_{\rm e}/\Omega_{\rm
e}\simeq 1.5\times 10^{-5}$ (electrons) or $2\times 10^{-2}$ (protons). Hence,
the azimuthal precession time scale is approximately $(\Omega_{\rm e}/\nu_{\rm
e})$ longer than the bounce time scale (i.e., on the scale of a few hours
compared to a few seconds). The bounce-center analysis Cary_Brizard_2009  of
particle motion in an axisymmetric dipole magnetic field was carried out in
Ref. Duthoit_2010 .

Figure 6: Guiding-center orbit shown over an entire precession period
($\tau_{\rm p}\simeq 115$) based on the numerical integration of Eqs.
(42)-(44) with normalized parameters $\overline{B}_{\rm t}={\cal E}/\mu B_{\rm
e}=2$ (i.e., $\overline{P}_{\|{\rm e}}=1$) and $\nu_{\rm e}/\Omega_{\rm
e}=1/50$ used in Eq. (43). Figure 7: Guiding-center parallel momentum
$-\,\overline{P}_{\|{\rm
e}}\leq\overline{P}_{\|}(\tau)\leq\overline{P}_{\|{\rm e}}$ as a function of
normalized time $\tau=\nu_{\rm e}t$ over one bounce period $\tau_{\rm b}\simeq
2.5$ for $\overline{P}_{\|{\rm e}}=1$.

Figure 6 shows a guiding-center orbit obtained by integrating Eqs. (42)-(44)
with the normalized parameters $\overline{B}_{\rm t}={\cal E}/\mu B_{\rm e}=2$
(i.e., $-1\leq\overline{P}_{\|}\leq 1$ with $\overline{P}_{\|{\rm e}}=1$) and
$\nu_{\rm e}/\Omega_{\rm e}=2\times 10^{-2}$ used in Eq. (43). With these
parameters, the normalized bounce period is $\tau_{\rm b}\simeq 2.5$ (see Fig.
7) and the normalized precession period is $\tau_{\rm p}\simeq 115\gg\tau_{\rm
b}$ (see Fig. 8). These guiding-center figures can be compared with the
particle figures 4, 3, and 2, respectively, where the fast gyromotion time
scale has been removed by the guiding-center transformation, while the bounce
and precession time scales remain intact.

Figure 8: Guiding-center azimuthal angle $\Phi(\tau)$ as a function of
normalized time $\tau=\nu_{\rm e}t$ (solid curve) over one bounce period
$\tau_{\rm b}\simeq 2.5$. The dashed curve is a bounce-averaged straight line
($2\pi\tau/\tau_{\rm p}$) whose slope is adjusted to a normalized precession
period $\tau_{\rm p}\simeq 115\gg\tau_{\rm b}$.

##  V Validity of the Guiding-center Approximation

Figure 9: Plot of the normalized magnetic moments $\overline{\mu}_{0}$ (light
curve) and $\overline{\mu}=\overline{\mu}_{0}+\epsilon\,\overline{\mu}_{1}$
(dark curve) over a bounce period for $\epsilon=1/50$.

We now proceed with a numerical comparison of the exact particle dynamics
based on Eqs. (10)-(12) with the guiding-center dynamics based on Eqs.
(42)-(44). In particular, we verify that the particle azimuthal canonical
momentum (III) and the guiding-center azimuthal canonical momentum (31) are
indeed exact invariants of their respective dynamics.

One of the important tests of the validity of the guiding-center approximation
is concerned with how well the guiding-center magnetic moment $\mu$ is
conserved by the particle dynamics. In Fig. 9, we see that the normalized
zeroth-order magnetic moment $\overline{\mu}_{0}=2\,\mu_{0}B_{\rm
e}/(\epsilon^{2}m\,r_{\rm e}^{2}\Omega_{\rm e}^{2})$ is not well conserved
during a bounce period. When first-order corrections are added to the magnetic
moment $\mu=\mu_{0}+\epsilon\,\mu_{1}$, however, we see that the normalized
magnetic moment
$\overline{\mu}=\overline{\mu}_{0}+\epsilon\,\overline{\mu}_{1}$ (dark curve)
shows very small oscillations, even when the particle crosses the equatorial
plane. The first-order correction $\mu_{1}$, which is calculated explicitly in
App. A for the magnetic dipole geometry, is seen to play a crucial role.

Figure 10: Plots of the normalized guiding-center pull-back ${\sf T}_{\rm
gc}\overline{P_{{\rm gc}\Phi}}={\sf T}_{\rm
gc}\overline{\Psi}-\epsilon^{2}\,(\overline{\mu}/2)\;{\sf T}_{\rm
gc}b_{z}^{*}$ (dark curve) and $\overline{P_{\varphi}}$ (black horizontal
line) with
$\overline{\Psi}(\overline{r},\lambda)=\cos^{2}\lambda/\overline{r}$ (light
curve).

Another important test explores the accuracy and faithfulness of the guiding-
center azimuthal canonical momentum (31) by transforming it back into particle
phase space and use the particle dynamics to assess its invariance. This task
allows us to determine the role of the higher-order curvature-corrected
function $b_{z}^{*}=b_{z}+\frac{1}{2}\,K$ Tronko_Brizard_2015 , defined in Eq.
(21), in contrast to the standard guiding-center term $b_{z}$ Littlejohn_1983
. We investigate how faithful the guiding-center azimuthal angular momentum
$P_{{\rm gc}\Phi}$, defined by Eq. (31), is to the exact azimuthal angular
momentum $P_{\varphi}$, defined by Eq. (III). To perform a valid comparison,
we need to transform the guiding-center azimuthal angular momentum $P_{{\rm
gc}\Phi}$ back into particle phase space, with the help of the guiding-center
pull-back operator ${\sf T}_{\rm gc}$. Here, using Eq. (19), the pull-back of
an arbitrary guiding-center function $F_{\rm gc}({\bf X})$ is defined on
particle phase space as

${\sf T}_{\rm gc}F_{\rm gc}({\bf x},{\bf p})\;=\;F_{\rm gc}\left({\bf
x}-\epsilon\mbox{\boldmath$\rho$}_{0}-\epsilon^{2}\mbox{\boldmath$\rho$}_{1}-\cdots\right),$
(45)

which replaces the guiding-center position ${\bf X}={\bf
x}-\epsilon\,\mbox{\boldmath$\rho$}_{0}-\epsilon^{2}\mbox{\boldmath$\rho$}_{1}-\cdots$
with the corrected particle position ${\bf x}$ and finite-Larmor-radius (FLR)
corrections involving the lowest-order gyroradius $\mbox{\boldmath$\rho$}_{0}$
and a correction $\mbox{\boldmath$\rho$}_{1}$ associated with magnetic-field
nonuniformity. Hence, the pull-back of the azimuthal guiding-center canonical
angular momentum (31) is expressed as

${\sf T}_{\rm gc}P_{{\rm gc}\Phi}\;=\;\frac{q}{c}{\sf T}_{\rm
gc}\Psi\;-\;\left\\{\begin{array}[]{lr}\epsilon^{2}\,J\;{\sf T}_{\rm
gc}b_{z}^{*}&({\sf A})\\\ &\\\ \epsilon^{2}\,J\;{\sf T}_{\rm gc}b_{z}&({\sf
B})\end{array}\right.$ (46)

where $b_{z}^{*}=b_{z}+\frac{1}{2}\,K$ adds the guiding-center polarization
(curvature) term $K$ to the standard guiding-center correction Littlejohn_1983
; Belova_2003 . First, we see in Fig. 10 that the normalized guiding-center
pull-back (solid)

${\sf T}_{\rm gc}\overline{P_{{\rm gc}\Phi}}\;=\;{\sf T}_{\rm
gc}\overline{\Psi}\;-\;\epsilon^{2}\,(\overline{\mu}/2)\;{\sf T}_{\rm
gc}b_{z}^{*}$

is very close to the exact (normalized) particle azimuthal canonical angular
momentum $\overline{P_{\varphi}}$ (solid horizontal line). The lowest-order
component $\overline{\Psi}(\overline{r},\lambda)=\cos^{2}\lambda/\overline{r}$
is also shown as a light curve, which highlights the effects of the guiding-
center pull-back operation.

Figure 11: Plot of the relative difference $({\sf T}_{\rm gc}P_{{\rm
gc}\Phi})/P_{\varphi}-1$ between the pull-back of the azimuthal guiding-center
canonical angular momentum (46) and the particle azimuthal canonical angular
momentum (III). The dark curve uses the guiding-center polarization correction
(A) in Eq. (46) while the light curve only shows the standard correction (B).

Second, the relative difference $({\sf T}_{\rm gc}P_{{\rm
gc}\Phi})/P_{\varphi}-1$ between the pull-back of the azimuthal guiding-center
canonical angular momentum (46) and the particle azimuthal canonical angular
momentum (III) is shown in Fig. 11, where excellent agreement is observed when
the guiding-center polarization curvature term $\frac{1}{2}\,K$ (dark curve)
is added to the standard correction $b_{z}$ (light curve). Hence, the guiding-
center polarization correction introduced in Ref. Tronko_Brizard_2015  yields
an azimuthal guiding-center canonical angular momentum that is faithful to the
particle azimuthal canonical angular momentum.

##  VI Summary

In previous work Brizard_2017 , we showed that the guiding-center
approximation was valid in a straight magnetic field with constant
perpendicular magnetic gradient, even in the presence of strong gradients. In
addition, based on the existence of an exact analytical solution for the
particle orbits, this work also confirmed that the guiding-center polarization
corresponded exactly with an orbit-averaged particle displacement.

In the present work, we extended our investigation of the validity of the
guiding-center approximation in describing charged single-particle motion in a
nonuniform magnetic field. Here, we considered particle motion in azimuthally
symmetric magnetic dipole geometry, in which the azimuthal angular canonical
momentum is conserved. We successfully validated the guiding-center
approximation in describing particle motion in an azimuthally symmetric
magnetic dipole field provided higher-order guiding-center corrections are
taken into account, which had already been noted for the case of an
axisymmetric tokamak magnetic field Belova_2003 . In particular, the guiding-
center polarization correction in the guiding-center azimuthal angular
canonical momentum, not taken into account in the standard guiding-center
approximation Littlejohn_1983 , proved crucial

Data Availability Statement

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

##  Appendix A Higher-order Guiding-center Approximation

In this Appendix, we will calculate explicit guiding-center corrections that
arise at first-order in the magnetic-field nonuniformity associated with
magnetic-field gradients $(\nabla\ln B)$ and magnetic-field curvature
$(\mbox{\boldmath$\kappa$}=\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf b}})$.

First, we explore the adiabatic invariance of the magnetic moment
$\mu_{0}=m|{\bf v}_{\bot}|^{2}/2B$. Indeed, we find that
$\dot{\mu}_{0}=-\,\mu_{0}\,{\bf v}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B+m\dot{\bf v}_{\bot}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}$, which can
ultimately be expressed as

$\displaystyle\dot{\mu}_{0}$ $\displaystyle=$ $\displaystyle-\,\mu_{0}\,{\bf
v}\,\mbox{\boldmath$\cdot$}\,\nabla\ln B+\frac{m}{B}\left(\dot{\bf
v}-\dot{v}_{\|}\,\widehat{{\sf b}}-v_{\|}\,{\bf
v}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\right)\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}$ (47) $\displaystyle=$
$\displaystyle-\,\mu_{0}\,{\bf v}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B-\frac{mv_{\|}}{B}\;{\bf v}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}$ $\displaystyle=$
$\displaystyle-\,\frac{v_{\|}}{B}\left(\mu_{0}\,\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla B\;+\frac{}{}m\,{\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}\right)$
$\displaystyle-\;\frac{{\bf
v}_{\bot}}{B}\,\mbox{\boldmath$\cdot$}\,\left(\mu_{0}\,\nabla
B\;+\;m\,v_{\|}^{2}\,\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf b}}\right).$

Although $\dot{\mu}_{0}\neq 0$ for general magnetic fields, its average over
the fast gyromotion time scale yields

$\langle\dot{\mu}_{0}\rangle\;=\;-\,\frac{\mu_{0}\,v_{\|}}{B}\left(\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla
B\;+\frac{}{}B\;\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\right)\;=\;-\,\frac{\mu_{0}\,v_{\|}}{B}\;(\nabla\,\mbox{\boldmath$\cdot$}\,{\bf
B}),$

where we used $\langle{\bf v}_{\bot}\rangle=0$ and $\langle{\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf
v}_{\bot}\rangle=(\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf b}})\;|{\bf
v}_{\bot}|^{2}/2$. Since magnetic fields are divergenceless, we immediately
find that $\langle\dot{\mu}_{0}\rangle=0$, i.e., $\mu_{0}$ is an invariant
over time scales that are long compared to the fast gyromotion time scale.

###  A.1 First-order correction to $\mu_{0}$

A guiding-center adiabatic invariant for the magnetic moment can be
constructed $\mu=\mu_{0}+\epsilon\,\mu_{1}+\cdots$ as an expansion in powers
of $\epsilon$, where $\mu_{1}$ is a correction that involves the nonuniformity
of the magnetic field. For a general magnetic field, the first-order
correction to the magnetic moment is Tronko_Brizard_2015

$\displaystyle\mu_{1}$ $\displaystyle=$
$\displaystyle\left(\mu_{0}\;\nabla_{\bot}\ln
B\;+\;\frac{p_{\|}^{2}\,\mbox{\boldmath$\kappa$}}{m\,B}\right)\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}$
(48)
$\displaystyle-\;\mu_{0}\;\varrho_{\|}\left(\tau\;+\frac{}{}\alpha_{1}\right),$

where
$\alpha_{1}=\frac{1}{2}\,\tau-\widehat{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\widehat{\rho}$. Since $\tau=\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\,\mbox{\boldmath$\times$}\,\widehat{{\sf
b}}=0$ and $\nabla_{\bot}\ln B=\mbox{\boldmath$\kappa$}$ in magnetic dipole
geometry, we find

$-\,\mu_{0}\varrho_{\|}\alpha_{1}\;=\;\frac{p_{\|}}{2B}\;\left(\frac{d\widehat{{\sf
b}}}{dt}\;-\;v_{\|}\;\mbox{\boldmath$\kappa$}\right)\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0},$

so that Eq. (48) becomes

$\mu_{1}\;=\;\left(\mu_{0}\;+\;\frac{p_{\|}^{2}}{2m\,B}\right)\mbox{\boldmath$\kappa$}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}\;+\;\frac{p_{\|}}{2B}\;\frac{d\widehat{{\sf
b}}}{dt}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0},$ (49)

where

$\mbox{\boldmath$\kappa$}\,\mbox{\boldmath$\cdot$}\,\mbox{\boldmath$\rho$}_{0}\;=\;-\,\frac{K(\lambda)}{\Omega}\;\dot{\varphi}\;=\;-\;\epsilon\,\frac{\varphi^{\prime}}{\overline{B}}\;K.$
(50)

We note here that the first-order correction (49) is specifically designed so
that the expression for $\dot{\mu}_{0}+\epsilon\dot{\mu}_{1}$ is given as

$\displaystyle\dot{\mu}_{0}+\epsilon\dot{\mu}_{1}$ $\displaystyle=$
$\displaystyle-\,\frac{v_{\|}}{B}\left(\mu_{0}\,\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla B\;+\frac{}{}m\,{\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}\right)\;-\;\frac{{\bf
v}_{\bot}}{B}\,\mbox{\boldmath$\cdot$}\,\left(\mu_{0}\,\nabla_{\bot}B\;+\;m\,v_{\|}^{2}\,\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf b}}\right)$ (51)
$\displaystyle+\;\left(\mu_{0}\;\nabla_{\bot}\ln
B\;+\;\frac{p_{\|}^{2}\,\mbox{\boldmath$\kappa$}}{m\,B}\right)\,\mbox{\boldmath$\cdot$}\,{\bf
v}_{\bot}\;+\;\frac{p_{\|}}{2B}\;\left({\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}-\frac{\partial{\bf
v}_{\bot}}{\partial\zeta}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\frac{\partial{\bf
v}_{\bot}}{\partial\zeta}\right)$ $\displaystyle=$
$\displaystyle-\,\frac{v_{\|}}{B}\;\mu_{0}\,\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla B\;-\;\frac{p_{\|}}{2B}\;\left({\bf
v}_{\bot}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,{\bf v}_{\bot}+\frac{\partial{\bf
v}_{\bot}}{\partial\zeta}\,\mbox{\boldmath$\cdot$}\,\nabla\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\frac{\partial{\bf
v}_{\bot}}{\partial\zeta}\right)$ $\displaystyle=$
$\displaystyle-\,\frac{\mu_{0}\,v_{\|}}{B}\left(\widehat{{\sf
b}}\,\mbox{\boldmath$\cdot$}\,\nabla
B\;+\frac{}{}B\;\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\right)\;=\;-\,\frac{\mu_{0}\,v_{\|}}{B}\;(\nabla\,\mbox{\boldmath$\cdot$}\,{\bf
B})\;=\;0,$

where we used $\epsilon\dot{\mu}_{1}=\Omega\,\partial\mu_{1}/\partial\zeta$ to
lowest order and
$\dot{\mbox{\boldmath$\rho$}}_{0}=\Omega\,\partial\mbox{\boldmath$\rho$}_{0}/\partial\zeta={\bf
v}_{\bot}$. Hence, the higher-order approximation
$\mu=\mu_{0}+\epsilon\,\mu_{1}$ of the guiding-center magnetic moment is
conserved up to first-order effects in magnetic-field nonunifornity, without
the need of gyroangle-averaging.

###  A.2 First-order correction to $\mbox{\boldmath$\rho$}_{0}$

When FLR corrections (19) need to be calculated up to first order in magnetic-
field nonuniformity, we must evaluate the pull-back of the guiding-center
position up to second order in $\epsilon$:

$\displaystyle{\sf T}_{\rm gc}{\bf X}$ $\displaystyle=$ $\displaystyle{\bf
x}+\epsilon\,G_{1}^{\bf x}+\epsilon^{2}\left(G_{2}^{\bf
x}\;+\;\frac{1}{2}\,{\sf G}_{1}\cdot{\sf d}G_{1}^{\bf x}\right)+\cdots$ (52)
$\displaystyle\equiv$ $\displaystyle{\bf x}-\mbox{\boldmath$\rho$}\;=\;{\bf
x}-\epsilon\,\mbox{\boldmath$\rho$}_{0}\;-\;\epsilon^{2}\,\mbox{\boldmath$\rho$}_{1}+\cdots$

where, using the components $(G_{2}^{\bf x},G_{1}^{\mu},G_{1}^{\zeta})$
calculated in Ref. Tronko_Brizard_2015 , we obtain the first-order gyroradius

$\displaystyle\mbox{\boldmath$\rho$}_{1}$ $\displaystyle=$
$\displaystyle-\,G_{2}^{\bf
x}\;-\;\frac{1}{2}\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\mbox{\boldmath$\rho$}_{0}\;+\;\frac{1}{2}\,\left(G_{1}^{\mu}\,\frac{\partial\mbox{\boldmath$\rho$}_{0}}{\partial\mu}+G_{1}^{\zeta}\,\frac{\partial\mbox{\boldmath$\rho$}_{0}}{\partial\zeta}\right)$
(53) $\displaystyle\equiv$ $\displaystyle\rho_{1\|}\,\widehat{{\sf
b}}\;+\;\frac{1}{2}\,\left(\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B\right)\mbox{\boldmath$\rho$}_{0}\;-\;\frac{\mu
B}{2m\Omega^{2}}\;\mbox{\boldmath$\kappa$}.$

Here, $\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\ln
B=-\,\dot{\varphi}\,K/\Omega$ and the parallel component is defined as

$\rho_{1\|}=\frac{5}{4}\left[\frac{\mu
B}{2m\Omega^{2}}\left(\nabla\,\mbox{\boldmath$\cdot$}\,\widehat{{\sf
b}}\right)-\frac{v_{\|}}{\Omega}\,K\,\mbox{\boldmath$\rho$}_{0}\,\mbox{\boldmath$\cdot$}\,\nabla\varphi\right]-\frac{3}{4\Omega^{2}}\;\frac{d\widehat{{\sf
b}}}{dt}\,\mbox{\boldmath$\cdot$}\,{\bf v}.$

This first-order correction can then be used in the guiding-center pull-back
(45).

We now calculate the guiding-center gyroradius $\mbox{\boldmath$\rho$}_{\rm
gc}\equiv{\sf T}_{\rm gc}^{-1}{\bf x}-{\bf X}={\sf T}_{\rm
gc}^{-1}\mbox{\boldmath$\rho$}$, which is defined in terms of the guiding-
center push-forward of the particle gyroradius into guiding-center phase
space, where $\mbox{\boldmath$\rho$}_{{\rm gc}0}=\mbox{\boldmath$\rho$}_{0}$
and $\mbox{\boldmath$\rho$}_{{\rm gc}1}=\mbox{\boldmath$\rho$}_{1}-{\sf
G}_{1}\cdot{\sf d}\mbox{\boldmath$\rho$}_{0}$. In Ref. Tronko_Brizard_2015 ,
it was shown that the guiding-center polarization can be constructed
explicitly from the guiding-center transformation, which is defined as the
gyroangle-averaged expression Brizard_2013

$\displaystyle\mbox{\boldmath$\pi$}_{\rm gc}$ $\displaystyle\equiv$
$\displaystyle q\,\langle\mbox{\boldmath$\rho$}_{{\rm
gc}1}\rangle\;-\;\nabla\,\mbox{\boldmath$\cdot$}\,\left\langle\frac{q}{2}\,\mbox{\boldmath$\rho$}_{0}\mbox{\boldmath$\rho$}_{0}\right\rangle$
$\displaystyle=$
$\displaystyle-\;\frac{q}{m\Omega^{2}}\left(\mu\,B\;\nabla_{\bot}\ln
B\;+\;\frac{p_{\|}^{2}}{m}\,\mbox{\boldmath$\kappa$}\right)\;=\;\frac{q\widehat{{\sf
b}}}{\Omega}\,\mbox{\boldmath$\times$}\,\frac{d{\bf X}}{dt},$

where the standard expression derived by Kaufman Kaufman_1986  is recovered.

## References

  * (1) R. G. Littlejohn, J. Plasma Phys. 29, 111 (1983).
  * (2) J. R. Cary and A. J. Brizard, Rev. Mod. Phys. 81, 693 (2009).
  * (3) J.W. Burby, J. Squire, and H. Qin, Phys. Plasmas 20, 072105 (2013).
  * (4) N. Tronko and A. J. Brizard, Phys. Plasmas 22, 112507 (2015).
  * (5) A. J. Brizard, Phys. Plasmas 24, 042115 (2017).
  * (6) E.V. Belova, N.N. Gorelenkov, and C.Z. Cheng, Phys. Plasmas 10, 3240 (2003).
  * (7) A.N. Kaufman, Phys. Fluids 29, 1736 (1986).
  * (8) A.J. Brizard, Phys. Plasmas 20, 092309 (2013).
  * (9) F.-X. Duthoit, A. J. Brizard, Y. Peysson, and J. Decker, Phys. Plasmas 17, 102903 (2010).


††institutetext: Institut für Theoretische Physik, Universität Tübingen,
Kepler Center for Astro and Particle Physics,
Auf der Morgenstelle 14, D-72076 Tübingen, Germany

# The massless non-adjacent double off-shell scalar box integral — branch cut
structure and all-order epsilon expansion

Juliane Haug 1 and Fabian Wunder111Corresponding author juliane-clara-
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We generalize the result of our recent paper on the massless single off-shell
scalar box integral to the case of two non-adjacent end points off the light
cone. An analytic result in $d=4-2\varepsilon$ dimensions is established in
terms of four Gauss hypergeometric ${{}_{2}}\mathrm{F}_{1}$ functions
respectively their single-valued counterparts. This allows for an explicit
splitting of real and imaginary parts, as well as an all-order
$\varepsilon$-expansion in terms of single-valued polylogarithms.

###### Keywords:

Feynman integrals, perturbative QCD, dimensional regularization, epsilon
expansion

††arxiv: 2302.01956

##  1 Introduction

The recent success of finding an all-order $\varepsilon$-expansion for the
single off-shell scalar box integral in terms of single-valued polylogarithms
(Haug2022, ) raised the question whether a similar result can be obtained in
the case of more than one end point off the light cone. Here, $\varepsilon$ is
the dimensional regularization parameter, with $d=4-2\varepsilon$ the number
of space-time dimensions. In this paper, we show that a direct generalization
of the method proposed in ref. (Haug2022, ) to two particles off the light
cone is indeed possible if they are on non-adjacent corners of the box. This
results in a hypergeometric representation which allows for an explicit
splitting of real and imaginary parts, as well as an all-order
$\varepsilon$-expansion. So far, the expansion was only available up to order
$\varepsilon^{0}$ in the literature (BernDixonKosower1, ; DuplancicNizic, ;
Ellis, ; Tarasov:2019mqy, ).

For adjacent particles off the light cone a similar straightforward
generalization does not seem feasible. Also in this case, results up to order
$\varepsilon^{0}$ have been known in the literature for some time
(BernDixonKosower1, ; DuplancicNizic, ; Ellis, ). To calculate higher orders
in the $\varepsilon$-expansion for the adjacent case, as well as for more
general kinematics, additional techniques will need to be employed. The
methods of differential equations (Gehrmann:1999as, ) or Mellin-Barnes
integrals (Smirnov:1999gc, ), which are highly successful in multi-loop
calculations, could be considered. A different approach towards the massless
box integral via negative dimensions is outlined in ref. (Anastasiou:1999cx,
). Yet another possibility are recurrence relations with respect to $d$, which
were employed in ref. (Fleischer:2003rm, ) to obtain a hypergeometric
representation of the box integral for general kinematics and masses.
Remarkably, a one-fold integral representation is presented there in eq. (96).
Using a wider set of functional equations, the box integral with massless
propagators and arbitrary external masses in terms of Appell functions
$\mathrm{F}_{1}$ and Gauss hypergeometric functions ${{}_{2}}\mathrm{F}_{1}$
was determined in ref. (Tarasov:2019mqy, ) (see eqs. (4.29) and (4.61) of that
reference). Comparing eqs. (4.29) and (4.32) of ref. (Tarasov:2019mqy, )
furthermore suggests that the adjacent double off-shell box integral is of the
same complexity as the general case. The method of functional equations has
recently been generalized to arbitrary internal masses (Tarasov:2022clb, ).
These results might be suitable starting points for generalizing the program
of the present paper towards an all-order $\varepsilon$-expansion with
explicit real and imaginary parts for kinematics beyond the scope of the
present work.

Since the calculation for two non-adjacent particles off the light cone is
very much analogous to the one presented in ref. (Haug2022, ) for the single
off-shell case, we kept the presentation more concise and refer the reader to
our previous paper for further details. This work is organized as follows. In
section 2, we show that we can keep two non-adjacent particles off the light
cone in the calculation of the box integral while closely following the steps
of our previous work. We obtain a representation in terms of four Gauss
hypergeometric functions ${{}_{2}}\mathrm{F}_{1}$, one more compared to the
single off-shell case. This procedure introduces spurious branch cuts. Since
the occurring ${{}_{2}}\mathrm{F}_{1}$ functions are exactly of the kind
present in the single-massive case, we can use their expansion known from ref.
(Haug2022, ) for the all-order $\varepsilon$-expansion. In section 3, we
establish the all-order $\varepsilon$-expansion of the non-adjacent double
off-shell box integral, with real and imaginary parts made explicit. The
result is free of spurious branch cuts due to the use of single-valued
polylogarithms.

##  2 Calculating the non-adjacent double off-shell scalar box integral

Figure 1: The scalar box diagram with two off-shell external particles
diagonally opposite from each other carrying momenta $p_{2}$ and $p_{4}$. It
is $p_{1}^{2}=p_{3}^{2}=0$ and $p_{2}^{2}=p_{4}^{2}\neq 0$. All external
momenta are taken to be incoming. Drawn with TikZ-Feynman (Ellis2017, ).

The general scalar box integral with massless internal lines is given by

$\displaystyle\mathrm{D}_{0}\,\equiv\,\frac{\mu^{4-d}}{\mathrm{i}\pi^{d/2}}\int\mathrm{d}^{d}l\,\frac{1}{\left[l^{2}+\mathrm{i}0\right]\left[(l+p_{2})^{2}+\mathrm{i}0\right]\left[(l+p_{2}+p_{3})^{2}+\mathrm{i}0\right]\left[(l-p_{1})^{2}+\mathrm{i}0\right]}\,,$
(1)

where the external momenta are labelled as indicated in the Feynman diagram
depicted in figure 1. We work in dimensional regularization with
$d=4-2\varepsilon$. While the starting point is the same as eq. (2.1) of our
recent paper (Haug2022, ), here we consider the possibility of more than one
external particle off the light cone. We have again explicitly kept the causal
$+\mathrm{i}0$ in the propagators, which is necessary to determine the
physical side of the branch cuts that appear in the final result.

After Feynman parametrization and evaluating the loop integral, we arrive at
the integral representation (see also eq. (7) of ref. (DuplancicNizic, ))

$\displaystyle D_{0}\,$
$\displaystyle=\,\mu^{2\varepsilon}\,\Gamma(2+\varepsilon)\int_{0}^{1}\frac{\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\mathrm{d}x_{3}\,\mathrm{d}x_{4}\,\delta(1-x_{1}-x_{2}-x_{3}-x_{4})}{\left[-x_{1}x_{2}p_{2}^{2}-x_{1}x_{3}s_{2}-x_{1}x_{4}p_{1}^{2}-x_{2}x_{3}p_{3}^{2}-x_{2}x_{4}s_{1}-x_{3}x_{4}p_{4}^{2}-\mathrm{i}0\right]^{2+\varepsilon}}\,,$
(2)

with the Mandelstam variables $s_{1}=(p_{1}+p_{2})^{2}$ and
$s_{2}=(p_{2}+p_{3})^{2}$. For convergence we require $\varepsilon<0$. In the
end, the result can be analytically continued to larger $\varepsilon$.

Using the substitution from our previous paper (see also refs.
(Smirnov_Feynman_Integrals, ; FabiValeryAlexey, )),

$\displaystyle x_{1}\,$ $\displaystyle=\,\eta_{1}\xi_{1}\,,\quad
x_{2}\,=\,\eta_{2}(1-\xi_{2})\,,\quad x_{3}\,=\,\eta_{2}\xi_{2}\,,\quad
x_{4}\,=\,\eta_{1}(1-\xi_{1})\,,$ (3)

the term in the denominator becomes

$\displaystyle-\eta_{1}\eta_{2}(1-\xi_{1})(1-\xi_{2})\,s_{1}\,-\,\eta_{1}\eta_{2}\xi_{1}\xi_{2}\,s_{2}\,-\,\eta_{1}\eta_{2}\xi_{1}(1-\xi_{2})\,p_{2}^{2}\,-\,\eta_{1}\eta_{2}\xi_{2}(1-\xi_{1})\,p_{4}^{2}$
$\displaystyle-\,\eta_{1}^{2}\xi_{1}(1-\xi_{1})\,p_{1}^{2}\,-\,\eta_{2}^{2}\xi_{2}(1-\xi_{2})\,p_{3}^{2}\,-\mathrm{i}0\,.$
(4)

Here, we see that the $\eta$-integrals factorize from the $\xi$-integrals if
$p_{1}^{2}=p_{3}^{2}=0$. Hence, we can generalize the calculational method put
forward in the single off-shell case to the case of two non-adjacent particles
off the light cone.

We note that adjacent particles off the light cone require a different
treatment. This is because setting for example $p_{1}^{2}=p_{2}^{2}=0$ in eq.
(2) leads to a less symmetric Feynman parameter integral, where $x_{1}$
appears only once, $x_{2}$ and $x_{4}$ twice, and $x_{3}$ thrice. Therefore a
substitution in the style of (4) does not factorize the integral. Presumably,
establishing an all order $\varepsilon$-expansion in the case of adjacent
particles off the light cone is considerably more difficult and requires
additional techniques such as functional equations (Tarasov:2019mqy, ),
differential equations (Gehrmann:1999as, ) or Mellin-Barnes integrals
(Smirnov:1999gc, ). In the literature, it was also noted that the adjacent
integral is “hard” in comparison to the “easy” non-adjacent integral
(BernDixonKosower1, ; DuplancicNizic, ). For this reason we only discuss the
non-adjacent case in this work.

By evaluating the $\eta$-integrals in terms of Gamma functions and
substituting $\xi_{1}\rightarrow 1-\xi_{1}$, we arrive at

$\displaystyle
D_{0}(s_{1},s_{2},p_{1}^{2}=0,p_{2}^{2},p_{3}^{2}=0,p_{4}^{2})\,=\,\mu^{2\varepsilon}\frac{\Gamma(2+\varepsilon)\Gamma^{2}(-\varepsilon)}{\Gamma(-2\varepsilon)}$
$\displaystyle\times\int_{0}^{1}\mathrm{d}\xi_{1}\int_{0}^{1}\mathrm{d}\xi_{2}\left[-p_{2}^{2}-(s_{1}-p_{2}^{2})\xi_{1}-(s_{2}-p_{2}^{2})\xi_{2}-(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})\xi_{1}\xi_{2}-\mathrm{i}0\right]^{-2-\varepsilon},$
(5)

which in the $p_{2}^{2}\rightarrow 0$ limit reduces to eq. (2.11) of our
previous work (Haug2022, ). Note that
$s_{3}\equiv(p_{1}+p_{3})^{2}=p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}$. Analogous to
the procedure there, we factor out
$\left(\frac{(s_{1}-p_{2}^{2})(s_{2}-p_{2}^{2})}{p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}}-\mathrm{i}0\right)^{\\!-2-\varepsilon}$.
This leads to

$\displaystyle
D_{0}(s_{1},s_{2},0,p_{2}^{2},0,p_{4}^{2})=\frac{\Gamma(2+\varepsilon)\Gamma^{2}(-\varepsilon)}{\Gamma(-2\varepsilon)}\left(\frac{\mu^{2}(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})}{(s_{1}-p_{2}^{2})(s_{2}-p_{2}^{2})}+\mathrm{i}0\right)^{\\!\varepsilon}\frac{I_{123}(x_{1},x_{2},x_{3})}{(s_{1}-p_{2}^{2})(s_{2}-p_{2}^{2})}\,,$
(6)

where we introduced the dimensionless variables222To keep the notation close
to our previous work (Haug2022, ), we use $x_{1,2,3}$ for the dimensionless
variables. They are not to be confused with the Feynman parameters, which only
appear in eqs. (2) and (3).

$\displaystyle
x_{1}\equiv-\frac{p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}}{s_{1}-p_{2}^{2}}\,,\quad
x_{2}\equiv-\frac{p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}}{s_{2}-p_{2}^{2}}\,,\quad
x_{3}\equiv-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})p_{2}^{2}}{(s_{1}-p_{2}^{2})(s_{2}-p_{2}^{2})}\,,$
(7)

and

$\displaystyle
I_{123}(x_{1},x_{2},x_{3})=x_{1}x_{2}\int_{0}^{1}\mathrm{d}\xi_{1}\int_{0}^{1}\mathrm{d}\xi_{2}\left[x_{3}+x_{2}\xi_{1}+x_{1}\xi_{3}-x_{1}x_{2}\xi_{1}\xi_{2}+\mathrm{i}0\,\mathrm{sg}\right]^{-2-\varepsilon}\,,$
(8)

with
$\mathrm{sg}\equiv\mathrm{sgn}\\!\left(\frac{(s_{1}-p_{2}^{2})(s_{2}-p_{2}^{2})}{p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}}\right)$.
The integral $I_{123}(x_{1},x_{2},x_{3})$ can be evaluated along the same
lines as $I_{12}(x_{1},x_{2})$ in the single off-shell case. As a
generalization of eq. (2.22) of ref. (Haug2022, ) we find

$\displaystyle
I_{123}(x_{1},x_{2},x_{3})=\frac{(1+x_{3}-\mathrm{i}0\,\mathrm{sg})^{-1-\varepsilon}}{1+\varepsilon}$
$\displaystyle\left[I\left(\frac{x_{1}+x_{3}}{1+x_{3}}\right)+I\left(\frac{x_{2}+x_{3}}{1+x_{3}}\right)\right.$
$\displaystyle\left.-I\left(\frac{x_{1}+x_{2}+x_{3}-x_{1}x_{2}}{1+x_{3}}\right)-I\left(\frac{x_{3}}{1+x_{3}}\right)\right],$
(9)

where

$\displaystyle
I(\chi)\equiv\int_{0}^{\chi}\frac{\mathrm{d}\zeta}{1-\zeta}\left(\left[\zeta-\mathrm{i}0\,\tilde{\mathrm{sg}}\right]^{-1-\varepsilon}-1\right),$
(10)

with $\tilde{\mathrm{sg}}=\mathrm{sgn}\\!\left(1+x_{3}\right)\,\mathrm{sg}$.
We can evaluate $I(\chi)$ in terms of a hypergeometric function and an
additional $\ln(1-\chi-\mathrm{i}\tilde{0})$ as in eq. (2.26) of ref.
(Haug2022, ),

$\displaystyle I(\chi)\,$
$\displaystyle=\,-\frac{1}{\varepsilon}\left[\chi-\mathrm{i}0\,\tilde{\mathrm{sg}}\right]^{-\varepsilon}{{}_{2}}\mathrm{F}_{1}\\!\left(1,-\varepsilon,1-\varepsilon;\chi+\mathrm{i}\tilde{0}\right)+\,\ln(1-\chi-\mathrm{i}\tilde{0})\,.$
(11)

This splitting introduces a new regulator $\mathrm{i}\tilde{0}$ (see ref.
(Haug2022, ) for details), independently for each $I(\chi)$. The real parts of
the logarithms cancel between the first two and second two $I(\chi)$. To
achieve cancellation of the imaginary parts as well, we have to choose the
$\mathrm{i}\tilde{0}_{i}$ of the four integrals such that

$\displaystyle 0\,\overset{!}{=}\,$
$\displaystyle\mathrm{i}\pi\left[-\mathrm{sgn}\\!\left(\tilde{0}_{1}\right)\Theta\left(\frac{x_{1}-1}{1+x_{3}}\right)-\mathrm{sgn}\\!\left(\tilde{0}_{2}\right)\Theta\left(\frac{x_{2}-1}{1+x_{3}}\right)\right.$
$\displaystyle\hphantom{\mathrm{i}\pi}\left.\;\,+\mathrm{sgn}\\!\left(\tilde{0}_{3}\right)\Theta\left(-\frac{(1-x_{1})(1-x_{2})}{1+x_{3}}\right)+\mathrm{sgn}\\!\left(\tilde{0}_{4}\right)\Theta\left(-\frac{1}{1+x_{3}}\right)\right],$
(12)

where $\Theta(x)$ is the Heaviside step function. This leads to the conditions
tabulated in table 1.

$x_{1}$ | $x_{2}$ | $x_{3}$ |  | Condition on regulators
---|---|---|---|---
$<1$ | $<1$ | $>-1$ |  |  $0\overset{!}{=}0$, i.e. no condition
$<1$ | $<1$ | $<-1$ |  | $0\overset{!}{=}-\tilde{0}_{1}-\tilde{0}_{2}+\tilde{0}_{3}+\tilde{0}_{4}$
$<1$ | $<1$ | $>-1$ |  | $0\overset{!}{=}-\tilde{0}_{2}+\tilde{0}_{3}$
$<1$ | $>1$ | $<-1$ |  | $0\overset{!}{=}-\tilde{0}_{1}+\tilde{0}_{4}$
$>1$ | $>1$ | $>-1$ |  | $0\overset{!}{=}-\tilde{0}_{1}+\tilde{0}_{3}$
$>1$ | $<1$ | $<-1$ |  | $0\overset{!}{=}-\tilde{0}_{2}+\tilde{0}_{4}$
$>1$ | $<1$ | $>-1$ |  | $0\overset{!}{=}-\tilde{0}_{1}-\tilde{0}_{2}$
$>1$ | $>1$ | $<-1$ |  | $0\overset{!}{=}\tilde{0}_{3}+\tilde{0}_{4}$
Table 1: Conditions for the regulators of the four $I(\chi)$ functions in eq.
(9) imposed by eq. (12) in the various kinematic regions.

We can satisfy these by choosing

$\displaystyle\mathrm{i}\tilde{0}_{1}\equiv\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(\frac{x_{1}-x_{2}}{1+x_{3}}\right)=\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})(s_{1}-s_{2})}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right),\quad\mathrm{i}\tilde{0}_{2}\equiv-\mathrm{i}\tilde{0}_{1}\,,$
$\displaystyle\mathrm{i}\tilde{0}_{3}\equiv\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(\frac{x_{1}x_{2}-x_{1}-x_{2}}{1+x_{3}}\right)=\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})(p_{4}^{2}-p_{2}^{2})}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right),\quad\mathrm{i}\tilde{0}_{4}\equiv-\mathrm{i}\tilde{0}_{3}\,.$
(13)

For brevity, the common factor
$\frac{p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}$ can be
dropped in the signum expressions, since only relative signs matter. Putting
everyting together analogous to the calculation of the single off-shell
integral in ref. (Haug2022, ), we find that the non-adjacent double off-shell
scalar box integral is given by

$\displaystyle D_{0}\\!$
$\displaystyle\left(s_{1},s_{2},0,p_{2}^{2},0,p_{4}^{2}\right)=\,\frac{1}{\varepsilon^{2}}\frac{\Gamma(1+\varepsilon)\,\Gamma^{2}(1-\varepsilon)}{\Gamma(1-2\varepsilon)}\,\frac{2}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}$
$\displaystyle\times\left\\{\left[\frac{\mu^{2}}{-s_{1}-\mathrm{i}0}\right]^{\varepsilon}{{}_{2}}\mathrm{F}_{1}\\!\left(1,-\varepsilon,1-\varepsilon;\,-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})s_{1}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}+\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(s_{1}-s_{2}\right)\right)\right.$
$\displaystyle\left.\;\;\,+\left[\frac{\mu^{2}}{-s_{2}-\mathrm{i}0}\right]^{\varepsilon}{{}_{2}}\mathrm{F}_{1}\\!\left(1,-\varepsilon,1-\varepsilon;\,-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})s_{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}+\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(s_{2}-s_{1}\right)\right)\right.$
$\displaystyle\left.\;\;\,-\left[\frac{\mu^{2}}{-p_{2}^{2}-\mathrm{i}0}\right]^{\varepsilon}{{}_{2}}\mathrm{F}_{1}\\!\left(1,-\varepsilon,1-\varepsilon;\,-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})p_{2}^{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}+\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(p_{4}^{2}-p_{2}^{2}\right)\right)\right.$
$\displaystyle\left.\;\;\,-\left[\frac{\mu^{2}}{-p_{4}^{2}-\mathrm{i}0}\right]^{\varepsilon}{{}_{2}}\mathrm{F}_{1}\\!\left(1,-\varepsilon,1-\varepsilon;\,-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})p_{4}^{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}+\mathrm{i}\tilde{0}\,\mathrm{sgn}\\!\left(p_{2}^{2}-p_{4}^{2}\right)\right)\right\\}.$
(14)

This generalizes eq. (2.36) of ref. (Haug2022, ), where the single off-shell
box integral is expressed in terms of three Gauss hypergeometric functions.
Remember that $s_{3}=p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}$. Taking the limit
$p_{2}^{2}\rightarrow 0$ or $p_{4}^{2}\rightarrow 0$ trivially reduces the
result to the single off-shell box integral.333The branch cut of each
hypergeometric function is the same in these limits if we restore the common
factor $\frac{p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}$
in the signum expressions. The result is explicitly symmetric under
simultaneously interchanging $s_{1}\leftrightarrow s_{2}$ and
$p_{2}^{2}\leftrightarrow p_{4}^{2}$, as well as under the interchange
$s_{1}\leftrightarrow p_{2}^{2}$, $s_{2}\leftrightarrow p_{4}^{2}$. The branch
cuts of the individual hypergeometric functions are all spurious, i.e. they
cancel in the sum.

While the hypergeometric representation of the single off-shell integral has
been known in the literature for some time Matsuura1989 ; BernDixonKosower1 ,
a similar representation of the non-adjacent double off-shell box integral in
terms of four Gauss hypergeometric functions was only recently found in ref.
(Tarasov:2019mqy, ) (see eqs. (4.34) and (4.66) of that reference). The main
difference to our result is that by a suitable choice of the explicit
$\mathrm{i}\tilde{0}_{i}$, we eliminated the logarithms which are present in
eq. (4.66) of ref. (Tarasov:2019mqy, ).

##  3 All order $\varepsilon$-expansion

To get rid of the spurious branch cuts introduced by the Gauss hypergeometric
functions, we follow the strategy of section 3 of ref. (Haug2022, ). We use
the single-valued counterparts of ${{}_{2}}\mathrm{F}_{1}$ introduced in eq.
(3.15) of ref. (Haug2022, ),

$\displaystyle\mathfrak{F}(\varepsilon;x)\,\equiv\,1\,+\,\ln\\!\left|\frac{x}{x-1}\right|\sum_{n=1}^{\infty}\frac{(-\varepsilon)^{n}}{n!}\ln^{n-1}|x|\,-\,\sum_{n=2}^{\infty}\varepsilon^{n}\mathcal{L}_{n}\\!\left(x\right)\,,$
(15)

with the single-valued polylogarithms (see Appendix C of ref. (Haug2022, ))
defined by

$\displaystyle\mathcal{L}_{n}\\!\left(x\right)=\sum_{k=0}^{n-1}\frac{(-1)^{k}}{k!}\,\ln^{k}|x|\,\mathrm{Li}_{n-k}\\!\left(x\right)+\frac{(-1)^{n-1}}{n!}\ln^{n-1}|x|\ln|1-x|\,.$
(16)

The function $\mathfrak{F}$ is defined directly by its all-order
$\varepsilon$-expansion. An alternative form useful for kinematic
regularization purposes444For an example see Appendix D of ref. (Haug2022, ),
which resums the $\ln|x|$ terms, reads

$\displaystyle\mathfrak{F}(\varepsilon;x)=|x|^{-\varepsilon}+\varepsilon\ln|1-x|-\sum_{n=2}^{\infty}\varepsilon^{n}\left[\frac{(-1)^{n}\ln|1-x|\ln^{n-1}|x|}{n!}+\mathcal{L}_{n}\\!\left(x\right)\right].$
(17)

Near $x=0$, $\mathfrak{F}(\varepsilon;x)$ behaves like $|x|^{-\varepsilon}$.
While the individual $\mathfrak{F}(\varepsilon;x)$ have a $\ln|1-x|$ term, the
divergence at $x=1$ cancels between the four $\mathfrak{F}(\varepsilon;x)$ in
the complete result.

Expressing the non-adjacent double off-shell scalar box integral in terms of
$\mathfrak{F}$, we find the representation analogous to eq. (3.26) of ref.
(Haug2022, ),

$\displaystyle\\!\\!\\!\\!\\!D_{0}\\!\left(s_{1},s_{2},p_{2}^{2},p_{4}^{2}\right)$
$\displaystyle=\,\frac{1}{\varepsilon^{2}}\frac{\Gamma(1+\varepsilon)\,\Gamma^{2}(1-\varepsilon)}{\Gamma(1-2\varepsilon)}\,\frac{2}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\left|\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})\mu^{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right|^{\varepsilon}$
$\displaystyle\times\left\\{\left(\Theta(-s_{1})+\Theta(s_{1})\mathrm{e}^{\mathrm{i}\pi\varepsilon}\right)\,\mathfrak{F}\\!\left(\varepsilon;-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})s_{1}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right)\right.$
$\displaystyle\phantom{\times}+\left(\Theta(-s_{2})+\Theta(s_{2})\mathrm{e}^{\mathrm{i}\pi\varepsilon}\right)\,\mathfrak{F}\\!\left(\varepsilon;-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})s_{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right)$
$\displaystyle\phantom{\times}-\left(\Theta(-p_{2}^{2})+\Theta(p_{2}^{2})\mathrm{e}^{\mathrm{i}\pi\varepsilon}\right)\,\mathfrak{F}\\!\left(\varepsilon;-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})p_{2}^{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right)$
$\displaystyle\left.\\!\phantom{\times}-\left(\Theta(-p_{4}^{2})+\Theta(p_{4}^{2})\mathrm{e}^{\mathrm{i}\pi\varepsilon}\right)\,\mathfrak{F}\\!\left(\varepsilon;-\frac{(p_{2}^{2}+p_{4}^{2}-s_{1}-s_{2})p_{4}^{2}}{s_{1}s_{2}-p_{2}^{2}p_{4}^{2}}\right)\right\\}.$
(18)

Real and imaginary parts in every kinematic region, as well as the all-order
$\varepsilon$-expansion, can conveniently be read off in this form. Note that
eq. (18) is free of spurious branch cuts due to the use of single-valued
polylogarithms. This representation constitutes the main new result of this
work.

##  4 Conclusion

We have generalized our recent result for the single off-shell scalar box
integral to the case of two non-adjacent external particles off the light cone
and established an all-order $\varepsilon$-expansion with explicit real and
imaginary parts. The solution turns out to be even more symmetric than the
single off-shell case to which it trivially reduces in the limit of either one
particle put on the light cone. However, for most applications the adjacent
double off-shell integral would also be needed. For that integral, a result in
the style found for the non-adjacent case does not seem to be accessible by a
straightforward generalization of the method proposed in ref. (Haug2022, ).
This is related to the lower degree of symmetry in its Feynman parameter
representation. The results of ref. (Tarasov:2019mqy, ) found through
functional equation techniques could be a promising starting point towards a
similar all-order $\varepsilon$-expansion of the adjacent double-off shell as
well as the triple and fully off-shell cases.

###### Acknowledgements.

We thank the anonymous referee of ref. (Haug2022, ) for inspiring this paper
by raising the question whether the result is generalizable to two end points
off the light cone. We also thank Oleg V. Tarasov for drawing our attention to
the similar hypergeometric representation he obtained for the non-adjacent
double off-shell box integral. We are grateful to Werner Vogelsang for helpful
comments. This study was supported in part by Deutsche Forschungsgemeinschaft
(DFG) through the Research Unit FOR 2926 (Project No. 409651613). J. H. is
grateful to the Landesgraduiertenförderung Baden-Württemberg for supporting
her research.

## References

  * (1) J. Haug and F. Wunder, _The massless single off-shell scalar box integral — branch cut structure and all-order epsilon expansion_ , _JHEP_ 02 (2023) 177 [2211.14110].
  * (2) Z. Bern, L. Dixon, and D. A. Kosower, _Dimensionally-regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751.
  * (3) G. Duplančić and B. Nižić, _Dimensionally-regulated one-loop scalar integrals with massless internal lines_ , _Eur. Phys. J. C_ 20 (2001) 357.
  * (4) R. K. Ellis and G. Zanderighi, _Scalar one-loop integrals for QCD_ , _JHEP_ 02 (2008) 002 [arxiV:0712.1851].
  * (5) O. V. Tarasov, _Functional reduction of Feynman integrals_ , _JHEP_ 02 (2019) 173 [1901.09442].
  * (6) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485 [hep-ph/9912329].
  * (7) V.A. Smirnov, _Analytical result for dimensionally regularized massless on shell double box_ , _Phys. Lett. B_ 460 (1999) 397 [hep-ph/9905323].
  * (8) C. Anastasiou, E. W. N. Glover, and C. Oleari, _Application of the negative dimension approach to massless scalar box integrals_ , _Nucl. Phys. B_ 565 (2000) 445 [hep-ph/9907523].
  * (9) J. Fleischer, F. Jegerlehner, and O. V. Tarasov, _A New hypergeometric representation of one loop scalar integrals in d dimensions_ , _Nucl. Phys. B_ 672 (2003) 303 [hep-ph/0307113].
  * (10) O. V. Tarasov, _Functional reduction of one-loop Feynman integrals with arbitrary masses_ , _JHEP_ 06 (2022) 155 [2203.00143].
  * (11) J. P. Ellis, _TikZ-Feynman: Feynman diagrams with TikZ_ , _Comp. Phys. Comm._ 210 (2017) 103 [arXiv:1601.05437].
  * (12) V. A. Smirnov, _Analytic Tools for Feynman Integrals_ , Springer (2012).
  * (13) V. E. Lyubovitskij, F. Wunder and A. S. Zhevlakov, _New ideas for handling of loop and angular integrals in D-dimensions in QCD_ , _JHEP_ 06 (2021) 066 [arXiv:hep-ph/2102.08943].
  * (14) T. Matsuura, S. C. van der Marck, and W. L. Van Neerven, _The calculation of the second order soft and virtual contributions to the Drell-Yan cross section_ , _Nucl. Phys. B_ 319 (1989) 570.