# Covariant Guiding Center Equations for Charged Particle Motions in General
Relativistic Spacetimes

Tyler Trent  School of Physics, Georgia Institute of Technology, 837 State St
NW, Atlanta, GA 30332, USA  Departments of Astronomy and Physics, University
of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA<EMAIL_ADDRESS>Karin Roley  School of Physics, Georgia Institute of Technology, 837 State St
NW, Atlanta, GA 30332, USA  Matthew Golden  School of Physics, Georgia
Institute of Technology, 837 State St NW, Atlanta, GA 30332, USA  Dimitrios
Psaltis  School of Physics, Georgia Institute of Technology, 837 State St NW,
Atlanta, GA 30332, USA  Feryal Özel  School of Physics, Georgia Institute of
Technology, 837 State St NW, Atlanta, GA 30332, USA

###### Abstract

Low density plasmas in curved spacetimes, such as those found in accretion
flows around black holes, are challenging to model from first principles,
owing to the large scale separation between the characteristic scales of the
microscopic processes and large mean-free-paths comparable to the system
sizes. Kinetic approaches become necessary to capture the relevant physics but
lack the dynamic range to model the global characteristics of the systems. In
this paper, we develop new covariant guiding center equations of motion for
charges in general relativistic spacetimes that are computationally tractable.
We decompose the particle motion into a fast gyration, which we integrate
analytically and a slow drift of the guiding center, which can be solved
numerically. We derive covariant conservation laws for the motions of the
guiding centers and show, through a number of limiting cases, that the
equations contain all known drift mechanisms. Finally, we present the general
relativistic expressions for the various drift velocities in Schwarzschild
spacetimes.

##  1 Introduction

Rarefied plasmas are found in a multitude of astrophysical environments
ranging from the heliosphere to the intracluster medium and the accretion
flows around black holes. They are often modeled as fluids using equations of
Magnetohydrodynamics (MHD) or, as needed, of General Relativistic
magnetohydrodynamics (GRMHD). The computational efficiency of a fluid approach
allows for global simulations and has proven useful for understanding the
overall dynamics of the plasma in these systems. As one particular example,
GRMHD simulations have guided the Event Horizon Telescope (EHT) imaging
observations of the black hole at the center of the M87 galaxy and Sagittarius
A* (SgrA*) and enabled the reconstruction and interpretation of these images
(Event Horizon Telescope Collaboration, 2019, 2022; Medeiros et al., 2023).

Despite their usefulness, fluid approaches are not formally valid in the
weakly collisional (or collisionless) conditions of these rarefied plasmas.
The long mean-free-paths of charged particles along field lines are expected
to lead to large pressure anisotropies and introduce non-local couplings of
the dynamics and thermodynamics of the plasmas. This shortcoming of fluid
approaches is especially important when modeling, e.g., the observable effects
of episodic magnetic reconnection events that accelerate electrons to ultra-
relativistic non-thermal energies and are believed to be responsible for the
X-ray flares observed from Sgr A*. Such events have been traditionally modeled
by injecting locally a distribution of ultra-relativistic electrons into the
background of a GRMHD simulation and then calculating their radiative
properties without accounting for their dynamics (see, e.g., Dodds-Eden et
al., 2010; Ball et al., 2016; Chatterjee et al., 2021; Scepi et al., 2022).
This local approach ignores the mean-free paths of the ultra-relativistic
particles along magnetic field lines, which can be significantly larger than
the macroscopic (i.e., horizon) scales of the systems. Perhaps more
importantly, it does not address the apparent tension with imaging
observations of the flare emission that appears to require these
ultrarelativistic particles to remain trapped in compact regions and orbit
around the black hole (GRAVITY Collaboration et al., 2018).

To overcome these shortcomings, kinetic approaches, such as Particle-in-Cell
(PIC), describe collisionless plasmas from first principles by self-
consistently solving for the individual motions of charged particles and the
fields they produce. These approaches are capable of resolving most of the
microphysics that fluid approaches cannot capture, such as the acceleration of
charged particles in magnetic reconnection events (see, e.g., Sironi &
Spitkovsky, 2014; Werner et al., 2018; Sironi & Beloborodov, 2020) and
instabilities that form from pressure anisotropies (Kunz et al., 2014;
Zhdankin et al., 2023). Hybrid particle-MHD methods have also been developed
to model the dynamics of relativistic particles within MHD backgrounds (see,
e.g., Bai et al., 2015; Mignone et al., 2018; van Marle et al., 2018). In
several situations of astrophysical interest, however, the large scale
separation between the kinetic scales, i.e., the radius of gyration around
local magnetic fields (gyroradius) and the plasma skin depth, and the
macroscopic length scales of astrophysical systems, such as the size of the
accretion flow, constrains the dynamic range that can be simulated
numerically. For this reason, current state-of-the-art kinetic and hybrid
simulations of plasmas around compact objects have largely been limited to
local simulations that study, e.g., magnetic reconnection rates and the
particles they accelerate (Werner et al., 2018; Ball et al., 2018, 2019;
Sironi et al., 2021; Hakobyan et al., 2023). At the opposite extreme, global
PIC simulations of black hole accretion disks use nonphysical parameters to
overcome the scale separation (Parfrey et al., 2019; Crinquand et al., 2020,
2021, 2022; Galishnikova et al., 2023).

Fortunately, the scale separation in the gyroradius and macroscopic length
scales puts these systems in a regime where the kinetic approach can be
employed in the the guiding center approximation. The guiding center formalism
decomposes the motion of charged particles in strong magnetic fields into a
fast gyration (sometimes referred to as gyromotion) about the local magnetic
field lines and the slow drift of the guiding center (the center of the
gyromotion, see e.g., Northrop 1963). This is justified in a strong and slowly
varying magnetic field, where the particle undergoes nearly uniform circular
motion with the center of that circle slowly drifting. Earlier work
analytically solved for the guiding center velocities (sometimes referred to
as the drift velocity) in flat spacetimes, which allows for studies of just
the drift motion of the particle without including the trivial and
uninformative gyromotion (see Cary & Brizard, 2009, for a review). Causes of
drifts perpendicular to the magnetic field lines include the presence of an
electric field $\vec{E}$ (the $\vec{E}\times\vec{B}$ drift), the presence of a
gravitational field (the $\vec{g}\times\vec{B}$ drift), and a gradient in the
magnetic field strength (the $\vec{\nabla}B\times\vec{B}$ drift). This
approach has been successfully implemented in a number of simulations where
test particle trajectories are calculated using the guiding center
approximation with a MHD background (see, e.g, Gordovskyy et al. 2014;
Threlfall et al. 2015; Ripperda et al. 2017; Borissov et al. 2020; Gordovskyy
et al. 2023; Mignone et al. 2023). The length scale in which the guiding
center evolves is much greater than the gyromotion, allowing for a coarser
numerical resolution and a significant numerical speed up in the calculation
of the particle motion. This in turn, enables tracking of these particles
throughout the global system, which otherwise would not have been
computationally possible.

To date, there has not been a general relativistic covariant formulation of
the guiding center approach. One attempt by Beklemishev & Tessarotto (2004)
uses the tetrad formalism and leaves the equations in a form that has not been
computationally tractable. In another drift-kinetic study in general
relativistic spacetimes (Bacchini et al., 2020), only the
$\vec{E}\times\vec{B}$ drift velocity was incorporated by using the special-
relativistic expression to evaluate it in a local Lorentz frame and then
boosting it to the coordinate frame. This method does not allow for magnetic
mirroring of charged particles and is not easily generalizable to include
drifts that cannot be readily written in terms of quantities evaluated in a
local Lorentz frame, such as those arising from nonuniform electromagnetic
fields or from gravity (which vanishes in the local Lorentz frame).

Trent et al. (2023) introduced a set of covariant, general relativistic
guiding center equations that can be used in arbitrary spacetimes and for
arbitrary particle velocities. The key insight in that treatment came from the
realization that it is possible to circumvent the difficulty of obtaining
covariant equations for the various drift velocities by deriving, instead, a
covariant acceleration equation for the motion of the gyrocenter. We then used
a set of numerical experiments in different configurations to demonstrate that
the new covariant formalism incorporates all of the drifts described by flat-
spacetime guiding-center equations.

In this paper, we present a detailed analytic exploration of the new covariant
guiding-center equations for charged particles in arbitrary spacetimes. In §2,
we outline the detailed derivation of the acceleration equation for the
guiding center in the presence of an arbitrary electromagnetic field. In §3,
we derive an integral of motion for the guiding center and show that it is an
analog to the unity norm of the 4-velocity of a relativistic particle. In §4,
we consider analytically a number of limiting cases and show that the guiding
center equations recover the standard drift velocities. We then generalize the
expressions for the drift velocities to the case of motion in a Schwarzschild
spacetime. Finally, in §5, we summarize our results and discuss possible
applications of these general relativistic guiding center equations.

##  2 Covariant Guiding Center Equations

In this section we derive the covariant guiding center equation of motion. We
start from the equation of motion for a charged particle in a general
spacetime with an arbitrary electromagnetic field

$\displaystyle\frac{d^{2}x^{\alpha}}{d\tau^{2}}=-\Gamma^{\alpha}_{\mu\nu}\frac{dx^{\mu}}{d\tau}\frac{dx^{\nu}}{d\tau}+\frac{q}{m}F^{\alpha}_{~{}\beta}\frac{dx^{\beta}}{d\tau},$
(1)

where $x^{\alpha}$ is the four-position of the particle, $\tau$ is the proper
time, $\Gamma^{\alpha}_{\mu\nu}$ are the Christoffel symbols, $q/m$ is the
charge-to-mass ratio of the particle, and $F^{\alpha}_{~{}\beta}$ is the
electromagnetic field tensor. Throughout this section, we use units such that
$G=c=1$.

For systems where the electromagnetic field is slowly varying in time and
space and the effects due to the curvature of spacetime are much weaker than
the effects from the electromagnetic force, the particle undergoes a nearly
uniform circular motion with a slow drift. We refer to this circular motion as
the gyromotion and to the motion of the guiding center as the drift. The goal
of the guiding center approximation is to decompose the motion of the particle
perpendicular to the magnetic field into a fast gyromotion and a slow drift of
the guiding center. This allows us to integrate out the gyromotion
analytically and formulate an equation of motion only for the guiding center.

We can carry out this approximation under the following conditions, (extended
from Vandervoort 1960):
1\. The gyroradius $\rho$ is significantly smaller than the characteristic
scale over which the electromagnetic field varies

$\rho\ll|F^{\alpha}_{~{}\beta}|/\left|\frac{\partial
F^{\alpha}_{~{}\beta}}{\partial x^{\mu}}\right|\;.$ (2)

2\. The particle drifts for many gyroperiods before the field changes
considerably,

$\frac{1}{\omega}\left|\frac{\partial\chi^{\nu}}{\partial\tau}\right|\ll|F^{\alpha}_{~{}\beta}|/\left|\frac{\partial
F^{\alpha}_{~{}\beta}}{\partial x^{\mu}}\right|\;,$ (3)

where $\omega$ is the gryofrequency and $d\chi/d\tau$ is the guiding center
four-velocity.
3\. The effect of the spacetime curvature on the motion of the particle is
weaker than that of the electromagnetic field,

$\left|\Gamma^{\alpha}_{\mu\nu}\frac{dx^{\mu}}{d\tau}\frac{dx^{\nu}}{d\tau}\right|\ll\frac{q}{m}\left|F^{\alpha}_{~{}\beta}\frac{dx^{\beta}}{d\tau}\right|\;.$
(4)

In the above expressions, we use the double vertical bar symbol $|~{}|$ to
denote the magnitude of a typical component. The above assumptions require
that the derivatives of the electromagnetic field along the path of the
guiding center and the Christoffel symbols are all considered first order
perturbations to the guiding center motion compared to the effects of the
electromagnetic field tensor, which are of zeroth order.

Because the final equation of motion will only describe the guiding center, we
evaluate all quantities at the position of the guiding center. The distance
between the position of the particle and that of the guiding center is equal
to the gyroradius. We expand the electromagnetic field tensor around the
guiding center to first order along coordinate lines, which involves only a
partial derivative. However, because we assume that the gravity term is much
smaller than the Lorentz force term, we evaluate the Christoffel symbols at
the guiding center. This amounts to a zeroth-order expansion of the
Christoffel symbols, given that the gravity term is already a first order term
in our approximation. The resulting equation becomes

$\displaystyle\underbrace{\frac{d^{2}x^{\alpha}}{d\tau^{2}}}_{\mathcal{O}(0)}=\underbrace{\left.-\Gamma^{\alpha}_{\mu\nu}\right|_{\chi}\frac{dx^{\mu}}{d\tau}\frac{dx^{\nu}}{d\tau}}_{\mathcal{O}(1)}+\underbrace{\frac{q}{m}\left.F^{\alpha}_{~{}\beta}\right|_{\chi}\frac{dx^{\beta}}{d\tau}}_{\mathcal{O}(0)}$
$\displaystyle+\underbrace{\frac{q}{m}\left.\frac{\partial
F^{\alpha}_{~{}\beta}}{\partial
x^{\mu}}\right|_{\chi}\left(x^{\mu}-\chi^{\mu}\right)\frac{dx^{\beta}}{d\tau}}_{\mathcal{O}(1)}.$
(5)

To derive the guiding center equation of motion, we take an iterative approach
by first solving the zeroth order system, which describes the motion of a
charged particle in a constant electromagnetic field and spacetime. The
solution to the zeroth order system then gives us the necessary building
blocks to construct our ansatz for the first order system. The first order
system includes a general spacetime with arbitrary electromagnetic field that
are consistent with equations (2)-(4). Ultimately, we solve for the
acceleration of the guiding center.

###  2.1 Zeroth Order Approximation

To zeroth order, the equation of motion (2) simplifies to that of a charged
particle in a constant spacetime metric with a constant electromagnetic field,
i.e.,

$\frac{d^{2}x^{\alpha}}{d\tau^{2}}=\frac{q}{m}F^{\alpha}_{~{}\beta}\frac{dx^{\beta}}{d\tau}.$
(6)

Since this is a homogeneous equation, the most convenient way to obtain an
analytic solution is by using the eigenvalues and eigenvectors of the tensor
$(q/m)\,F^{\alpha}_{~{}\beta}$ (Vandervoort, 1960). In the next subsection, we
introduce the eigensystem of the electromagnetic field tensor before deriving
the solution to the differential equation.

####  2.1.1 Eigensystem of the Electromagnetic Field Tensor

We denote the eigenvectors of $q/m\,F^{\alpha}_{~{}\beta}$ by
$\\{\sigma^{\beta},\delta^{\beta},\psi^{\beta},\Upsilon^{\beta}\\}$ and their
eigenvalues by $\\{i\omega,-i\omega,\lambda,-\lambda\\}$, respectively. Here,
$\omega$ and $\lambda$ are real numbers, the eigenvectors $\sigma$ and
$\delta$ are complex, and the eigenvectors $\psi$ and $\Upsilon$ are real. We
conveniently write the eigenvalues in terms of the field invariants, as
(Fradkin, 1978)

$\displaystyle\omega$
$\displaystyle=\frac{q}{2m}\biggl{\\{}F^{\alpha\beta}F_{\alpha\beta}+\Bigl{[}(F^{\alpha\beta}F_{\alpha\beta})^{2}+(F^{\alpha\beta}F^{D}_{\alpha\beta})^{2}\Bigr{]}^{1/2}\biggr{\\}}^{1/2},$
(7) $\displaystyle\lambda$
$\displaystyle=\frac{q}{2m}\biggl{\\{}-F^{\alpha\beta}F_{\alpha\beta}+\Bigl{[}(F^{\alpha\beta}F_{\alpha\beta})^{2}+(F^{\alpha\beta}F^{D}_{\alpha\beta})^{2}\Bigr{]}^{1/2}\biggr{\\}}^{1/2},$
(8)

where $F^{D}_{\alpha\beta}$ stands for the dual of $F_{\alpha\beta}$.

We construct the eigenvectors by invoking the Cayley-Hamilton Theorem, which
shows that every square matrix satisfies its own characteristic polynomial,
i.e., that

$\displaystyle(F^{\alpha}_{~{}\beta}-i\omega
I^{\alpha}_{\beta})(F^{\beta}_{~{}\zeta}+i\omega
I^{\beta}_{\zeta})(F^{\zeta}_{~{}\mu}-\lambda
I^{\zeta}_{\mu})(F^{\mu}_{~{}\nu}+\lambda I^{\mu}_{\nu})=0,$ (9)

where $I$ is the four-by-four identity matrix. The non-normalized eigenvector
for a given eigenvalue is then the subset of the product of the three
binomials from the characteristic polynomial that correspond to the other
eigenvalues. We write the non-normalized eigenvectors as

$\displaystyle\tilde{\sigma}^{\alpha}=$
$\displaystyle(F^{\alpha}_{~{}\beta}+i\omega
I^{\alpha}_{\beta})(F^{\beta}_{~{}\mu}-\lambda
I^{\beta}_{\mu})(F^{\mu}_{~{}\nu}+\lambda I^{\mu}_{\nu})e^{\nu},$ (10)
$\displaystyle\tilde{\delta}^{\alpha}=$
$\displaystyle(F^{\alpha}_{~{}\beta}-i\omega
I^{\alpha}_{\beta})(F^{\beta}_{~{}\mu}-\lambda
I^{\beta}_{\mu})(F^{\mu}_{~{}\nu}+\lambda I^{\mu}_{\nu})e^{\nu},$ (11)
$\displaystyle\tilde{\psi}^{\alpha}=$
$\displaystyle(F^{\alpha}_{~{}\beta}-i\omega
I^{\alpha}_{\beta})(F^{\beta}_{~{}\mu}+i\omega
I^{\beta}_{\mu})(F^{\mu}_{~{}\nu}+\lambda I^{\mu}_{\nu})e^{\nu},$ (12)
$\displaystyle\tilde{\Upsilon}^{\alpha}=$
$\displaystyle(F^{\alpha}_{~{}\beta}-i\omega
I^{\alpha}_{\beta})(F^{\beta}_{~{}\mu}+i\omega
I^{\beta}_{\mu})(F^{\mu}_{~{}\nu}-\lambda I^{\mu}_{\nu})e^{\nu},$ (13)

where $e^{\nu}$ is any non-zero four-vector.

When there are no degenerate eigenvalues, we normalize the eigenvectors such
that $\sigma^{\alpha}\delta_{\alpha}=1$ and
$\psi^{\alpha}\Upsilon_{\alpha}=1$. When an eigenvalue is degenerate, we
normalize the corresponding eigenvector to 1 by contracting it with itself
instead.

####  2.1.2 Zeroth-order Solution of Charged Particle Trajectory

To obtain a solution to equation (6), we start by projecting the four-velocity
onto the normalized eigenbasis111If there are degenerate eigenvalues, given
the symmetry of the eigensystem, the eigenvectors can always be constructed to
be linearly independent.. We then insert the expanded four-velocity into
equation (6) and obtain

$\displaystyle\frac{d}{d\tau}\biggl{(}\frac{dx^{\mu}}{d\tau}\delta_{\mu}\sigma^{\alpha}+\frac{dx^{\mu}}{d\tau}\sigma_{\mu}\delta^{\alpha}+\frac{dx^{\mu}}{d\tau}\Upsilon_{\mu}\psi^{\alpha}+\frac{dx^{\mu}}{d\tau}\psi_{\mu}\Upsilon^{\alpha}\biggr{)}=$
$\displaystyle
F^{\alpha}_{~{}\beta}\biggl{(}\frac{dx^{\mu}}{d\tau}\delta_{\mu}\sigma^{\beta}+\frac{dx^{\mu}}{d\tau}\sigma_{\mu}\delta^{\beta}+\frac{dx^{\mu}}{d\tau}\Upsilon_{\mu}\psi^{\beta}+\frac{dx^{\mu}}{d\tau}\psi_{\mu}\Upsilon^{\beta}\biggr{)}.$
(14)

Note that we perform the expansion onto the eigenbasis in the above equation
because the eigenvectors are null vectors.

On the right-hand side, we contract the eigenvectors with the field tensor to
obtain the eigenvalues. On the left-hand-side, we move the eigenvectors
outside of the derivative, since the field tensor is constant to zeroth order,
and obtain

$\displaystyle\frac{d}{d\tau}\biggl{(}\frac{dx^{\mu}}{d\tau}\delta_{\mu}\biggr{)}\sigma^{\alpha}+\frac{d}{d\tau}\biggl{(}\frac{dx^{\mu}}{d\tau}\sigma_{\mu}\biggr{)}\delta^{\alpha}$
$\displaystyle+\frac{d}{d\tau}\biggl{(}\frac{dx^{\mu}}{d\tau}\Upsilon_{\mu}\biggr{)}\psi^{\alpha}+\frac{d}{d\tau}\biggl{(}\frac{dx^{\mu}}{d\tau}\psi_{\mu}\biggr{)}\Upsilon^{\alpha}$
(15)
$\displaystyle=i\omega\biggl{(}\frac{dx^{\mu}}{d\tau}\delta_{\mu}\biggr{)}\sigma^{\alpha}\,-\,i\omega\biggl{(}\frac{dx^{\mu}}{d\tau}\sigma_{\mu}\biggr{)}\delta^{\alpha}~{}+\,\lambda\,\biggl{(}\frac{dx^{\mu}}{d\tau}\Upsilon_{\mu}\biggr{)}\psi^{\alpha}~{}-\,\lambda\,\biggl{(}\frac{dx^{\mu}}{d\tau}\psi_{\mu}\biggr{)}\Upsilon^{\alpha}.$

Using the linear independence of the eigenvectors, we turn this equation into
four separate differential equations, each of which has an exponential
solution. The general solution is the sum of all four, i.e.,

$\displaystyle\frac{dx^{\alpha}}{d\tau}(\tau)=i\omega\rho_{0}e^{i\omega\tau}\sigma^{\alpha}-i\omega\rho_{0}^{*}e^{-i\omega\tau}\delta^{\alpha}$
$\displaystyle+\Upsilon_{\beta}\frac{dx^{\beta}}{d\tau}\bigg{|}_{\tau=0}e^{\lambda\tau}\psi^{\alpha}+\psi_{\beta}\frac{dx^{\beta}}{d\tau}\bigg{|}_{\tau=0}e^{-\lambda\tau}\Upsilon^{\alpha},$
(16) Integrating this once over proper time, we obtain $\displaystyle
x^{\alpha}(\tau)=\underbrace{\rho_{0}e^{i\omega\tau}\sigma^{\alpha}+\rho_{0}^{*}e^{-i\omega\tau}\delta^{\alpha}}_{\text{gyromotion}}$
$\displaystyle+\underbrace{\frac{1}{\lambda}\Upsilon_{\beta}\frac{dx^{\beta}}{d\tau}\bigg{|}_{\tau=0}e^{\lambda\tau}\psi^{\alpha}-\frac{1}{\lambda}\psi_{\beta}\frac{dx^{\beta}}{d\tau}\bigg{|}_{\tau=0}e^{-\lambda\tau}\Upsilon^{\alpha}+\kappa^{\alpha}}_{\text{drift}}.$
(17)

Here, $\kappa^{\alpha}$ is a constant determined by the initial conditions and
we defined the constant

$\rho_{0}\equiv\biggl{(}\frac{i\sigma_{\beta}}{\omega}\frac{dx^{\beta}}{d\tau}\biggr{)}\biggl{|}_{\tau=0}\;.$
(18)

We can interpret the solution in equation (17) as consisting of two parts: a
gyromotion and a drift motion. The oscillatory terms represent the gyromotion,
with $\omega$ as the gyrofrequency and $\sigma$ and $\delta$ as four-vectors
that live in the gyromotion plane. The complex quantity $\rho_{0}$ defined
above has a magnitude that is equal to the gyroradius of motion divided by
$\sqrt{2}$ (because the gyromotion is decomposed into two orthogonal
gyromotions along the eigenvectors $\sigma^{\alpha}$ and $\delta^{\alpha}$)
and a phase that corresponds to the initial azimuthal position of the particle
on the circle of gyration. The exponential terms represent the drift (or
guiding-center) motion at zeroth order.

###  2.2 First Order Approximation

When the requirements in equations (2)-(4) are satisfied, we can decompose the
first-order solution to equation (1) into a gyromotion component, which
remains of zeroth order, and a drifting guiding center component, which will
be of first order. Retaining only the gyromotion terms from the zeroth-order
equation (17) and allowing for the drift term to be arbitrary, we construct
the following ansatz

$\displaystyle
x^{\alpha}(\tau)=\rho(\tau)\sigma^{\alpha}+\rho^{*}(\tau)\delta^{\alpha}+\chi^{\alpha}\;,$
(19)

where $\chi^{\alpha}$ is the four-position of the guiding center. Here,
$\rho(\tau)$ is the instantaneous complex gyroradius. Its magnitude is
determined by the conservation of the magnetic moment (see §3)

$\mu\equiv\rho_{0}^{2}\omega_{0}=\rho^{2}\omega$ (20)

such that

$\rho(\tau)=\rho_{0}\sqrt{\frac{\omega_{0}}{\omega}}e^{-i\omega\tau}\;.$ (21)

Taking the time derivative of the ansatz gives the four-velocity

$\displaystyle\frac{dx^{\alpha}}{d\tau}=$
$\displaystyle\,i(\omega+\frac{d\omega}{d\tau}\tau)\rho\sigma^{\alpha}+i\omega\rho\frac{d\sigma^{\alpha}}{d\tau}$
(22)
$\displaystyle-i(\omega+\frac{d\omega}{d\tau}\tau)\rho^{*}\delta^{\alpha}-i\omega\rho^{*}\frac{d\delta^{\alpha}}{d\tau}+\frac{d\chi}{d\tau}.$

Our second assumption, equation (3), implies that derivatives of the
electromagnetic field along the path of the particle are small. This is
equivalent to assuming that total derivatives with respect to proper time of
the eigenvalues and eigenvectors are of first order in our approximation.
Dropping all terms in equation (22) that are of second order, we obtain

$\displaystyle\frac{dx^{\alpha}}{d\tau}\approx\,i\omega\rho\sigma^{\alpha}-i\omega\rho^{*}\delta^{\alpha}+\frac{d\chi^{\alpha}}{d\tau}.$
(23)

Similarly, we find the four-acceleration to first order as

$\displaystyle\frac{d^{2}x^{\alpha}}{d\tau^{2}}\approx\,-\omega^{2}\rho\sigma^{\alpha}-\omega^{2}\rho^{*}\delta^{\alpha}+\frac{d^{2}\chi^{\alpha}}{d\tau^{2}}.$
(24)

Inserting equation (19) for the position of the particle, equation (23) for
the four-velocity, and equation (24) for the four-acceleration into equation
(2), we find

$\displaystyle-\omega^{2}\rho\sigma^{\alpha}-\omega^{2}\rho^{*}\delta^{\alpha}+\frac{d^{2}\chi^{\alpha}}{d\tau^{2}}=$
$\displaystyle-\Gamma^{\alpha}_{\mu\nu}\biggl{(}i\omega\rho\sigma^{\mu}-i\omega\rho^{*}\delta^{\mu}+\frac{d\chi^{\mu}}{d\tau}\biggr{)}\biggl{(}i\omega\rho\sigma^{\nu}-i\omega\rho^{*}\delta^{\nu}+\frac{d\chi^{\nu}}{d\tau}\biggr{)}$
(25)
$\displaystyle\quad+\frac{q}{m}F^{\alpha}_{~{}\beta}\biggl{(}i\omega\rho\sigma^{\beta}-i\omega\rho^{*}\delta^{\beta}+\frac{d\chi^{\beta}}{d\tau}\biggr{)}+\frac{q}{m}\frac{\partial
F^{\alpha}_{~{}\beta}}{\partial
x^{\mu}}\left(\rho\sigma^{\mu}+\rho^{*}\delta^{\mu}\right)\biggl{(}i\omega\rho\sigma^{\beta}-i\omega\rho^{*}\delta^{\beta}+\frac{d\chi^{\beta}}{d\tau}\biggr{)}.$

We integrate equation (25) over one gyroperiod, which removes all oscillatory
terms. As a result, all terms that are linear in $\rho$ and $\rho^{*}$ or
contain the products $\rho\rho$ and $\rho^{*}\rho^{*}$ vanish, leaving only
terms with $\rho\rho^{*}$ (since $e^{i\omega\tau}e^{-i\omega\tau}=1$) or
$d\chi/d\tau$. The resulting equation is

$\displaystyle\frac{d^{2}\chi^{\alpha}}{d\tau^{2}}=$
$\displaystyle-\Gamma^{\alpha}_{\mu\nu}\biggl{(}\omega\omega_{0}\rho_{0}^{2}\sigma^{\mu}\delta^{\nu}+\omega\omega_{0}\rho_{0}^{2}\sigma^{\nu}\delta^{\mu}+\frac{d\chi^{\mu}}{d\tau}\frac{d\chi^{\nu}}{d\tau}\biggr{)}$
$\displaystyle+\frac{q}{m}F^{\alpha}_{~{}\beta}\frac{d\chi^{\beta}}{d\tau}+i\omega_{0}\rho_{0}^{2}\frac{q}{m}\frac{\partial
F^{\alpha}_{~{}\beta}}{\partial
x^{\mu}}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta}).$ (26)

We can now combine the first two terms and the last term of equation (2.2) to
form a covariant derivative of the field tensor. To show this, we start with
the expression for the covariant derivative of the electromagnetic field
tensor

$\displaystyle\frac{q}{m}\nabla_{\mu}F^{\alpha}_{~{}\beta}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta})=$
$\displaystyle\frac{q}{m}\frac{\partial F^{\alpha}_{~{}\beta}}{\partial
x^{\mu}}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta})+\frac{q}{m}\Gamma^{\alpha}_{\mu\nu}F^{\nu}_{~{}\beta}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta}).$
(27)

In the last term, we contract the $\beta$ index of the eigenvectors with the
field tensor to get the corresponding eigenvalues

$\displaystyle\frac{q}{m}\nabla_{\mu}F^{\alpha}_{~{}\beta}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta})=$
$\displaystyle\frac{q}{m}\frac{\partial F^{\alpha}_{~{}\beta}}{\partial
x^{\mu}}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta})\quad$
$\displaystyle+i\omega\Gamma^{\alpha}_{\mu\nu}(\sigma^{\nu}\delta^{\mu}+\sigma^{\mu}\delta^{\nu}).$
(28)

Comparing equations (2.2) and (2.2) allows us to write the equation of motion
for the guiding center in a manifestly covariant form

$\displaystyle\frac{d^{2}\chi^{\alpha}}{d\tau^{2}}=$
$\displaystyle-\Gamma^{\alpha}_{\mu\nu}\frac{d\chi^{\mu}}{d\tau}\frac{d\chi^{\nu}}{d\tau}+\frac{q}{m}F^{\alpha}_{~{}\beta}\frac{d\chi^{\beta}}{d\tau}$
$\displaystyle+i\omega_{0}\rho_{0}^{2}\frac{q}{m}\nabla_{\mu}F^{\alpha}_{~{}\beta}(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta}).$
(29)

Equation (2.2) is the equation of motion presented in Trent et al. (2023).
Even though the last term has an imaginary coefficient, the complex nature of
the eigenvectors render it real. We can simplify it further and put it in a
form that is manifestly real and will allow us in §3 to derive a conservation
law by using the homogeneous Maxwell equation

$\nabla_{\mu}F_{\alpha\beta}+\nabla_{\alpha}F_{\beta\mu}+\nabla_{\beta}F_{\mu\alpha}=0\;.$
(30)

Substituting this into the last term of equation (2.2), we obtain

$\displaystyle(-\frac{q}{m}\nabla_{\alpha}F_{\beta\mu}-\frac{q}{m}\nabla_{\beta}F_{\mu\alpha})(\sigma^{\beta}\delta^{\mu}-\sigma^{\mu}\delta^{\beta})$
$\displaystyle\qquad=\frac{q}{m}\nabla_{\alpha}F_{\mu\beta}\sigma^{\beta}\delta^{\mu}.$
(31)

The $\alpha$ index on the covariant derivative allows us to fully contract the
field tensor with the eigenvectors by using the product rule:

$\displaystyle\frac{q}{m}\nabla^{\alpha}F^{\mu}_{~{}\beta}\sigma^{\beta}\delta_{\mu}$
$\displaystyle=\frac{q}{m}\nabla^{\alpha}(F^{\mu}_{~{}\beta}\sigma^{\beta})\delta_{\mu}-\frac{q}{m}F^{\mu}_{~{}\beta}\nabla^{\alpha}(\sigma^{\beta})\delta_{\mu}$
$\displaystyle=\nabla^{\alpha}(i\omega\sigma^{\mu})\delta_{\mu}-i\omega\delta^{\beta}\nabla^{\alpha}(\sigma_{\beta})$
$\displaystyle=i\nabla^{\alpha}\omega.$ (32)

Finally, substituting equation (32) into equation (2.2), we obtain

$\displaystyle\frac{d^{2}\chi^{\alpha}}{d\tau^{2}}=$
$\displaystyle-\Gamma^{\alpha}_{\mu\nu}\frac{d\chi^{\mu}}{d\tau}\frac{d\chi^{\nu}}{d\tau}+\frac{q}{m}F^{\alpha}_{~{}\beta}\frac{d\chi^{\beta}}{d\tau}-\omega_{0}\rho_{0}^{2}\nabla^{\alpha}\omega.$
(33)

This equation shows that, to first order, the covariant guiding center
equation of motion differs from the one that describe the full motion of the
particle only by a single term involving the derivative of the scalar
gyrofrequency.

##  3 Conservation Laws

In flat spacetime, integrating the gyromotion of a charged particle results in
the adiabatic conservation of the magnetic moment, usually defined as
$\mu=v_{\perp}^{2}/2B$, where $v_{\perp}$ is the velocity of the charge in the
direction perpendicular to the magnetic field $B$. Fradkin (1978) showed that
a covariant generalization of the definition of the magnetic moment can be
derived as

$\mu=\frac{1}{2m\omega}F^{\alpha}_{\;\beta}{\cal A}^{\beta}_{\;\alpha}\;,$
(34)

where $m$ is the mass of the charge, $\omega$ is the appropriate eigenvalue of
$(q/m)F^{\alpha}_{\;\beta}$, as before, and the antisymmetric electromagnetic
moment is defined as

${\cal
A}^{\alpha\beta}\equiv\frac{1}{2}\frac{q}{c}\left(x^{\alpha}\frac{dx^{\beta}}{d\tau}-x^{\beta}\frac{dx^{\alpha}}{d\tau}\right)\;.$
(35)

In this expression, the position vector $x^{\alpha}$ is understood to include
only the gyration and not the motion of the guiding center.

Writing, in general, the equation of the gyromotion as

$x^{\alpha}=\hat{\rho}(\tau)e^{i\omega\tau}\sigma^{\alpha}+\hat{\rho}^{*}(\tau)e^{-i\omega\tau}\delta^{\alpha}$
(36)

and inserting this expression into equations (34)-(35), we obtain
$\mu=\hat{\rho}^{2}\omega$. This is the origin of equation (21) that we used
before.

The guiding center is a virtual position of a particle and, therefore, does
not have to obey any of the conservation laws that apply to actual particles.
The most trivial conservation law for massive particles is the norm of the
four-velocity, $u^{\alpha}u_{\alpha}=-1$. Interestingly, the guiding center
has an analog to this conservation law, as we derive in this section.

We contract the acceleration equation of the guiding center velocity, equation
(33), with the four-velocity of the guiding center to obtain

$\displaystyle\frac{d\chi_{\alpha}}{d\tau}\frac{d^{2}\chi^{\alpha}}{d\tau^{2}}=$
$\displaystyle-\frac{d\chi_{\alpha}}{d\tau}\Gamma^{\alpha}_{\mu\nu}\frac{d\chi^{\mu}}{d\tau}\frac{d\chi^{\nu}}{d\tau}$
$\displaystyle\quad+\frac{d\chi_{\alpha}}{d\tau}\frac{q}{m}F^{\alpha}_{~{}\beta}\frac{d\chi^{\beta}}{d\tau}-\frac{d\chi_{\alpha}}{d\tau}\omega_{0}\rho_{0}^{2}\nabla^{\alpha}\omega.$
(37)

The term containing the field tensor is zero due to its anti-symmetry. We then
simplify the left hand side using the product rule,

$\displaystyle\frac{1}{2}\frac{d}{d\tau}\biggl{(}\frac{d\chi^{\alpha}}{d\tau}\frac{d\chi_{\alpha}}{d\tau}\biggr{)}=$
$\displaystyle-\frac{d\chi_{\alpha}}{d\tau}\Gamma^{\alpha}_{\mu\nu}\frac{d\chi^{\mu}}{d\tau}\frac{d\chi^{\nu}}{d\tau}-\frac{d\chi_{\alpha}}{d\tau}\omega_{0}\rho_{0}^{2}\nabla^{\alpha}\omega,$
(38)

and combine the left hand side and the term with the Christoffel symbols to
form the directional covariant derivative such that

$\displaystyle\frac{d\chi^{\alpha}}{d\tau}\nabla_{\alpha}\biggl{(}\frac{d\chi^{\beta}}{d\tau}\frac{d\chi_{\beta}}{d\tau}+2\omega_{0}\rho_{0}^{2}\omega\biggr{)}=0.$
(39)

This equation implies that the quantity inside the derivative is constant
along the path of the guiding center, i.e., that

$\frac{d\chi^{\beta}}{d\tau}\frac{d\chi_{\beta}}{d\tau}+2\omega_{0}\rho_{0}^{2}\omega=-1.$
(40)

Using equation (21) we can write this conserved quantity as

$\frac{d\chi^{\beta}}{d\tau}\frac{d\chi_{\beta}}{d\tau}+2\rho^{2}\omega^{2}=-1.$
(41)

Note that we could have obtained the same result by taking directly the norm
of the ansatz for the four-velocity, equation (23).

In the Newtonian limit, the second term in the above equation is just the
magnitude of the gyrating velocity, which initially is the perpendicular
velocity of the particle with respect to the magnetic field. The conservation
law then simply states that the sum of the square of the norm of the guiding
center velocity and the square of the gyrating velocity is equal to $-1$ one
or, more precisely, that it converges to $-1$ to first order in gyroradius.
Surprisingly, this is the same integral of motion found by Vandervoort (1960)
in flat spacetime.

##  4 Limiting Cases

In this section, we examine the behavior of the guiding center motion in
various limiting cases in order to demonstrate explicitly that the covariant
guiding center equations contain all drift velocities that were known
previously. In particular, we consider: (i) a flat spacetime with a non-
uniform electromagnetic field, to recover the $\vec{E}\times\vec{B}$ and
$\vec{\nabla}B\times\vec{B}$ drifts and (ii) a Schwarzschild spacetime to
zeroth Newtonian order to recover the classical gravitational drifts. We then
explore the equation of motion for the full Schwarzschild metric to analyze
the drift velocities of particles in curved spacetime. Throughout this
section, we reintroduce to the equations the various physical constants in
order to compare the drift velocities in curved spacetime with the familiar
drift velocities established in flat spacetime.

###  4.1 Drift Velocities in Flat Spacetime

We first consider the motion of a particle with mass $m$ and charge $q$ in a
flat spacetime under the influence of an electric field $E_{z}$ in the
$\hat{z}$-direction and a magnetic field $B_{x}$ in the $\hat{x}$-direction,
with the latter varying as a function of $z$. We assume that the particle
velocity $u_{p}^{\beta}\equiv dx^{\beta}/d\tau$ at proper time $\tau=0$ is
perpendicular to the direction of the magnetic field, i.e., that $u_{\rm
p}^{x}=0$.

In this case, the electromagnetic tensor is

$F^{\alpha\beta}=\left[\begin{array}[]{cccc}0&0&0&-E_{z}\\\ 0&0&0&0\\\
0&0&0&-B_{x}\\\ E_{z}&0&B_{x}&0\end{array}\right]$ (42)

and the gyro-frequency $\omega_{0}$ is (see eq. [7])

$\displaystyle\omega_{0}=\frac{qB_{x}}{mc}\left(1-\frac{E_{z}^{2}}{B_{x}^{2}}\right)^{1/2}\;,$
(43)

which is the familiar result. To calculate the generalized radius of gyration
$\rho_{0}$, we first evaluate the corresponding eigenvector of the
electromagnetic field using equation (9), i.e.,

$\displaystyle\sigma^{\alpha}=\frac{1}{\sqrt{2}}\biggl{(}\frac{-iE_{z}/B_{x}}{\sqrt{1-E_{z}^{2}/B_{x}^{2}}},0,\frac{-i}{\sqrt{1-E_{z}^{2}/B_{x}^{2}}},1\biggr{)},$
(44)

and then evaluate equation (18) to obtain

$\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\frac{mc}{\sqrt{2}qB_{x}(1-E_{z}^{2}/B_{x}^{2})}\times$ (45)
$\displaystyle\qquad\left[-\frac{E_{z}}{B_{x}}u_{p}^{t}+u_{p}^{y}\right.\left.+i\left(1-\frac{E_{z}^{2}}{B_{x}^{2}}\right)^{1/2}u_{p}^{z}\right].$

Denoting now the 4-velocity of the guiding center by $u^{\beta}$, we can use
equation (33) to write the guiding-center acceleration equations as

$\displaystyle\frac{du^{t}}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{q}{m\,c}\,E_{z}\,u^{z},$ (46)
$\displaystyle\frac{du^{x}}{d\tau}$ $\displaystyle=$ $\displaystyle 0,$ (47)
$\displaystyle\frac{du^{y}}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{q}{c\,m}B_{x}\,u^{z},$ (48)
$\displaystyle\frac{du^{z}}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{q}{m\,c}(E_{z}\,u^{t}-B_{x}\,u^{y})$ (49)
$\displaystyle\qquad-\frac{q}{m\,c}\frac{\rho_{0}^{2}\,\omega_{0}\,}{\sqrt{1-E_{z}^{2}/B_{x}^{2}}}\frac{\partial
B}{\partial z}\;.$

Following the traditional approach, we assume that the acceleration terms on
the left-hand-sides of these equations are of higher order and can, therefore,
be neglected. Equations (46)-(48) immediately lead to $u^{z}=0$. Equation (49)
can then be solved for the only non-zero component of the drift velocity of
the guiding center, $u^{y}$.

In order to recover the $\vec{E}\times\vec{B}$ drift velocity, we will first
consider the case in which $\partial B/\partial z$ can also be neglected.
Under these conditions, equation (49) simplifies to

$\displaystyle 0=\frac{q}{mc}(E_{z}\,u^{t}-B_{x}\,u^{y}).$ (50)

For non-relativistic particle velocities, $u_{\rm p}^{t}=c$. Solving equation
(50) for $u^{y}$ recovers the standard expression for the
$\vec{E}\times\vec{B}$ drift in flat spacetime, i.e.,

$\displaystyle
u^{y}=c\frac{E_{z}}{B_{x}}=\frac{c|\vec{E}\times\vec{B}|}{B^{2}}.$ (51)

Returning to equation (49), we now reintroduce the terms that allow us to
recover the $\vec{\nabla}B\times\vec{B}$ drift. As described in Vandervoort
(1960), this drift has a correction term for particles traveling at
relativistic speeds under the influence of a perpendicular electric field. As
a result, for this particular case, we assume that the 3-velocity of the
particle speed may be comparable to the speed of light. In order to simplify
the expression, we also assume that the particle velocity is perpendicular to
both the electric and magnetic fields, i.e., that $u_{p}^{x}=u_{p}^{y}=0$.

As before, we examine equation (49) after neglecting the acceleration term,
i.e., we set

$\displaystyle
0=\frac{q}{mc}(E_{z}u^{t}-B_{x}u^{y})-\frac{q}{m\,c}\frac{\rho_{0}^{2}\,\omega_{0}\,}{\sqrt{1-E_{z}^{2}/B_{x}^{2}}}\frac{\partial
B}{\partial z}\;.$ (52)

Using equations (43), (45), and (52), we can solve for the guiding center
velocity in the $\hat{y}$-direction. Note that all velocities in this
expression are in terms of the proper time $\tau$. In order to put this
expression in the traditional form, we write the particle 3-velocity
perpendicular to the magnetic field as

$v_{\perp}=\frac{dz}{dt}=\frac{dz}{d\tau}\frac{d\tau}{dt}=u_{p}^{z}\left(\frac{u_{p}^{t}}{c}\right)=\gamma
u_{p}^{z}\;,$ (53)

where $\gamma$ is the Lorentz factor corresponding to the particle velocity.
Similarly, expressing the drift velocity in terms of coordinate time, we write
$u^{y}=\gamma v^{y}$. The net result is

$\displaystyle v^{y}=$ $\displaystyle-\frac{cE_{z}}{B_{x}}+\frac{\gamma
mcv_{\perp}^{2}}{2q}\frac{E_{z}^{2}}{B_{x}^{4}}\nabla
B(1-\frac{E_{z}^{2}}{B_{x}^{2}})^{-2}$ $\displaystyle+\frac{\gamma
mcv_{\perp}^{2}}{2q}\frac{\nabla B}{2B_{x}(1-\frac{E_{z}^{2}}{B_{x}^{2}})}\;.$
(54)

The first term in the equation above is the $\vec{E}\times\vec{B}$ drift we
calculated previously. The second and third terms describe the
$\vec{\nabla}B\times\vec{B}$ drift velocity, along with their relativistic
correction terms. Simplifying the second and third terms further, we recover
the familiar form, i.e.,

$\displaystyle v^{y}$ $\displaystyle=-c\frac{|\vec{E}\times\vec{B}|}{B^{2}}$
$\displaystyle\qquad+\frac{\gamma
mcu_{\perp}^{2}}{2qB(1-\frac{E_{z}^{2}}{B_{x}^{2}})}\left|\vec{B}\times\nabla\biggl{[}B\biggl{(}1-\frac{E_{z}^{2}}{B_{x}^{2}}\biggr{)}\biggr{]}\right|\;.$
(55)

as also presented in Vandervoort (1960).

###  4.2 Newtonian Gravitational Drift

Next, we examine the drift introduced by a central gravitational field, in the
Newtonian limit, in the presence of a constant magnetic field in the
$\hat{\theta}$-direction. To simplify the algebra, we assume that there is no
electric field. To describe the gravitational field, we use the Schwarzschild
metric in the zeroth Newtonian order

$ds^{2}=-(1-2GM/c^{2}r)dt^{2}+dr^{2}+r^{2}(d\theta^{2}+\sin{\theta}^{2}d\phi^{2})\;.$
(56)

We will also assume that the particle is non-relativistic and traveling on the
equatorial plane at $\theta=\pi/2$.

The gyro-frequency of such a system (see eq. [7]) is similar to the one found
in flat spacetime, but without the electric field contributions, i.e.,

$\displaystyle\omega_{0}=\frac{qB_{\theta}}{mc}\;.$ (57)

The corresponding eigenvector is

$\displaystyle\sigma^{\alpha}=\frac{1}{\sqrt{2}}\biggl{(}0,\,i,\,0,\,\frac{1}{r}\biggr{)},$
(58)

and the complex gyroradius is (see eq. [18])

$\displaystyle\rho_{0}=\frac{mc}{\sqrt{2}qB_{\theta}}(-u_{p}^{r}+iru_{p}^{\phi})\;.$
(59)

Evaluating equation (33) in this configuration, we obtain the following set of
equations for the acceleration of the guiding center:

$\displaystyle\frac{du^{t}}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{-1}{(1-\frac{2GM}{rc^{2}})}\frac{2GM}{r^{2}c^{2}}u^{r}u^{t},$
(60) $\displaystyle\frac{du^{r}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{q}{m\,c}B_{\theta}r\,u^{\phi}-\frac{GM}{c^{2}r^{2}}(u^{t})^{2}$
(61) $\displaystyle\qquad+r\,(u^{\theta})^{2}+r\,(u^{\phi})^{2},$
$\displaystyle\frac{du^{\theta}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{2}{r}u^{r}u^{\theta},$ (62)
$\displaystyle\frac{du^{\phi}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{2}{r}u^{r}u^{\phi}+\frac{q}{m\,c}\frac{B_{\theta}}{r}u^{r}.$
(63)

As before, we neglect the acceleration terms such that equations (60)-(62)
result in $u^{r}=0$. We then set $u^{t}=c$ in equation (61) and solve it for
the linear velocity of the guiding center in the $\hat{\phi}$-direction to
find

$\displaystyle
ru^{\phi}=-\frac{m\,c}{qB_{\theta}}\frac{GM}{r^{2}}=\frac{m\,c}{q}\frac{|\vec{g}\times\vec{B}|}{B^{2}}\;.$
(64)

In the last term we have set $\vec{g}\equiv-(GM/r^{2})\hat{r}$, recovering the
standard expression for the gravitational drift velocity.

###  4.3 Gravitational Drift in the Schwarzschild Spacetime

We now consider relativistic charged particle motion in the Schwarzschild
spacetime

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle-(1-2GM/rc^{2})dt^{2}+(1-2GM/rc^{2})^{-1}dr^{2}$ (65)
$\displaystyle\qquad+r^{2}d\theta^{2}+r^{2}\sin{\theta}^{2}d\phi^{2}\;,$

in the presence of a constant magnetic field in the $\hat{\theta}$-direction.
As before, we assume that there is no electric field.

The gyro-frequency for such a system is given by eq. (7), i.e.,

$\displaystyle\omega_{0}=\sqrt{1-\frac{2GM}{rc^{2}}}\frac{qB_{\theta}}{mc}.$
(66)

This is analogous to the gyrofrequency in the Newtonian limit (cf. eq. [57])
but with the addition of an associated redshift factor. Similarly, the
associated eigenvector is

$\displaystyle\sigma^{\alpha}=\frac{1}{\sqrt{2}}\Biggl{(}0,\,-i\sqrt{1-\frac{2GM}{c^{2}r}},\,0,\,\frac{1}{r}\Biggr{)}$
(67)

and the complex gyro-radius $\rho_{0}$ is (see eq. [18])

$\displaystyle\rho_{0}=\frac{1}{\sqrt{1-\frac{2GM}{c^{2}r}}}\frac{mc}{\sqrt{2}qB_{\theta}}\biggl{(}\frac{1}{\sqrt{1-\frac{2GM}{c^{2}r}}}u_{p}^{r}-iru_{p}^{\phi}\biggr{)}\;.$
(68)

The acceleration equation (33) for this configuration becomes

$\displaystyle\frac{du^{r}}{d\tau}=$
$\displaystyle-\frac{q}{m\,c}B_{\theta}r\biggl{(}1-\frac{2GM}{c^{2}r}\biggr{)}u^{\phi}$
$\displaystyle-\frac{GM}{c^{2}r^{2}}\biggl{(}1-\frac{2GM}{c^{2}r}\biggr{)}(u^{t})^{2}-\frac{GM}{c^{2}r^{2}}\frac{1}{1-\frac{2GM}{c^{2}r}}(u^{r})^{2}$
$\displaystyle\qquad\biggl{(}1-\frac{2GM}{c^{2}r}\biggr{)}r(u^{\theta})^{2}+\biggl{(}1-\frac{2GM}{c^{2}r}\biggr{)}r(u^{\phi})^{2}$
$\displaystyle\qquad\qquad-\frac{1}{2}\frac{GM}{c^{2}r}\biggl{[}\frac{(u_{p}^{r})^{2}}{(1-\frac{2GM}{rc^{2}})}+(r\,u_{p}^{\phi})^{2}\biggr{]}\;,$
(69)

where, for brevity, we only display the important $r-$component. As in the
Newtonian case, the other components of the acceleration equation lead to
$u^{r}=0$. Setting the radial acceleration also equal to zero and solving the
resulting equation to first order in the linear velocity in the
$\hat{\phi}$-direction, we obtain

$\displaystyle
ru^{\phi}=-\frac{m\,c}{qB_{\theta}}\frac{GM}{c^{2}r^{2}}(u^{t})^{2}=\frac{m}{qc}\frac{|\vec{g}\times\vec{B}|}{B^{2}}(u^{t})^{2}.$
(70)

Here, we again defined $\vec{g}\equiv-(GM/r^{2})\hat{r}$, for direct
comparison to the Newtonian case, even though this is not the effective
“gravitational” acceleration experienced by the particle in the Schwarzschild
spacetime.

We now use the conservation law (40) to write

$\displaystyle
u^{\beta}u_{\beta}+\biggl{(}-\frac{1}{\sqrt{1-\frac{2GM}{rc^{2}}}}u_{p}^{r}+iru_{p}^{\phi}\biggr{)}^{2}$
$\displaystyle=$ $\displaystyle-c^{2}\Rightarrow$ $\displaystyle
u^{t}u_{t}+u^{i}u_{j}+\left(-\sqrt{g_{rr}}u_{rmp}^{r}+i\sqrt{g_{\phi\phi}}u^{\phi}_{\rm
p}\right)^{2}$ $\displaystyle=$ $\displaystyle-c^{2}\;,$ (71)

where, in the last equation, we used $g_{rr}$ and $g_{\phi\phi}$ to denote
explicitly some of the components of the Schwarzschild metric. Solving this
equation for the time component of the guiding center velocity, we find

$\displaystyle u^{t}=\frac{\gamma c}{\sqrt{1-\frac{2GM}{c^{2}r}}}\;,$ (72)

where

$\gamma=\left(1-u^{i}u_{i}-g_{rr}u_{p}^{r}u_{p}^{r}-g_{\phi\phi}u_{p}^{\phi}u_{p}^{\phi}\right)^{1/2}\;.$
(73)

Here the $i-$index is running only over the spatial components of the
4-vectors. The 3-velocity of the guiding center $u^{i}$ will have a component
along the magnetic field lines, which is equal to the same component of the
particle velocity, i.e., $u^{\theta}=u^{\theta}_{p}$ in this case, and
therefore, of zeroth order. It will also have a component perpendicular to the
magnetic field lines, i.e., the drift velocity, which is of first order.
Because we will be inserting the expression for $u^{t}$ into equation (70),
which is already of first order, we keep only the zeroth order terms in
equation (73) to write

$\gamma=\left(1-g_{ij}u_{p}^{i}u_{p}^{i}\right)^{1/2}\;.$ (74)

This is simply the Lorentz factor associated with the particle velocity.

Substituting these into equation (70), we obtain the gravitational drift
velocity equation in the Schwarzschild spacetime

$\displaystyle
ru^{\phi}_{\vec{g}\times\vec{B}}=\frac{m\,c}{q}\frac{\gamma^{2}}{(1-\frac{2GM}{c^{2}r})}\frac{|\vec{g}\times\vec{B}|}{{B}^{2}}.$
(75)

In order to convert this into an equation for the 3-velocity of the guiding
center written in terms of coordinate time, we use equation (72) to obtain

$\displaystyle
rv^{\phi}_{\vec{g}\times\vec{B}}=ru^{\phi}_{\vec{g}\times\vec{B}}\left(\frac{u^{t}}{c}\right)=\frac{m\,c}{q}\frac{\gamma}{(1-\frac{2GM}{c^{2}r})^{1/2}}\frac{|\vec{g}\times\vec{B}|}{{B}^{2}}.$
(76)

##  5 Conclusions

GRMHD has been the method of choice for simulating the global plasma dynamics
around radiatively inefficient accreting black holes, despite the very large
mean-free paths of particles in these systems. Kinetic approaches, on the
other hand, are able to describe collisionless plasmas from first principles
but are limited in the dynamic range they can simulate due to the large scale
separation between the microscopic kinetic scales, i.e., the gyroradius and
the plasma skin depth, and the macroscopic length scales of the systems. In
this paper, we presented a set of covariant, guiding center equations of
motion applicable to charged particles in general spacetimes and arbitrary
electromagnetic fields. We derived conservation laws for the guiding center
motion and showed that, in various limits, the equations simplify to the
traditional expressions for the drift velocities.

Assuming that the electric field is sufficiently weaker than the magnetic
field, we also derived the fully relativistic expressions for the various
drift velocities in the Schwarzschild spacetime. For brevity, we present here
only the final result

$\displaystyle ru^{\phi}_{\vec{g}\times\vec{B}}$ $\displaystyle=$
$\displaystyle\frac{m\,c}{q}\frac{\gamma^{2}}{(1-\frac{2GM}{c^{2}r})}\frac{|\vec{g}\times\vec{B}|}{{B}^{2}},$
(77) $\displaystyle ru^{\phi}_{\nabla\vec{B}\times\vec{B}}$ $\displaystyle=$
$\displaystyle\frac{m\,c\,v_{\perp}^{2}}{2qB}\frac{\gamma^{2}}{(1-\frac{2GM}{c^{2}r})}\frac{|\nabla
B\times\vec{B}|}{B^{2}},$ (78) $\displaystyle
ru^{\phi}_{\vec{E}\times\vec{B}}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{(1-\frac{2GM}{c^{2}r})^{1/2}}\frac{c|\vec{E}\times\vec{B}|}{B^{2}}.$
(79)

As expected, the drift velocities in the Schwarzschild spacetime follow
similar relations with the counterparts in flat spacetimes, with the addition
of the associated redshift and Lorentz factors.

Finally, writing all equations for the drift velocities in terms of the
corresponding 3-velocities using equation (72), we obtain

$\displaystyle rv^{\phi}_{\vec{g}\times\vec{B}}$ $\displaystyle=$
$\displaystyle\frac{m\,c}{q}\frac{\gamma}{(1-\frac{2GM}{c^{2}r})^{1/2}}\frac{|\vec{g}\times\vec{B}|}{{B}^{2}},$
(80) $\displaystyle rv^{\phi}_{\nabla\vec{B}\times\vec{B}}$ $\displaystyle=$
$\displaystyle\frac{m\,c\,v_{\perp}^{2}}{2qB}\frac{\gamma}{(1-\frac{2GM}{c^{2}r})^{1/2}}\frac{|\nabla
B\times\vec{B}|}{B^{2}},$ (81) $\displaystyle
rv^{\phi}_{\vec{E}\times\vec{B}}$ $\displaystyle=$
$\displaystyle\frac{c|\vec{E}\times\vec{B}|}{B^{2}}.$ (82)

It is interesting to note that the redshift factor appears in the
gravitational drift and in the $\nabla B\times\vec{B}$ drift but not in the
$\vec{E}\times\vec{B}$ drift. In the first case, the redshift factor corrects
the Newtonian gravitational acceleration so that the latter becomes infinite
at the horizon. In the second case, the redshift factor corrects the spatial
derivative of the magnetic field strength, which needs to be expressed in
terms of proper and not coordinate distances. Since the $\vec{E}\times\vec{B}$
drift depends on neither the gravitational acceleration nor on derivatives, it
does not contain any redshift corrections.

This work has been supported by NSF PIRE award OISE-1743747. T.T. acknowledges
support from the Alfred P. Sloan Foundation and the Ford Foundation.

## References

  * Bacchini et al. (2020) Bacchini, F., Ripperda, B., Philippov, A. A., & Parfrey, K. 2020, The Astrophysical Journal Supplement Series, 251, 10, doi: 10.3847/1538-4365/abb604
  * Bai et al. (2015) Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55, doi: 10.1088/0004-637X/809/1/55
  * Ball et al. (2018) Ball, D., Sironi, L., & Özel, F. 2018, The Astrophysical Journal, 862, 80, doi: 10.3847/1538-4357/aac820
  * Ball et al. (2019) —. 2019, The Astrophysical Journal, 884, 57, doi: 10.3847/1538-4357/ab3f2e
  * Ball et al. (2016) Ball, D., Özel, F., Psaltis, D., & kwan Chan, C. 2016, The Astrophysical Journal, 826, 77, doi: 10.3847/0004-637X/826/1/77
  * Beklemishev & Tessarotto (2004) Beklemishev, A., & Tessarotto, M. 2004, A&A, 428, 1, doi: 10.1051/0004-6361:20034208
  * Borissov et al. (2020) Borissov, A., Neukirch, T., Kontar, E. P., Threlfall, J., & Parnell, C. E. 2020, A&A, 635, A63, doi: 10.1051/0004-6361/201936977
  * Cary & Brizard (2009) Cary, J. R., & Brizard, A. J. 2009, Reviews of Modern Physics, 81, 693, doi: 10.1103/RevModPhys.81.693
  * Chatterjee et al. (2021) Chatterjee, K., Markoff, S., Neilsen, J., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 5281, doi: 10.1093/mnras/stab2466
  * Crinquand et al. (2022) Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. 2022, Phys. Rev. Lett., 129, 205101, doi: 10.1103/PhysRevLett.129.205101
  * Crinquand et al. (2020) Crinquand, B., Cerutti, B., Philippov, A., Parfrey, K., & Dubus, G. 2020, Phys. Rev. Lett., 124, 145101, doi: 10.1103/PhysRevLett.124.145101
  * Crinquand et al. (2021) Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. 2021, A&A, 650, A163, doi: 10.1051/0004-6361/202040158
  * Dodds-Eden et al. (2010) Dodds-Eden, K., Sharma, P., Quataert, E., et al. 2010, The Astrophysical Journal, 725, 450, doi: 10.1088/0004-637X/725/1/450
  * Event Horizon Telescope Collaboration (2019) Event Horizon Telescope Collaboration. 2019, The Astrophysical Journal Letters, 875, L1, doi: 10.3847/2041-8213/ab0ec7
  * Event Horizon Telescope Collaboration (2022) —. 2022, The Astrophysical Journal Letters, 930, L12, doi: 10.3847/2041-8213/ac6674
  * Fradkin (1978) Fradkin, D. M. 1978, Journal of Physics A: Mathematical and General, 11, 1069, doi: 10.1088/0305-4470/11/6/010
  * Galishnikova et al. (2023) Galishnikova, A., Philippov, A., Quataert, E., et al. 2023, Phys. Rev. Lett., 130, 115201, doi: 10.1103/PhysRevLett.130.115201
  * Gordovskyy et al. (2014) Gordovskyy, M., Browning, P. K., Kontar, E. P., & Bian, N. H. 2014, A&A, 561, A72, doi: 10.1051/0004-6361/201321715
  * Gordovskyy et al. (2023) Gordovskyy, M., Browning, P. K., Kusano, K., Inoue, S., & Vekstein, G. E. 2023, The Astrophysical Journal, 952, 75, doi: 10.3847/1538-4357/acdb4d
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 618, L10, doi: 10.1051/0004-6361/201834294
  * Hakobyan et al. (2023) Hakobyan, H., Ripperda, B., & Philippov, A. A. 2023, The Astrophysical Journal Letters, 943, L29, doi: 10.3847/2041-8213/acb264
  * Kunz et al. (2014) Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014, Phys. Rev. Lett., 112, 205003, doi: 10.1103/PhysRevLett.112.205003
  * Medeiros et al. (2023) Medeiros, L., Psaltis, D., Lauer, T. R., & Özel, F. 2023, The Astrophysical Journal Letters, 947, L7, doi: 10.3847/2041-8213/acc32d
  * Mignone et al. (2018) Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, ApJ, 859, 13, doi: 10.3847/1538-4357/aabccd
  * Mignone et al. (2023) Mignone, A., Haudemand, H., & Puzzoni, E. 2023, Computer Physics Communications, 285, 108625, doi: 10.1016/j.cpc.2022.108625
  * Northrop (1963) Northrop, T. G. 1963, The Adiabatic Motion of Charged Particles (Interscience Publishers).  https://books.google.com/books?id=pPNQAAAAMAAJ&printsec=frontcover&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false
  * Parfrey et al. (2019) Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101, doi: 10.1103/PhysRevLett.122.035101
  * Ripperda et al. (2017) Ripperda, B., Porth, O., Xia, C., & Keppens, R. 2017, Monthly Notices of the Royal Astronomical Society, 471, 3465, doi: 10.1093/mnras/stx1875
  * Scepi et al. (2022) Scepi, N., Dexter, J., & Begelman, M. C. 2022, Monthly Notices of the Royal Astronomical Society, 511, 3536, doi: 10.1093/mnras/stac337
  * Sironi & Beloborodov (2020) Sironi, L., & Beloborodov, A. M. 2020, The Astrophysical Journal, 899, 52, doi: 10.3847/1538-4357/aba622
  * Sironi et al. (2021) Sironi, L., Rowan, M. E., & Narayan, R. 2021, The Astrophysical Journal Letters, 907, L44, doi: 10.3847/2041-8213/abd9bc
  * Sironi & Spitkovsky (2014) Sironi, L., & Spitkovsky, A. 2014, The Astrophysical Journal Letters, 783, L21, doi: 10.1088/2041-8205/783/1/L21
  * Threlfall et al. (2015) Threlfall, J., Neukirch, T., Parnell, C. E., & Eradat Oskoui, S. 2015, A&A, 574, A7, doi: 10.1051/0004-6361/201424366
  * Trent et al. (2023) Trent, T., Christian, P., Chan, C.-k., Psaltis, D., & Özel, F. 2023, ApJ, 959, L6, doi: 10.3847/2041-8213/acf8c8
  * van Marle et al. (2018) van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394, doi: 10.1093/mnras/stx2509
  * Vandervoort (1960) Vandervoort, P. O. 1960, Annals of Physics, 10, 401, doi: https://doi.org/10.1016/0003-4916(60)90004-X
  * Werner et al. (2018) Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840, doi: 10.1093/mnras/stx2530
  * Zhdankin et al. (2023) Zhdankin, V., Kunz, M. W., & Uzdensky, D. A. 2023, The Astrophysical Journal, 944, 24, doi: 10.3847/1538-4357/acaf54