# Populations of Magnetized Filaments in the
Intracluster Medium and the Galactic Center

F. Yusef-Zadeh Dept Physics and Astronomy, CIERA, Northwestern University,
2145 Sheridan Road, Evanston , IL 60207, USA<EMAIL_ADDRESS> R. G.
Arendt Code 665, NASA/GSFC, 8800 Greenbelt Road, Greenbelt, MD 20771, USA
UMBC/CRESST 2<EMAIL_ADDRESS> M. Wardle School of Mathematical
and Physical Sciences, Research Centre for Astronomy, Astrophysics and
Astrophotonics, Macquarie University, Sydney NSW 2109, Australia,
<EMAIL_ADDRESS>
###### Abstract

Magnetized radio filaments are found in abundance in the inner few hundred pc
of our Galaxy. Progress in understanding this population of filaments has been
slow over the last few decades, in part due to a lack of detection elsewhere
in the Galaxy or in external galaxies. Recent highly sensitive radio continuum
observations of radio galaxies in galaxy clusters have revealed remarkable
isolated filamentary structures in the ICM that are linked to radio jets,
tails and lobes. The origin of this class of filaments is not understood
either. Here, we argue that the underlying physical mechanisms responsible for
the creation of the two populations are the same because of their similarity
in morphology, spacing between the filaments, aspect ratio, magnetic energy
densities to thermal pressure of the medium, and that both populations have
undergone synchrotron aging. These similarities provide an opportunity to
investigate the physical processes in the ISM and ICM for the first time. We
consider that the origin of the filaments in both the GC and ICM is a result
of the interaction of a large-scale wind with clouds, or the filaments arise
through the stretching and collection of field lines by turbulence in weakly
magnetized medium. We examine these ideas toward four radio galaxy filaments
associated with four radio galaxies IC 40B, IC 4496, J1333–3141, ESO137-006
and argue that much can be understood in the future by comparing these two
populations of filaments.

##  1 Introduction

In the last few years, sensitive high resolution radio continuum observations
of galaxy clusters have begun to reveal a population of magnetized filamentary
structures in their intra-cluster medium (ICM) (Shimwell et al., 2016;
Ramatsoku et al., 2020; van Weeren et al., 2021; Fanaroff et al., 2021; Condon
et al., 2021; Brienza et al., 2022; Giacintucci et al., 2022; Knowles et al.,
2022; Rudnick et al., 2022). Discoveries of new structures in the ICM
surrounding radio galaxies (RGs) are in large part due to recent advances in
broad-band correlators installed on radio telescopes, providing remarkable
brightness sensitivity and dynamic range with exquisite details at $\sim
6^{\prime\prime}$ spatial resolution. A striking aspect of the appearance of
these long and narrow isolated filaments is that they are not located within
the radio lobes, the tails of radio galaxies or associated with the radio core
of the Galaxy (Lane et al., 2002; Hardcastle et al., 2019). The filaments are
detected outside radio galaxies but linked to radio lobes and jets with a wide
range of angles to the orientation of radio jets or the lobes. The nature of
this new phenomenon, which appears to be common in high dynamic range images
of radio galaxies, had not been appreciated or predicted in previous studies
due to lack of sufficient surface brightness sensitivity and angular
resolution, particularly at lower frequencies.

Closer to home, in the vastly different environment of the Galactic center
(GC), radio observations have uncovered a population of magnetized radio
filaments with linearly polarized synchrotron emission tracing cosmic ray
activity throughout the inner few hundred pc of the Galaxy (Yusef-Zadeh,
Morris & Chance, 1984; Liszt, 1985; Bally & Yusef-Zadeh, 1989; Gray, et al.,
1991; Haynes, et al., 1992; Staguhn, et al., 1998; Lang, Morris & Echevarria,
1999; LaRosa, Lazio & Kassim, 2001; LaRosa, et al., 2004; Yusef-Zadeh et al.,
2004; Nord et al., 2004; Law, Yusef-Zadeh & Cotton, 2008; Pound & Yusef-Zadeh,
2018; Staguhn et al., 2019; Arendt et al., 2019). MeerKAT observations have
provided a remarkable mosaic of the inner few degrees of the Galactic center
revealing hundreds of radio filaments, housed within a $\sim 400$pc bipolar
Radio Bubble filled with thermal X-ray gas that surrounds the Galactic center
(Heywood, et al., 2019; Heywood et al., 2022; Yusef-Zadeh et al., 2022a, b,
c). $H_{3}^{+}$ measurements toward the Galactic center have quantifield the
cosmic ray ionization rate to be 102 to $10^{3}$ times greater than across the
rest of the Galaxy (Oka, et al., 2005; Oka et al., 2019). This high ionization
rate led to the interpretion that cosmic-rays drive a wind, producing the
Radio Bubble (Yusef-Zadeh & Wardle, 2019).

In spite of decades of studying the GC filaments, a good understanding of
these structures is still lacking. The emerging population of extragalactic
radio filaments appears to have remarkably similar properties to the radio
filaments in the Galactic center. The comparison of generic properties of both
populations provide an opportunity to gain insight into the origin of both
populations. Here, side-by-side images of four examples of GC and ICM
filaments are displayed in §2, to illustrate their similar morphology,
followed by a discussion of the properties that are common in two populations
of GC and ICM filaments in §3. We then compare the physical parameters of the
two populations of filaments in §4, and discuss their implications for models
of their structure and origin in §5.

##  2 GC versus ICM filaments: four examples

The GC filaments are rich in their morphological variety and details as well
as being abundant in the inner few hundred pc of the Galactic center. ICM
filaments are being revealed in high dynamic images of radio galaxies and
their surrounds. We have selected a subset of GC filaments based on their
similar morphology to ICM filaments. The comparison of limited number of ICM
filaments with a large number of GC filaments offers an opportunity to advance
our understanding of the origin of the filaments and to probe the magneto-
ionized media of the GC and ICM. We examine below the morphology and
properties of filaments associated with four RGs and compare them with seven
GC filaments.

###  2.1 3C40B in Abell 194 vs G359.132-0.296 (the Snake), G359.808+0.130 (the
Flamingo), and G359.717+0.228 (the Eyebrow)

Figure 1 displays a comparison of the filamentary structures associated with
the radio jet of 3C40B and three morphologically similar groups of filaments
identified in the Galactic center (Yusef-Zadeh et al., 2022b). 3C40B is a
radio galaxy at a distance of $\sim 75.5$ Mpc, embedded in a poor cluster
environment Abell 194 (Knowles et al., 2022) which has recently been studied
in detail (Rudnick et al., 2022). Two remarkable filaments emerge from a
deviated region of the jet and run parallel to each other as they bend
together with two different curvatures leading to diffuse ends. The angular
spacing between the filaments is $\sim 30^{\prime\prime}$ ($\sim$11 kpc),
their widths range between 0.85 and 2.8 kpc and both have an extent of $\sim
220$ kpc (Rudnick et al., 2022). The filaments show signs of physical
interaction with a region of the the jet that is highly bent, suggesting a
scenario in which a cloud is crossing the jet and the filaments are then
dragged from the interaction site (Rudnick et al., 2022).

The Snake has a length of $\sim 30^{\prime}$ ($\sim 70$ pc), is one of the
longest in the Galactic center, and is characterized by two kinks and three
different curvatures along its length. Another two examples of GC filaments
are the groups known as the Flamingo and the Eyebrow. The Flamingo shows
multiple filaments with similar curvatures and the Eyebrow consists of two
parallel filaments with the same curvature suggesting that they are associated
with each other.

The morphology of the four filament groups shown in Figure 1 are similar to
each other to the extent that it is difficult to distinguish which one is in a
galaxy cluster based on appearance alone. The systems of filaments in the GC
and 3C40B have similar morphologies, surface brightnesses, length-to-width
ratios and angular spacings. What is remarkable is that the ICM filaments in
3C40B are $\sim 7000$ times longer than nonthermal filaments (aka NTFs) in the
GC with coherent magnetic structure, steepening spectral index, streaming
cosmic-ray particles over such a large distance (Rudnick et al. (2022)) and
yet they show similar physical characteristics.

The ranges of the spectral indices are $\alpha$ = [-0.60,-0.79] for the Snake,
$\alpha$ = [-0.80,-1.05] for the Flamingo, $\alpha$ = [-0.51,-0.76] for the
Eyebrow (Yusef-Zadeh et al., 2022b). These estimates of spectral index are
made along individual filaments. The GC spectral index values are all flatter
than the spectral indices of the filaments in 3C40B which have a mean value of
$\alpha\sim-2$ and steepen from -1.3 to -2.5 with distance from the jet.

The equipartition magnetic field of GC filaments ranges between 20-100 $\mu$G
(for a cosmic ray proton to electron ratio of $p/e=1$) which is much higher
than that of 3C40B with the magnetic field of 2-7$\mu$G. The equipartition
magnetic field of the Snake has a maximum of $\sim 0.1$ mG and decreases away
from the Galactic plane (Yusef-Zadeh et al., 2022a). Furthermore, the
direction of the magnetic field runs along the 3C40B and the Snake filaments.
(The magnetic field directions of the Flamingo and Eyebrow are unknown).

###  2.2 Radio galaxy IC 4296 vs G359.992-0.572 (The Comet Tail)

IC 4296 is a low-luminosity elliptical galaxy, in the cluster A3566 at a
distance of 50 Mpc (Condon et al., 2021). A sensitive radio image of this
galaxy shows a remarkable group of three filaments with angular spacing of
$\sim 20^{\prime\prime}$ (called threads in Condon et al. (2021)) that appear
to converge to a point. The projected lengths and widths of the filaments are
$\sim 50$ and $\sim 2$ kpc, respectively, with surface brightness of $\sim
0.1$ mJy per $\sim 7.6^{\prime\prime}$ beam (Condon et al., 2021). This group
of filaments lies to the south of the western radio jet where the jet shows a
wiggle, thought to be due to Kelvin-Helmholtz (KH) instability (Condon et al.,
2021). These authors suggest that relativistic electrons escape from the jet
where there is KH instability along the radio jets (Condon et al., 2021).
There is also a single isolated filament with different position angle that
arises from the eastern jet. One puzzle is that the spectral index of the
filaments steepens from $\alpha\sim-1$ to -1.5 closer to the injection point
of the radio jet.

The GC grouping of the Comet Tail consists of six linear structures which
resemble a fragmented comet tail. Some of the filaments branch out into
fainter components. The spectral index of individual filaments vary and show
steep spectra in the range [-1,-1.8]. There is a trend that the GC filaments
with steeper spectral indices lie at high latitudes (Yusef-Zadeh et al.,
2022a). The comparison indicates similar morphology and steep spectral indices
in both the Comet Tail and IC 4296 filaments.

###  2.3 The core of A3562 galaxy cluster vs G359.411-0.709 (The Feather)

Figure 3a shows a radio image of the core of the A3562 galaxy cluster at a
distance of $\sim 200$ Mpc (Giacintucci et al., 2022). A linearly polarized
filament with a length of 50 kpc runs perpendicular to the end of the tail of
the radio galaxy J1333-3141. The filament that arises from the region where
the jet terminates is distorted before it splits into two components
(Giacintucci et al., 2022). The splitting occurs 60 kpc away from the radio
tail. The mean spectral index of the filament is $\alpha\sim-1.5$.

For comparison we show in Figure 3b, a striking GC filament, the Feather which
splits into a two-pronged forked filament with the junction at the location of
a compact source G359.416-0.706 with possible stellar counterpart detected at
IR (Yusef-Zadeh et al., 2022b, c). The compact source is suggested to be an
obstacle that splits the bright filament into two parallel, fainter components
with steeper spectral index values. The comparison noted in Figure 3 provides
another example that radio filaments in both classes of Galactic and
extragalactic filaments break up into multiple components, suggestive of a
flow along the filaments that splits into two components as it runs into an
obstacle.

###  2.4 ESO137-006 vs G359.484+0.1122 (The Bent Harp)

Figure 4a shows a remarkable group of filaments linked to the radio lobes of
ESO137-006, a luminous radio galaxy in the Norma galaxy cluster at 70 Mpc from
us (Ramatsoku et al., 2020). The wide-angle tail morphology of this radio
galaxy includes multiple straight and bent filaments extending from the
eastern lobe toward the western lobe. The longest straight filaments act as a
bridge connecting the two radio lobes. In addition, a number of semi-
regularly-spaced filaments are noted to the south of the lobes. The length and
width of the longest filament are 80 and 1.2 kpc, respectively, with an aspect
ratio (length:width) of 67 (Ramatsoku et al., 2020). The mean spectral index
of the filaments, $\alpha\simeq-2$, is similar to the spectral indices of the
lobe from which they emerge. The western half of the galaxy is completely
surrounded by diffuse X-ray emission. The western lobe is also highly
distorted where diffuse and filamentary structures are present. The filaments
bridging the lobes, with similar morphology to the above three examples, arise
from a region that the jet is expanding into the lobe.

G359.484+0.1122 consists of five regularly spaced filaments with the
appearance of a bent harp (Thomas, Pfrommer, & Enßlin, 2020; Yusef-Zadeh et
al., 2022b). A trend was noted between the length of the filaments and their
brightness as well as with the steepening of the spectral index in the range
between [-1.3,-1.0] (Yusef-Zadeh et al., 2022b).

The longest and straightest filament in ESO137-006 (bridging the radio lobes)
appears like many straight GC filaments running perpendicular to the Galactic
plane (Yusef-Zadeh et al., 2022b). The Bent Harp is one group of GC filaments
with a resemblance to the ICM structures in ESO137-006. Unlike the filaments
of ESO137+006, the Bent Harp is not connected to any obvious sources of cosmic
rays.

##  3 General physical properties

Galactic and extra-galactic filaments are produced in totally different
environments, namely, the nucleus of a normal galaxy and the intracluster
medium. Their length scales, energy densities, and other characteristic
parameters differ by orders of magnitude. In spite of these differences, we
argue that similar processes operate in both systems motivated by their
similar morphology as well as similar dimensionless ratios of their physical
parameters. We will examine and compare four examples of GC and RGs in the
next section, following the common physical characteristics, as described
below.

$\bullet$ There is a remarkable similarity in morphological details of the GC
and ICM filaments despite their distances differing by a factor of $\sim
10^{3-4}$. Both classes of filaments are narrow and long with typical aspect
ratios ranging between 10 to 50, assuming that the widths of the filaments are
resolved. The filaments show bending, and wiggles along their lengths and have
similar surface brightness roughly ranging between 0.01 and 0.1 mJy beam-1.

$\bullet$ Both classes of filaments appear either as single isolated filaments
or multiple filaments grouped together. Some of the apparently single
filaments of the GC are found to consist of bundles of filaments in higher
$\sim 1^{\prime\prime}$ resolution VLA observations. Many filaments in
groupings are approximately equally separated from each other and run parallel
to each other, giving a harp-like appearance. Filaments in some groupings
converge to a point or to an extended source. Groups of filaments shift
sideways together, changing direction coherently, implying that they are parts
of the same system of filaments with similar origin.

$\bullet$ The mean spacings between parallel filaments in a given GC grouping
peaks at $\sim 16^{\prime\prime}$ (Yusef-Zadeh et al., 2022b). There is a
limited number of groups of ICM filaments but they show regular spacing
between parallel filaments, as noted in numerous groups of GC filaments. In
the case of 3C40B, the spacing between two parallel filaments is $\sim
30^{\prime\prime}$ (Rudnick et al., 2022). We expect that the spacing to
length ratio of ICM and GC filament to be similar to each other if similar
processes operate in both systems.

$\bullet$ Both systems of GC and ICM filaments trace synchrotron emission from
cosmic-ray electrons. They both show a high percentage of polarization, $\sim
50\%$, revealing ordered magnetic fields (B) running along the filaments. The
orientation of the field is expected from velocity shear at the interface
between the ICM filaments and the ambient gas in the ICM (Condon et al.,
2021).

$\bullet$ The mean magnetic field strengths along the GC filaments range
between $\sim 100$ to 400 $\mu$G depending on the assumed ratio of cosmic-ray
protons to electrons (p/e). The magnetic field of ICM filaments is typically
estimated to be a few $\mu$G, 50 to 100 times smaller. The low magnetic field
of ICM filaments is expected because the filaments have large depth L (see
equation 1), and old if steep spectral index, where flux density
$S_{\nu}\,\propto\,\nu^{\alpha}$, is an indication of age of cosmic-ray
electrons. The weaker magnetic field in the ICM filaments imply long cooling
time scales which scales as $B^{-3/2}$ for observations at a fixed frequency.

$\bullet$ The mean spectral indices of the GC filaments are steeper than
supernova remnants (SNRs, $\alpha\sim-0.62$) with a value of $\alpha\sim-0.8$.
Furthermore, a trend is noted in steepening of the spectral index as a
function of absolute Galactic latitude but not necessarily for typical GC
filaments (Yusef-Zadeh et al., 2022a; LaRosa et al., 2000). With limited
number of ICM filaments, the typical spectral indices appear steeper than GC
filaments, which is not surprising given the vastly longer time scales. In
addition, the GC filaments show a number of filaments, such as the radio Arc
near l$\sim 0.2^{\circ}$, with relatively flat spectral index (Heywood et al.,
2022; Yusef-Zadeh et al., 2022a).

$\bullet$ In some of ICM filaments, radio lobes and tails of radio galaxies
appear to provide injection of cosmic ray particles to illuminate radio
filaments. However, the sources of particle injection in the GC filaments has
not been well established. Recent MeerKAT measurements indicate a number of
compact radio sources located at one or other end of the GC filaments. This
suggests that compact radio sources are potentially responsible for injection
of cosmic ray particles into radio filaments or acting as interfaces to
collimate cosmic-ray driven winds.

$\bullet$ Diffuse thermal X-ray gas and diffuse nonthermal radio emission in
the ISM of the GC and ICM of RGs suggest another similarity in the physical
conditions of these environment. In addition, there is evidence of high
cosmic-ray ionization rate from $H_{3}^{+}$ studies of the Galactic center
region.

##  4 Origin of two populations of filaments

The previous sections described the overall physical properties of a subset of
ICM and GC filaments, as members of two populations of magnetized radio
filaments. There are hundreds of radio filaments found in the Galactic center
at similar distance. A statistical study characterizing the mean properties of
the spectral index, the equipartition magnetic field, the spacing between the
filaments and surface brightness, is based on MeerKAT data at 1.28 GHz (Yusef-
Zadeh et al., 2022a, b). The physical parameters of GC and ICM filaments are
listed in columns 2 and 3, respectively, of Table 1. In this table, the range
and mean values of the length, width, magnetic field strength, surface
brightness and spacing between the filaments are given. The lengths of the GC
filaments follow a power-law distribution (Yusef-Zadeh et al. 2022, in
preparation). At present only a handful of ICM filaments are known, thus there
are not sufficient measurements to do statistics beyond simple means. This is
because ICM filaments require high dynamic range images to uncover faint
magnetized filaments and only limited observations are presently available.
The GC filaments have a wide range of surface brightnesses with many more
fainter filaments. The faint end is cutoff by limited sensitivity. We
recognize that the ICM filaments detected so far quite likely are the most
prominent members of a larger population, with many fainter filaments yet to
be detected.

Table 1: The physical characteristics of the Galactic Center and radio galaxy filaments  Physical Parameters | Galactic Center Filaments | Radio Galaxy Filaments
---|---|---
length (pc) | [4,60] | 103-5
typical width (pc) | 0.5 | few$\times 10^{3}$
typical aspect ratio | [10,100] | [10,70]
magnetic field strength (mG) | [0.1,0.6] | few$\times 10^{-3}$
spectral index ($\alpha$) |  $\sim-0.8$ (mean) [-2,0] | [-2,-1]
surface brightness (mJy beam-1) | [0.01-10] | [0.01,10]
spacing between filaments (pc) |  $\sim$0.7 mean [0.4,1.2] | [10-20]$\times 10^{3}$

###  4.1 Physical parameters of GC and ICM filaments

We compare the measurements and estimates of the physical characteristics of
the two populations of filaments at the GC and in the ICM of galaxy clusters,
as tabulated in Table 1, to explore the implications for the formation of the
magnetic structures and the source of the ultra-relativistic electrons
responsible for the observed radio emission. While the filament length scales
and magnetic field strengths differ by orders of magnitude, it is their
respective dimensionless ratios of time scales and length scales that should
be compared. This allows us to examine whether the underlying physical
mechanisms responsible for the creation of the two populations are similar. In
particular we are interested in the mechanisms responsible for creating
magnetized filaments, and for accelerating the relativistic electrons.

First, the morphological parameters of the two populations of filament are
similar: the aspect ratios are $\sim 30:1$, multiple parallel filaments,
splitting into two pronged fork, bending and spacing between the filaments.
The filaments appear in parallel bundles or converging toward a point in both
populations with a spacing of a few filament widths. These similarities are
hardly surprising given that the motivation behind making this comparison is
that structures show a remarkable resemblance in the two environments, as
shown in Figures 1-4.

Second, the filaments’ magnetic field are estimated to be of order 100 and
$3\,\mu$G in the GC and ICM, respectively. Their corresponding magnetic energy
densities, $B^{2}/8\pi\approx 3\times 10^{6}$ and $3\times 10^{3}$ cm-3 K, are
comparable with or lower than the thermal pressure of the medium in which they
are embedded. These characteristics suggest that the filaments are confined by
the thermal pressure of the external medium. This is supported in the case of
the filaments in Abell 194 and IC 4296 (Rudnick et al. (2022); Condon et al.
(2021)), where the X-ray surface brightness of the intracluster medium is
reduced along the filament, consistent with exclusion of hot gas from the
filament’s interior. In the case of GC filaments, high external cosmic ray
energy density inferred from H${}_{3}^{+}$ measurements contributes to
confinement of the GC filaments (Yusef-Zadeh & Wardle, 2019).

Third, the spectral index values are steeper for ICM filaments than the mean
spectral index of the GC filaments, as given in Table 1. Both systems of
filaments also show steepening of the spectral index along individual
filaments. The relativistic electron populations responsible for the observed
synchrotron emission at 1.2 GHz, have steep spectral indices and show high
percentage of linear polarization based on sensitive observations for both the
GC and ICM filament populations.

For the equipartition field strengths mentioned above, the characteristic
energy of the electrons dominating the emission at 1.2 GHz is $E=(4\pi
m_{e}c\nu/3eB)^{1/2}\,m_{e}c^{2}\approx 0.9$ and 5 GeV, respectively. The
synchrotron power radiated by an ultra-relativistic electron of energy
$E=\gamma m_{e}c^{2}$ is $\dot{E}=2\sigma_{T}c\gamma^{2}B^{2}/4\pi$ where for
simplicity we have assumed a $90^{\circ}$ pitch angle. Then the synchrotron
loss time scale $t_{s}=E/\dot{E}\approx 0.9$ or 200 Myr for the GC or IC
filaments, respectively. As electron acceleration mechanisms tend to produce
electrons with $E^{-2}$ energy spectrum and hence a $\nu^{-0.5}$ synchrotron
spectrum, the steepness of the observed spectra imply that the filament ages
are longer than the synchrotron life times given above. In addition the ages
cannot be too much more than this because the filaments would rapidly become
faint. Thus the synchrotron life times give a rough idea of the age of the
filaments, and impose a requirement on the cosmic-ray electron transport
mechanism which must be able to distribute relativistic electrons along the
observed filament lengths on this time scale.

###  4.2 Transport of cosmic-ray electrons: diffusion versus streaming

The disparate synchrotron loss time scales need to be placed in context. For
source models involving acceleration or injection of relativistic electrons at
a localized site, they impose a requirement to be able to transport the
electrons via diffusion or streaming from the injection site throughout the
observed filament lengths $L\sim 30\,$pc and $100\,$kpc for the GC and ICM
populations, respectively.

If transport is by diffusion, then the diffusion time scale
$L^{2}/\kappa\lesssim t_{s}$ so the diffusion coefficient $\kappa\gtrsim
L^{2}t_{s}\sim 3\times 10^{26}$ and $2\times 10^{31}$ cm2 s-1, respectively,
implying mean free paths $l=3c/\kappa\sim 0.01$ and 570 pc. The required
diffusion is plausible for GC filaments with $L\lesssim 30$ pc but seems
rather large in the case of the ICM filaments (Rudnick et al., 2022) as well
as GC filaments, such as the Snake, with lengths $L\gtrsim 30$ pc. (see Fig.
1).

The alternative to diffusion is that the cosmic ray pressure is sufficient to
allow them to stream along the magnetic field lines, in which case scattering
by self-generated waves limits the streaming to the Alfvén speed, $v_{A}$. The
requirement to populate the entire filament length then implies that
$L/v_{A}\lesssim t_{s}$, i.e. $v_{A}\gtrsim L/t_{s}\sim 30$ and $\sim 600$ km
s-1 for the GC and ICM filaments, respectively. The density of the ionized
thermal plasma present within the filaments is unknown, but we can usefully
estimate the maximum density able to yield the required Alfvén speeds, using
the characteristic field magnetic strengths $B\sim 100$ and 3 $\mu$G for the
two populations. This implies that the thermal electron density $n_{e}\lesssim
3\times 10^{11}$ and $8\times 10^{5}$ cm-3 for the GC and ICM filaments,
respectively. These upper limits clearly exceed any plausible value in each of
the environments by orders of magnitude, so we can conclude another similarity
in parameter space between the two populations that streaming would be able to
transport electrons from their injection site along the entire length of the
observed filaments.

###  4.3 ICM and GC filament models

Models of ICM filaments benefit from the fact that there is an obvious
reservoir of cosmic-ray particles in jets, lobes, mini radio halos (Brunetti &
Jones, 2014), that can be injected into filaments. For example, the filaments
of radio galaxy IC 4296 are thought to originate where helical KH instability
disrupts the flow of the radio jet followed by the relativistic particles
feeding the filaments (Condon et al., 2021). One puzzle is why the spectrum
flattens instead of steepening moving away from the presumed injection point
at the end of the cometary structure’s tail towards the convergence point of
the filaments. One possible explanation for this is that the injected
particles’ pitch angle would tend to increase to conserve the magnetic moment
$p_{\perp}/B$ as they spiral along the converging field lines. This mechanism
relies on the initial pitch angle of the injected particles being anti-
correlated with energy, and pitch-angle scattering to be ineffective in
isotropizing the pitch-angle distribution.

In the case of the A3562 galaxy cluster, it has been suggested that the tail
of the radio galaxy is interacting with a tangential wind blown in the
direction along the filaments (Giacintucci et al., 2022). An alternative to
this picture is a scenario in which the cosmic-ray electrons injected from the
tail, feed cluster magnetic field lines (Giacintucci et al., 2022). In another
model, cosmic-ray electrons associated with two parallel filaments of 3C40B
are injected from a region where the jet is highly distorted (Rudnick et al.,
2022). In particular, an impact between a moving dense cloud and the jet is
expected to distort the shape of the jet as the magnetic field lines are
dragged away from the site of the interacting jet (Rudnick et al., 2022).
Another remarkable object is source C in A2256 with its extremely narrow and
long tail (owen14). This unusual source is considered to be a one-sided radio
jet. However, it has a very similar morphology to GC radio filaments. Lastly,
the collimated synchrotron filaments in ESO 137-006 could be due to the
interaction of the magnetic fields of the radio lobes with the magneto-ionic
plasma of the intracluster medium (Ramatsoku et al., 2020).

###  4.4 Scenarios for the origin of the ICM and GC filaments

In GC and ICM populations, there are two broad scenarios that could
potentially explain how magnetic filaments are formed. Both scenarios posit
amplification and organization of an initially weak magnetic field in the
surrounding medium. In the first scenario, the filaments are formed by an
initially weak field that is amplified and structured by subsonic turbulent
motions in a weakly magnetized medium, a scenario that has been proposed both
for GC (Boldyrev & Yusef-Zadeh (2006)) and the ICM environments (Porter,
Jones, & Ryu, 2015; Vazza et al., 2018; Rudnick et al., 2022). The challenge
in this picture is that simulations show the overall orientation of most
filaments is random (Porter, Jones, & Ryu, 2015; Vazza et al., 2018). The GC
filaments are mainly directed in the direction away from the Galactic plane.
It is possible that alignment of GC filaments is generated by cosmic-ray
driven wind running away from the Galactic plane, thus creating a velocity
shear in the turbulent medium that aligns the filaments.

In the second scenario the filaments are created by subsonic flow of a weakly
magnetized medium past an obstacle. The passing flow drapes the obstacle with
field lines that are then compressed and stretched back into comet-like tail,
creating a magnetic filament with magnetic pressure of order the thermal ram
pressure in the medium. The obstacle could be any relatively dense ionized
clump that creates a bow shock in the flow such as a stellar wind bubble, a
red giant, an expanding HII region, wind tubes associated with red giants in
AGN relativistic bubbles (Chugai, Churazov, & Sunyaev, 2011), or a disturbed
region along jets or lobes showing size scales similar to the separation
between parallel filaments. In this picture, the magnetized filaments are
produced at the interaction sites of a large-scale driven wind outflowing away
from the Galactic plane and embedded stellar wind bubbles (Rosner & Bodo,
1996; Shore & LaRosa, 1999; Yusef-Zadeh & Wardle, 2019). The large-scale Radio
Bubble filled with X-ray gas in the Galactic center is considered to be
arising from cosmic-ray driven outflow (Yusef-Zadeh & Wardle, 2019). In this
picture, the obstacle sets the length scale of the separation between the
filaments. This is analogous to the interaction of fast-moving mass-losing
stars whose winds interact with the ISM and create cometary tails (Martin et
al. 2007). In addition to cometary model, numerous other models have been
proposed in the past to explain the origin of the GC filaments and the sources
in which cosmic-ray electrons are accelerated. Unlike ICM filaments, it is not
clear what feeds cosmic-ray electrons into filaments (Nicholls & Le Strange,
1995; Rosner & Bodo, 1996; Shore & LaRosa, 1999; Bicknell & Li, 2001;
Dahlburg, et al., 2002; Yusef-Zadeh, 2003; Boldyrev & Yusef-Zadeh, 2006;
Ferrière, 2009; Banda-Barragán, et al., 2016; Yusef-Zadeh & Wardle, 2019;
Thomas, Pfrommer, & Enßlin, 2020; Sofue, 2020; Coughlin, Nixon, & Ginsburg,
2021).

In the context of the ICM filaments, the relative motion with the external
medium is created by the orbital motion of the obstacle or a head-tail radio
galaxy through the ICM around the center of the cluster. Additional
contribution for slow-moving galaxies close to the center of the cluster could
be moving through denser ICM (Toniazzo & Schindler, 2001; Douglass et al.,
2008). Other sources of external pressure could be due ”weather” motion within
the ICM (Brunetti & Jones, 2014). On the other hand, in the Galactic center
ISM, there is contribution due to a large-scale nuclear wind generated by high
cosmic-ray pressure (Yusef-Zadeh & Wardle, 2019).

One interesting aspect of this picture when applied to ICM filaments is that
the filaments are linked to jets or lobes of radio galaxies. The locations
where the filaments cross radio galaxies are highly distorted and cosmic-ray
particles are likely to escape from the jet (Condon et al., 2021). This
supports the picture in which there is a concentration of cosmic-ray particles
in the jet or the lobe that injects relativistic plasma into the ICM
filaments. As described in §4.3, the origin of filamentation is interaction
with an obstacle and splitting of the filament into multiple filaments.
Another possibility for the origin of the filamentation is synchrotron cooling
instability formed as a result of the interaction of cosmic-ray generated by
the winds (e.g., Simon & Axford, 1967). The mean spacing and the mean width of
the filaments are expected to be similar to each other (Yusef-Zadeh et al.,
2022a). In the case of ICM filaments, the spacing is estimated to be the
product of the cooling time scale and the Alfvén speed which is of the order
of kpc.

Assuming that the GC and ICM filaments are formed with the same mechanism in
the context of an origin involving an obstacle collimating the magnetized
wind, the lack of compact radio sources associated with some of the GC
filaments is likely to be related to the lifetime of the obstacle. In
particular, if an obstacle is an expanding HII region or a planetary nebula,
then it may have already dissipated as simulations of cloud-wind interactions
indicate (Banda-Barragán, et al., 2016). Another possibility is that the
source of injection has a very steep spectrum and would be visible only at
very low frequencies.

The above cometary and turbulent models are technically different than a
scenario in which there is a strong pre-existing organized field dominating
the region, with only some field lines being lit up by an injected population
of cosmic ray electrons. This model was one of the first models proposed for
the GC which was originally thought to be threaded by a milligauss poloidal
field (Morris, 2007). However, the diffuse radio emission from inner few
hundred pc of the galaxy suggests a much weaker global field, $\sim
10-20\,\mu$G (Yusef-Zadeh et al., 2013). This scenario is also inapplicable to
the ICM environment in which the field strengths are much less than several
$\mu$G inferred for the magnetized filaments.

##  5 Summary

We compared the two populations of magnetized filaments in the Galactic Centre
and ICM. They differ vastly in their physical properties, such as lengths,
widths and magnetic field strengths. However, we argued that nevertheless,
they are analogous to each other as might be anticipated based on their
similar morphologies.

In both cases, the filaments are in rough pressure equilibrium with their
surroundings, but are more strongly magnetized, consistent with scenarios in
which they are formed by dynamical processes at work in their surroundings. In
one scenario, which has previously been discussed in both the GC and ICM
contexts, the filaments arise through the stretching and collection of field
lines by turbulence in weakly magnetized medium. An alternative scenario that
has been proposed for the GC filaments is the collection and draping of field
lines by a moving stellar wind source or similar obstacle with respect to the
medium. It is not clear that this mechanism can create the large aspect ratios
of the ICM filaments, because any orbital motion is mildly transonic at best
and so is unlikely to lead to a tight magnetized tail or lead to sufficient
particle acceleration to energize it.

The source of the relativistic particles still remains to be established. In
the GC, a few filaments have associated compact radio sources, and in the ICM
the filaments are associated with radio jets and lobes, but the exact nature
of any connection in both cases remains unclear. Both filament populations
have steep synchrotron spectra indicating an aged cosmic-ray electron
population, so it is possible that any injection event happened long ago and
the link to the injecting source is no longer apparent.

The striking similarities between these enigmatic populations, despite their
very different environments, suggests the possibility that one or the other
may be amenable to observational probes that shed light in the physical
processes at work in both populations. Some examples include: (i) more
sensitive MeerKAT observations of ICM sources to look for fainter as well as
shorter (young and bright) filaments, (ii) high resolution $1^{\prime\prime}$
VLA study of ICM and GC filaments, (iii) X-ray studies searching for
depressions in X-ray surface brightness overlapping ICM and GC filaments.

## Acknowledgments

Work by R.G.A. was supported by NASA under award number 80GSFC21M0002. We are
grateful to W. Cotton, S. Giacintucci, M. Ramatsoku and L. Rudnick for
providing their FITS images of their published data. We thank L. Rudnick and
the referee for useful comments. The MeerKAT telescope is operated by the
South African Radio Astronomy Observatory, which is a facility of the National
Research Foundation, an agency of the Department of Science and Innovation.
The National Radio Astronomy Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc.

## References

  * Arendt et al. (2019) Arendt R. G., Staguhn J., Dwek E., Morris M. R., Yusef-Zadeh F., Benford D. J., Kovács A., et al., 2019, ApJ, 885, 71. doi:10.3847/153
  * Bally & Yusef-Zadeh (1989) Bally J., Yusef-Zadeh F., 1989, ApJ, 336, 173\. doi:10.1086/167003
  * Banda-Barragán, et al. (2016) Banda-Barragán W. E., Parkin E. R., Federrath C., Croc ker R. M., Bicknell G. V., 2016, MNRAS, 455, 1309
  * Bicknell & Li (2001) Bicknell G. V., Li J., 2001, ApJ, 548, L69
  * Boldyrev & Yusef-Zadeh (2006) Boldyrev S., Yusef-Zadeh F., 2006, ApJL, 637, L101. doi:10.1086/500411
  * Brienza et al. (2022) Brienza M., Lovisari L., Rajpurohit K., Bonafede A., Gastaldello F., Murgia M., Vazza F., et al., 2022, A&A, 661, A92. doi:10.1051/0004-6361/202142579
  * Brunetti & Jones (2014) Brunetti G., Jones T. W., 2014, IJMPD, 23, 1430007-98. doi:10.1142/S0218271814300079
  * Chuard et al. (2018) Chuard D., Terrier R., Goldwurm A., Clavel M., Soldi S., Morris M. R., Ponti G., et al., 2018, A&A, 610, A34. doi:10.1051/0004-6361/20173
  * Chugai, Churazov, & Sunyaev (2011) Chugai N. N., Churazov E. M., Sunyaev R. A., 2011, MNRAS, 414, 879. doi:10.1111/j.1365-2966.2011.18274.x
  * Comisso (2020) Comisso L., Sobacchi E., Sironi L., 2020, ApJL, 895, L40
  * Condon et al. (2021) Condon J. J., Cotton W. D., White S. V., Legodi S., Goedhart S., McAlpine K., Ratcliffe S. M., et al., 2021, ApJ, 917, 18. doi:10.3847/1538-4357/ac0880
  * Coughlin, Nixon, & Ginsburg (2021) Coughlin E. R., Nixon C. J., Ginsburg A., 2021, MNRAS, 501, 1868. doi:10.1093/mnras/staa3771
  * Dahlburg, et al. (2002) Dahlburg R. B., Einaudi G., LaRosa T. N., Shore S. N., 2002, ApJ, 568, 220
  * Douglass et al. (2008) Douglass E. M., Blanton E. L., Clarke T. E., Sarazin C. L., Wise M., 2008, ApJ, 673, 763. doi:10.1086/523886
  * Drake, Swisdak, & Fermo (2013) Drake J. F., Swisdak M., Fermo R., 2013, ApJL, 763, L5. doi:10.1088/2041-8205/763/1/L5 Shay M. A., Drake J. F., Karimabadi H., 2013, AGUFMSM13 Zweibel E. G., 713, 592
  * Fanaroff et al. (2021) Fanaroff B., Lal D. V., Venturi T., Smirnov O. M., Bondi M., Thorat K., Bester L. H., et al., 2021, MNRAS, 505, 6003. doi:10.1093/mnras/stab1540
  * Ferrière (2009) Ferrière K., 2009, A&A, 505, 1183 T., 1999, ApJ, 510, 251
  * Goto, et al. (2014) Goto M., Geballe T. R., Indriolo N., Yusef-Zadeh F., Usuda T., Henning T., Oka T., 2014, ApJ, 786, 96
  * Giacintucci et al. (2022) Giacintucci S., Venturi T., Markevitch M., Bourdin H., Mazzotta P., Merluzzi P., Dallacasa D., et al., 2022, ApJ, 934, 49. doi:10.3847/1538-4357/ac7805
  * Gray, et al. (1991) Gray A. D., Cram L. E., Ekers R. D., Goss W. M., 1991, Natur, 353, 237
  * Green (1988) Green D. A., 1988, Ap&SS, 148, 3. doi:10.1007/BF00646462
  * Guo et al. (2020) Guo F., Liu Y.-H., Li X., Li H., Daughton W., Kilian P., 2020, PhPl, 27, 080501. doi:10.1063/5.0012094 1981, A&A, 103, 197
  * Haynes, et al. (1992) Haynes R. F., Stewart R. T., Gray A. D., Reich W., Reich P., Mebold U., 1992, A&A, 264, 500 Ginsburg A., Sormani M. C., Walker
  * Hardcastle et al. (2019) Hardcastle M. J., Croston J. H., Shimwell T. W., Tasse C., Gürkan G., Morganti R., Murgia M., et al., 2019, MNRAS, 488, 3416. doi:10.1093/mnras/stz1910
  * Heywood, et al. (2019) Heywood I. et al. 2019, Nature, 10.1038 DOI number 10.1038/s41586-019-1532-5
  * Heywood et al. (2022) Heywood I., Rammala I., Camilo F., Cotton W. D., Yusef-Zadeh F., Abbott T. D., Adam R. M., et al., 2022, ApJ, 925, 165. doi:10.3847/1538
  * Knowles et al. (2022) Knowles K., Cotton W. D., Rudnick L., Camilo F., Goedhart S., Deane R., Ramatsoku M., et al., 2022, A&A, 657, A56. doi:10.1051/0004-6361/202141488
  * Lang, Morris & Echevarria (1999) Lang C. C., Morris M., Echevarria L., 1999, ApJ, 526, 727
  * Lane et al. (2002) Lane W. M., Kassim N. E., Ensslin T. A., Harris D. E., Perley R. A., 2002, AJ, 123, 2985. doi:10.1086/340359
  * LaRosa, Lazio & Kassim (2001) LaRosa T. N., Lazio T. J. W., Kassim N. E., 2001, ApJ, 563, 163
  * LaRosa, et al. (2004) LaRosa T. N., Nord M. E., Lazio T. J. W., Kassim N. E., 2004, ApJ, 607, 302
  * LaRosa et al. (2000) LaRosa T. N., Kassim N. E., Lazio T. J. W., Hyman S. D., 2000, AJ, 119, 207. doi:10.1086/301168
  * Law, Yusef-Zadeh & Cotton (2008) Law C. J., Yusef-Zadeh F., Cotton W. D., 2008, ApJS, 177, 515
  * Liszt (1985) Liszt H. S., 1985, ApJ, 293, L65
  * Morris (2007) Morris, M. 2007, astro-ph/0701050
  * Nord et al. (2004) Nord M. E., Lazio T. J. W., Kassim N. E., Hyman S. D., LaRosa T. N., Brogan C. L., Duric N., 2004, AJ, 128, 1646
  * Nicholls & Le Strange (1995) Nicholls J., Le Strange E. T., 1995, ApJ, 443, 638
  * Oka, et al. (2005) Oka T., Geballe T. R., Goto M., Usuda T., McCall B. J., 2005, ApJ, 632, 882
  * Oka et al. (2019) Oka T., Geballe T. R., Goto M., Usuda T., Benjamin, McCall J., Indriolo N., 2019, ApJ, 883, 54. doi:10.3847/1538-4357/ab3647
  * Owen et al. (2014) Owen F. N., Rudnick L., Eilek J., Rau U., Bhatnagar S., Kogan L., 2014, ApJ, 794, 24. doi:10.1088/0004-637X/794/1/24 H., Mao S. A., 2019, ApJ, 884, 170.
  * Ponti et al. (2019) Ponti, G., Hofmann, F., Churazov, E., et al. 2019, Nature, 567, 347
  * Porter, Jones, & Ryu (2015) Porter D. H., Jones T. W., Ryu D., 2015, ApJ, 810, 93. doi:10.1088/0004-637X/810/2/93
  * Pound & Yusef-Zadeh (2018) Pound M. W., Yusef-Zadeh F., 2018, MNRAS, 473, 2899. doi:10.1093/mnras/stx2490
  * Royster et al. (2011) Royster M. J., Hewitt J. W., Roberts D. A., Yusef-Zadeh F., 2011, ASPC, 439, 31
  * Rosner & Bodo (1996) Rosner R., Bodo G., 1996, ApJ, 470, L49
  * Ramatsoku et al. (2020) Ramatsoku M., Murgia M., Vacca V., Serra P., Makhathini S., Govoni F., Smirnov O., et al., 2020, A&A, 636, L1. doi:10.1051/0004-6361/202037800
  * Rudnick et al. (2022) Rudnick L., Bruggen M., Brunetti G., Cotton W., Forman W., Jones T. W., Nolting C., et al., 2022, arXiv, arXiv:2206.14319
  * Shimwell et al. (2016) Shimwell T. W., Luckin J., Brüggen M., Brunetti G., Intema H. T., Owers M. S., Röttgering H. J. A., et al., 2016, MNRAS, 459, 277. doi:10.1093/mnras/stw661
  * Shore & LaRosa (1999) Shore S. N., LaRosa T. N., 1999, ApJ, 521, 587
  * Simon & Axford (1967) Simon M., Axford W. I., 1967, ApJ, 150, 105. doi:10.1086/149316
  * Staguhn et al. (2019) Staguhn J., Arendt R. G., Dwek E., Morris M. R., Yusef-Z adeh F., Benford D. J., Kovács A., et al., 2019, ApJ, 885, 72. doi:10.3847/1538-4357/ab451b
  * Staguhn, et al. (1998) Staguhn J., Stutzki J., Uchida K. I., Yusef-Zadeh F., 1998, A&A, 336, 290
  * Sofue (2020) Sofue Y., 2020, PASJ, 72, L4. doi:10.1093/pasj/psaa011
  * Thomas, Pfrommer, & Enßlin (2020) Thomas T., Pfrommer C., Enßlin T., 2020, ApJL, 890, L18. doi:10.3847/2041-8213/ab7237
  * Toniazzo & Schindler (2001) Toniazzo T., Schindler S., 2001, MNRAS, 325, 509. doi:10.1046/j.1365-8711.2001.04385.x T., 1995, PASJ, 47, 829
  * van Weeren et al. (2021) van Weeren R. J., Shimwell T. W., Botteon A., Brunetti G., Brüggen M., Boxelaar J. M., Cassano R., et al., 2021, A&A, 651, A115. doi:10.1051/0004-6361/20203982
  * Vazza et al. (2018) Vazza F., Brunetti G., Brüggen M., Bonafede A., 2018, MNRAS, 474, 1672. doi:10.1093/mnras/stx2830
  * Yusef-Zadeh et al. (2021) Yusef-Zadeh F., Wardle M., Heinke C., Heywood I., Arendt R., Royster M., Cotton W., et al., 2021, MNRAS, 500, 3142. doi:10.1093/mnras/staa3257
  * Yusef-Zadeh (2003) Yusef-Zadeh, F. 2003, ApJ, 598, 325.
  * Yusef-Zadeh & Wardle (2019) Yusef-Zadeh F., Wardle M., 2019, MNRAS, 490, L1
  * Yusef-Zadeh, Morris & Chance (1984) Yusef-Zadeh F., Morris M., Chance D., 1984, Natur, 310, 557
  * Yusef-Zadeh et al. (2022a) Yusef-Zadeh F., Arendt R. G., Wardle M., Heywood I., Cotton W., Camilo F., 2022, ApJL, 925, L18. doi:10.3847/2041-8213/ac4802
  * Yusef-Zadeh et al. (2022b) Yusef-Zadeh F., Arendt R. G., Wardle M., Boldyrev S., Heywood I., Cotton W., Camilo F., 2022, MNRAS, 515, 3059. doi:10.1093/mnras/stac1696
  * Yusef-Zadeh et al. (2022c) Yusef-Zadeh F., Arendt R. G., Wardle M., Heywood I., Cotton W., 2022, arXiv, arXiv:2208.11589 Parastaran P., 1997, ApJ, 475, L119
  * Yusef-Zadeh et al. (2004) Yusef-Zadeh, F., Hewitt, J. W., & Cotton, W. 2004, ApJS, 155, 421.
  * Yusef-Zadeh et al. (2013) Yusef-Zadeh F., Hewitt J. W., Wardle M., Tatischeff V., Roberts D. A., Cotton W., Uchiyama H., et al., 2013, ApJ, 762, 33. doi:10.1088/0004-637X/762/1/33
  * Yusef-Zadeh, et al. (2007) Yusef-Zadeh F., Muno M., Wardle M., Lis D. C., 2007, ApJ, 656, 847

Figure 1:  Left A grayscale MeerKAT image of the radio galaxy 3C40B and its
long filaments (Rudnick et al., 2022). The filaments are connected to a
portion of the jet shown at the top of the panel. Right The three panels to
the right display three groups of GC radio filaments based on 4′′ resolution
MeerKAT mosaic of the Galactic center at 1.28 GHz (Heywood et al., 2022),
discussed in detail in (Yusef-Zadeh et al., 2022b). The intensity I (Jy
beam-1) scales are: 3C40B $-4<\log(I)<-1.4$ only for the fainter emission;
Snake $-4<\log(I+3\times 10^{-4})<-2.33$; Flamingo $(-4.3<\log(I+3\times
10^{-4})<-2.3$; Eyebrow $(-4.3<\log(I+5\times 10^{-4})<-3.$ Figure 2:  Left A
grayscale MeerKAT image of the core and its radio jet (top of the panel) and
the filaments (center of the panel) of IC 4296, an elliptical radio galaxy,
are displayed (Condon et al., 2021). Because of the strong surface brightness
of the core and the jet compared to the filaments, different color tables were
applied. Right A MeerKAT filtered image of the Comet tail as discussed in
(Yusef-Zadeh et al., 2022b). The intensity I (Jy beam-1) scales are: IC4296
$-4<\log(I)<-1.2(I>3.6\times 10^{-4})$; the Comet tail: $-5<\log(I+3\times
10^{-4})<-4$.  Figure 3:  Top (a) A MeerKAT image of the galaxy cluster A 3562
(Giacintucci et al., 2022) convolved with a Gaussian kernel with FWHM $\sim
10.3^{\prime\prime}$. Bottom (b) A portion of the Feather filament where the
filament splits into two fainter components at a compact source (Yusef-Zadeh
et al., 2022b). The intensity I (Jy beam-1) scales are listed as: Feather:
$-4.7<\log(I+3\times 10^{-4})<-4.3$; A3562: $-4<\log(I)<-2(I>7\times
10^{-5})$.  Figure 4:  Top A MeerKAT image of ESO137-006 galaxy in Norma
galaxy cluster (Ramatsoku et al., 2020). The beam FWHM is $\sim
10^{\prime\prime}$. $7.62^{\prime\prime}\times 7.29^{\prime\prime}$. Bottom A
portion of the Bent Harp filament showing multiple parallel filaments (Yusef-
Zadeh et al., 2022b). The intensity I (Jy beam-1) scale for the Bent Harp is:
$-3.8<\log(I+3\times 10^{-4})<-3.$


# Smooth Time-Optimal Trajectories for Steering Drones

Srinath Tankasala1, Can Pehlivanturk1, Efstathios Bakolas1, Mitch Pryor1

###### Abstract

In this paper, we address a minimum-time steering problem for a drone modeled
as point mass with bounded acceleration, across a set of desired waypoints in
the presence of gravity. We first provide a method to solve for the minimum-
time control input that will steer the point mass between two waypoints based
on a continuous-time problem formulation which we address by using
Pontryagin’s Minimum Principle. Subsequently, we solve for the time-optimal
trajectory across the given set of waypoints by discretizing in the time
domain and formulating the minimum-time problem as a nonlinear program (NLP).
The velocities at each waypoint obtained from solving the NLP in the
discretized domain are then used as boundary conditions to extend our two-
point solution across those multiple waypoints. We apply this planning
methodology to execute a surveying task that minimizes the time taken to
completely explore a target area or volume. Numerical simulations and
theoretical analyses of this new planning methodology are presented. The
results from our approach are also compared to traditional polynomial
trajectories like minimum snap planning.

##  I Introduction

The problem of steering a particle using an optimal controller with bounded
inputs has received great attention in the literature, [1, 2, 3], for example.
Several variants of the problem, such as steering a particle in the presence
of a flow field have also been well-studied, [4, 5]. However, previous works
don’t address the case of continuous-time optimal steering of the particle in
the presence of gravity. More importantly, most steering algorithms only
consider the case of steering a particle between two points. There is no
framework to extend the two-point solution to find a continuous-time optimal
input for steering the particle across an arbitrary number of waypoints. In
this work, we solve a minimum time optimal problem between two waypoints while
taking gravity into account. We then provide a way to extend this solution to
include multiple waypoints using a discrete-time problem formulation. By using
this extension of the two point solution, we address the problem of steering a
particle through a set of prescribed waypoints in minimum time. ††1 Cockrell
School of Engineering, The University of Texas at Austin
301 E. Dean Keeton St. C2100 Austin, Texas 78712-2100 It is assumed that the
magnitude of the control input of the particle is upper bounded by an a priori
fixed constant. Hence the thrust input value set for the particle is a
‘hypersphere’. The problem considered here belongs to a class of minimum-time
problems for linear systems where the control input is constrained to a convex
set [6]. When the input set is constrained to a hypersphere, one can
characterize the time optimal control input law. In our case, there always
exists a time-optimal control input that attains its values in the boundary of
the input value set. This is commonly known as the boundary control strategy,
[7], or bang-bang control in the special case when the input set is a
hypercube. Several approaches for generating time-optimal trajectories for
UAVs are available in the literature: Lai et al [8] uses a discretization in
the time domain to formulate the quadrotor steering problem as an non-linear
program (NLP) that is solved using a genetic algorithm. The obtained solution
is time optimal in the discretized input space, but that does not guarantee
that as the time resolution of the discretization increases that one obtains
the optimal solution to the continuous-time problem.

In Loock et al, [9], the trajectory is parameterized based on a set of
reference geometric curves and a non-linear optimal control problem is
formulated to calculate the minimum-time path. In [10] the authors simplify
the quadrotor dynamics to a 2D model and apply Pontryagin’s Minimum Principle
(PMP) to solve for the minimum-time solution. They assume a hypercube domain
and use a bang-bang control strategy. The bang-bang strategy provides a
solution where the system uses maximum actuator input along the boundary of
the admissible control envelope. However, most actuators cannot
instantaneously switch inputs and thus bang-bang solutions are difficult to
track. Hence, a continuous time optimal input would make the generated
trajectory smoother and easier to track.

The bang-bang control strategy in particular has found widespread application
in the field of UAVs as a means to compute minimum-time trajectories [10, 11].
There are many applications where steering the drone through a desired set of
waypoints in minimum time is essential, such as inspection surveys, drone
racing and package delivery, [12, 13]. Quadrotors are however underactuated
systems where some of the rotation and translation degrees of freedom are
coupled. This means that maximum actuator input is not always used for
accelerating the vehicle and some effort is spent on orienting the vehicle
thrust (z-axis) correctly. Nonetheless, the minimum-time solution based on the
point mass simplification of the drone dynamics serves as a starting point/low
fidelity result in many planning strategies [14, 15]. It allows us to evaluate
how much of the theoretical speed potential does a certain control strategy
utilize and identify maneuvers that can improve it. Hence the problem of time-
optimal steering of a point mass with bounded thrust is now of great interest
in the UAV community. To this end, we demonstrate our proposed trajectory
generation strategy by applying it to solve a minimum-time surveying task. The
trajectory generation method can be applied to other tasks as well such as
drone racing, etc.

A common goal of a surveying task is to reconstruct a representation of the 3D
scene (point cloud, textured map, etc.). To complete a survey of a given
Region of Interest (RoI), the agent needs to pass through a known set of
waypoints generated from a global planner. Here it is assumed that the drone
does not need to halt at any of the intermediate waypoints. Determining the
set of waypoints for complete coverage of the RoI depends on the information
available about the terrain. Coverage path planning is well studied in the
literature. Cabreira et al, [16] presents a comprehensive study on the state
of the art in coverage path planning methods. Standard coverage planning
techniques apply cellular decomposition to divide the survey area into evenly
sized cells (rectangles). Given the flight pattern, camera field-of-view
(FoV), required ground resolution, and the desired image overlap, one can
compute the cellular decomposition of the RoI (‘A’) and corresponding
waypoints that need to be visited for ensuring 100% coverage of the desired
area, as shown in Fig 1.

The scope of the surveying task presented here is limited to offline planning,
i.e. when information about the terrain layout (flat, hilly, urban, etc.) is
available a priori. The global planner is also responsible for calculating the
optimal touring path that connects these waypoints, i.e. the order in which
they will be visited. Depending on the user’s objective (minimum time, energy,
etc.) an appropriate touring path can be determined. An example of the
cellular decomposition for our surveying task is shown in Fig 1. The agent is
required to pass through the center of the overlapping rectangles in a known
order.

In this paper, we make three main contributions. Firstly, we present a
systematic method to compute the time optimal trajectory for steering a drone
of point mass from an initial position and velocity to a desired endpoint and
final velocity in minimum time in the presence of gravity and bounded thrust.
We do this by reducing the time-optimal control problem to a set of algebraic
eqs that can be solved to obtain the minimum-time trajectory. Secondly, we
provide a numerical method to address the minimum-time surveying problem by
discretizing the problem in time, similar to [8]. The discretized time domain
formulation is framed as a non-linear program (NLP) that can be solved using
standard optimization solvers. Thirdly, we lay out a framework to extend the
solution from our synthesized time-optimal steering problem to a generalized
waypoint navigation problem. We do this by using the waypoint velocities
obtained from the NLP solution as a boundary condition for our trajectory
generator. By applying this method, we achieve reduced survey times and
smoother inputs making it easier to track on a real drone. The trajectories
generated were faster than traditional planners, like minimum snap planning
[17], completing the survey in $\approx$ 20% lesser time.

Figure 1: Representative pic of the drone and discretization of of survey area
with defined image overlap (yellow rectangles) for stitching

In Section II, we present the dynamic model and solution to the time-optimal
steering problem in the presence of gravity. In Section III, we discuss the
minimum-time surveying problem as an NLP and extend our time-optimal solution
to multiple waypoints. In Section IV, we perform simulations to demonstrate
the efficiency of the trajectory by comparing the NLP solution and the minimum
time solution obtained by PMP. We also show comparisons with traditional
minimum snap solutions. Section V summarizes our conclusions and discusses
avenues for future extension in the planning algorithm.

##  II Problem Formulation

###  II-A Analysis of the optimal steering problem under gravity

We consider the dynamics of the drone as a double-integral plant traveling in
3D space in a constant gravitational field, g. The magnitude of control input
(acceleration) for steering the drone is upper bounded, $\bar{u}$. We assume
the control input has no lower bound, i.e. zero control input is allowed. To
simplify, we take the initial time, $t_{0}=0,\ r(0)=r_{0}$, and $v(0)=v_{0}$.
Then we can assume the motion of the drone is described by the following set
of eqs:

$\begin{split}&\dot{{r}}(t)={v(t)},\quad{\dot{v}(t)}={u(t)+g},\\\
&{\left\lVert u(t)\right\rVert}\leq\bar{{u}}\\\ \end{split}$ (1)

where ${r(t)}\in\mathcal{R}^{3}({r_{0}}\in\mathcal{R}^{3})$ and
${v(t)}\in\mathcal{R}^{3}({v}_{0}\in\mathcal{R}^{3})$ are the position and
velocities of the drone respectively at time ${t}$ (time ${t=0}$), ${u(t)}$ is
the control input at time ${t}$ and $\bar{{u}}$ and ${g}$ are non-negative
constants.
Definition 1 We refer to the composite state of the drone as $z$, which is a
concatenation of the position and velocity vectors,
${z(t)}:=\mathrm{col}({r(t)},{v(t)})$. The initial state of the drone is given
by ${z(0)=z_{0}}=\mathrm{col}({r_{0}},{v_{0}})$.

We rewrite the system in (1) according to the convention in Athans & Falb,
[18] as follows:

$\begin{split}&{\dot{z}(t)}={f(z)+Bu(t)}\\\ &{\left\lVert
u(t)\right\rVert\leq\bar{u}}\\\ \text{where, }&{f(z)}=\begin{bmatrix}{v}\\\
{g}\end{bmatrix};{B=\begin{bmatrix}{0_{3\times 3}}\\\ {I_{3\times
3}}\end{bmatrix}};\\\ &\text{and, }{z(0)=z_{0}}\end{split}$ (2)

We wish to steer the drone from $z_{0}$ at time $t=0$ to the final state
${z_{f}=(r(t_{f}),v(t_{f}))=(r_{f},v_{f})}$ at time ${t=t_{f}}$ for a given
input ${u(t)}$.
Definition 2 Given a initial state $z_{0}$ and final state $z_{f}$ of the
drone, the set of all admissible control input functions that steers the drone
from start is given by $\mathcal{U}$. In other words:
$\mathcal{U}$ is the set of all piecewise continuous functions of time,
$u(t)$, s.t. $||{u(t)}||\leq\bar{{u}}\ \forall\ t\in[0,t_{f}]$ and
${u(t):z_{0}\rightarrow z_{f}}$.

For the system given in eq (2), the Hamiltonian for the time optimal steering
problem can be written as:

$\displaystyle\mathcal{H}{(t,z,u,\lambda)=1+\lambda^{T}(f(v)+Bu(t))}$ (3)

where
$\mathcal{H}{(\cdot):[0,\infty)}\times\mathcal{R}^{6}\times\mathcal{R}^{3}\times\mathcal{R}^{6}\mapsto\mathcal{R}$,
and ${\lambda:[0,t_{f}]}\mapsto\mathcal{R}^{6}$ is the co-state vector. For
simplification purposes, we rewrite ${\lambda(\cdot)}$ as:

$\displaystyle{\lambda=\mathrm{col}(p,q)};\ \ \text{where}\ \
{p,q:[0,t_{f}]}\mapsto\mathcal{R}^{3}$

Thus the Hamiltonian for the system in eq (2) can be rewritten as:

$\displaystyle\mathcal{H}={1+p^{T}v+q^{T}(u+g)}.$ (4)

Next, we apply Pontryagin’s minimum principle (PMP), [18], to characterize the
time-optimal steering input. Let
${z^{*}(\cdot):[0,t_{f}^{*}]}\mapsto\mathcal{R}^{6}$, where
${z^{*}(\cdot)}=\mathrm{col}({r^{*}(\cdot),v^{*}(\cdot))}$ denote the minimum-
time trajectory generated by applying the time optimal control input
${u^{*}(\cdot)}$ over the interval ${[0,t_{f}^{*}]}$. Then, the existence of a
continuous co-state function
${\lambda^{*}(\cdot):[0,t_{f}^{*}]}\mapsto\mathcal{R}^{6}$ corresponding to
the optimal trajectory ${z^{*}(\cdot)}$ is guaranteed and it neccessarily
satisfies the following system of differential eqs:

$\displaystyle\begin{split}&{\dot{\lambda}^{*}}=-\left(\nabla_{z}\mathcal{H}\right)^{T}=\begin{bmatrix}-\left(\frac{{\partial}\mathcal{H}}{{\partial
x}}\right)^{T}\\\ \\\ -\left(\frac{{\partial}\mathcal{H}}{{\partial
v}}\right)^{T}\end{bmatrix}\\\ \end{split}$ (5)

Let ${\lambda^{*}}=\mathrm{col}{(p^{*},q^{*})}$, then from eq (5), it follows
that:

$\displaystyle\begin{bmatrix}{\dot{p}^{*}}\\\
{\dot{q}^{*}}\end{bmatrix}=\begin{bmatrix}0\\\ {-p^{*}}\end{bmatrix},$ (6)
$\displaystyle\therefore\ $
$\displaystyle{\lambda^{*}(t)}=\begin{bmatrix}{p_{c}}\\\
{-p_{c}t+q_{c}}\end{bmatrix},\ \ {\forall\ t\in[0,t_{f}^{*}]}$ (7)

where $p_{c},q_{c}\in\mathcal{R}^{3}$ are constants. The time-optimal control
$u^{*}(\cdot)$ necessarily minimizes the Hamiltonian evaluated along the
optimal state ($z^{*}(\cdot)$) and costate ($\lambda^{*}(\cdot)$)
trajectories, that is,

${\langle u^{*}(t),B^{T}\lambda^{*}(t)\rangle\leq\langle
v,B^{T}\lambda^{*}(t)\rangle\ \forall\ \\{v\in\mathcal{R}^{3}\ \mathrm{s.t.}\
\left\lVert v\right\rVert\leq\bar{u}\\}}$ (8)

where ${u^{*}(t)}$ is the optimal control input, and ${\lambda^{*}(t)}$ is the
solution of eq (6) that corresponds to the optimal trajectory ${z^{*}(t)}$ and
optimal input ${u^{*}(t)}$. Using eq (8), the optimal input is given by,

$\displaystyle{u^{*}(t)}=\begin{cases}-\frac{{B^{T}\lambda^{*}(t)}}{\left\lVert{B^{T}\lambda^{*}(t)}\right\rVert}\
\quad\text{if}\ {B^{T}\lambda^{*}(t)\neq 0}\\\
{v}\in\mathcal{U}\quad\quad\quad\quad\text{otherwise}\end{cases}$ (9)

In view of eq (7), the optimal input can be written as:

$\displaystyle{u^{*}(t)}=\begin{cases}{\bar{u}\frac{-p_{c}t+q_{c}}{\left\lVert-
p_{c}t+q_{c}\right\rVert}}\ \quad\text{if}\ {-p_{c}t+q_{c}\neq 0}\\\
{v}\in\mathcal{U}\quad\quad\quad\quad\text{otherwise}\end{cases}$ (10)

For the complete derivation of eq (10) please refer to [18].

###  II-B Navigation formula for time-optimally steering a point mass to a
specified position and velocity

Given the form of the minimum-time control input in eq (10) which was derived
using PMP, the corresponding time-optimal trajectory can be found for a given
set of boundary conditions. For convenience, the eqs in (1) can be non-
dimensionalized w.r.t. ${\bar{u}}$, i.e. $\mathrm{v}=v/\bar{u}$ and
$\mathrm{r}=r/\bar{u}$. Thus the non-dimensionalized set of eqs becomes:

$\displaystyle\begin{split}&\mathrm{\dot{r}(t)=v(t)};\quad\mathrm{\dot{v}(t)=u(t)+g}\\\
&\mathrm{\left\lVert u(t)\right\rVert\leq 1}\end{split}$ (11)

where, $\mathrm{g}=g/\bar{u}$ is the non-dimensionalized gravity.
The form of $u^{*}(t)$ is unaffected by gravity acceleration as shown in eq
(10). However, analytically solving for the optimal input $u^{*}(t)$ that
satisfies a given boundary condition, $\mathrm{col}(r_{0},v_{0}),\
\mathrm{col}(r_{f},v_{f})$ is non-trivial because of its form. This is why
most minimum-time solutions adopt alternative methods (numerical), [8, 10].
The analytical solution for the case of time optimal steering without gravity
has been solved in Akulenko et al, [1]. Here these eqs are extended to account
for gravity. For completeness, all eqs with the gravity terms below use the
same notation as [1]. We define the following parameters in order to solve the
optimal input for the system in eq (11):

$\displaystyle u:=\frac{-p_{0}t+q_{0}}{\left\lVert-p_{0}t+q_{0}\right\rVert}$
$\displaystyle\xi:=\frac{p_{0}}{\left\lVert
q_{0}\right\rVert};\quad\eta:=\frac{q_{0}}{\left\lVert q_{0}\right\rVert}$
$\displaystyle\rho:=\left\lVert\xi\right\rVert\neq
0;\quad\sigma:=\frac{\xi\cdot\eta}{\left\lVert\xi\right\rVert\left\lVert\eta\right\rVert}=\frac{\xi\cdot\eta}{\left\lVert\xi\right\rVert}\
.$ (12)

where $p_{0}$, $q_{0}\in\mathcal{R}^{3}$ are constant vectors. Using the
definitions in eqs (II-B), the optimal input from eq (10) is written as:

$\displaystyle\quad\quad\quad\quad\ u(t)=\frac{Q(t)}{R(t)},$ (13)
$\displaystyle\text{where }\ Q(t):=-\xi t+\eta;\quad
R(t):=\sqrt{\rho^{2}t^{2}-2\sigma\rho t+1}.$

Integrating the RHS of eqs in (11) results in:

$\displaystyle v(t)=v_{0}+\int_{0}^{t}\frac{Q(\tau)}{R(\tau)}d\tau+gt,$ (14)
$\displaystyle
r(t)=r_{0}+v_{0}t+\frac{1}{2}gt^{2}+\int_{0}^{t}\frac{(t-\tau)Q(\tau)}{R(\tau)}d\tau.$
(15)

For any instance of time $t$, the drone’s velocity and position in presence of
gravity can be written as:

$\begin{split}&v(t)=v_{0}+gt+V_{\xi}(t)\xi+V_{\eta}(t)\eta,\\\
&r(t)=r_{0}+v_{0}t+\frac{1}{2}gt^{2}+X_{\xi}(t)\xi+X_{\eta}(t)\eta,\\\
\end{split}$ (16)

where

$\begin{split}&S(t):=\log\left(k+\sqrt{1+k^{2}}\right);\ k=\frac{(\rho
t-\sigma)}{\sqrt{1-\sigma^{2}}},\\\
&V_{\xi}(t):=-\frac{1}{\rho^{2}}\left[\sigma
S(\tau)+R(\tau)\right]_{0}^{t};V_{\eta}(t)=\frac{1}{\rho}S(\tau)|_{0}^{t},\\\
&X_{\eta}(t):=\frac{1}{\rho^{2}}\left[-R(\tau)+(\rho\tau-\sigma)S(\tau)\right]_{0}^{t}-\frac{t}{\rho}S(0),\\\
&X_{\xi}(t):=\frac{1}{2\rho^{3}}\left[(\rho\tau+3\sigma)R(\tau)\right]_{0}^{t}+\frac{t}{\rho^{2}}(1+\sigma
S(0))\\\
&\quad\quad\quad\quad+\left[(-2\rho\sigma\tau+3\sigma^{2}-1)S(\tau)\right]_{0}^{t}.\\\
\end{split}$ (17)

The trajectory must start at time $t=0$ from $r_{0}$ with velocity $v_{0}$ in
the non-dimensional space and end at $r_{f}$ with non-dimensional velocity
$v_{f}$. Assuming the time taken to complete this path is given by $t_{f}$,
these boundary conditions can be substituted to obtain the eqs to be solved,
i.e.:

$\displaystyle-
r_{f}+r_{0}+v_{f}t_{f}-g\frac{t_{f}^{2}}{2}=t_{f}^{2}(a_{\zeta}\zeta+a_{\eta}\eta),$
(18) $\displaystyle v_{0}-v_{f}+gt_{f}=t_{f}(b_{\zeta}\zeta+b_{\eta}\eta),$
(19) where $\displaystyle\zeta:=\xi
t_{f},\quad\left\lVert\zeta\right\rVert:=\mu=\rho t_{f},$ $\displaystyle
a_{\zeta}=-\frac{(\mu+3\sigma)a+\mu+(3\sigma^{2}-1)b}{2\mu^{3}},$
$\displaystyle a_{\eta}=b_{\zeta}=\frac{a+\sigma b}{\mu^{2}};\
b_{\eta}=-\frac{b}{\mu},\ a=\sqrt{\mu^{2}-2\mu\sigma+1}-1,$ $\displaystyle
b=\log\left(\frac{\mu-\sigma+\sqrt{\mu^{2}-2\mu\sigma+1}}{1-\sigma}\right)\ .$

In order to solve eqs (18) and (19), we convert them into 3 scalar eqs by
taking dot product of the vectors with respect to each other, i.e.
LHS((18))$\cdot$LHS((18)) $=$ RHS((18)) $\cdot$ RHS((18)). Similar dot
products are applied between ((19),(19)) and ((18),(19)) giving:

$\displaystyle t_{f}^{2}f_{vv}=\left\lVert
v_{0}-v_{f}\right\rVert^{2}+\left\lVert
g\right\rVert^{2}t_{f}^{2}+2(v_{0}-v_{f})\cdot gt_{f}$ (20)
$\displaystyle\begin{split}&t_{f}^{4}f_{xx}=\left\lVert
r_{0}\right\rVert^{2}+\left\lVert
v_{f}\right\rVert^{2}t_{f}^{2}+\frac{\left\lVert
g\right\rVert^{2}}{4}t_{f}^{4}\\\ &\quad\quad\quad-g\cdot
v_{f}t_{f}^{3}-g\cdot r_{0}t_{f}^{2}+2x_{0}\cdot v_{f}t_{f}\\\ \end{split}$
(21)
$\displaystyle\begin{split}&f_{xv}t_{f}^{3}=r_{0}\cdot(v_{0}-v_{f})+(r_{0}\cdot
g)t_{f}\\\ &\quad\quad\quad+(v_{f}\cdot(v_{0}-v_{f}))t_{f}+(v_{f}\cdot
g)t_{f}^{2}\\\
&\quad\quad\quad-g\cdot(v_{0}-v_{f})\frac{t_{f}^{2}}{2}-\left\lVert
g\right\rVert^{2}\frac{t_{f}^{3}}{2}\\\ \end{split}$ (22)

By simplifying eqs (20)-(22), we obtain the following set of 3 scalar eqs that
need to be solved in terms of $\eta,\xi,t_{f}$ that determine the time-optimal
trajectory:

$\displaystyle\begin{split}&r_{0}\cdot(v_{0}-v_{f})=-(r_{0}\cdot
g)t_{f}-(v_{f}\cdot(v_{0}-v_{f}))t_{f}\\\
&\quad\quad\quad\quad\quad+g\cdot(v_{0}-3v_{f})\frac{t_{f}^{2}}{2}\\\
&\quad\quad\quad\quad\quad+\left(f_{xv}+\frac{\left\lVert
g\right\rVert^{2}}{2}\right)t_{f}^{3}\ \\\ \end{split}$ (23)
$\displaystyle\left\lVert
v_{0}-v_{f}\right\rVert^{2}=t_{f}^{2}(f_{vv}-\left\lVert
g\right\rVert^{2})-2((v_{0}-v_{f})\cdot g)t_{f}\ $ (24)
$\displaystyle\begin{split}&\left\lVert
r_{0}\right\rVert^{2}=t_{f}^{4}(f_{xx}-\frac{\left\lVert
g\right\rVert^{2}}{4})+(g\cdot v_{f})t_{f}^{3}\\\
&\quad\quad\quad-\left(\left\lVert v_{f}\right\rVert^{2}-g\cdot
r_{0}\right)t_{f}^{2}-2(r_{0}\cdot v_{f})t_{f}\end{split}$ (25)

where

$\displaystyle
f_{vv}=(b^{2}_{\zeta}\mu^{2}+b_{\eta}^{2}\sigma^{2}+2b_{\eta}b_{\zeta}\mu\sigma),$
$\displaystyle
f_{xx}=(a^{2}_{\zeta}\mu^{2}+a_{\eta}^{2}\sigma^{2}+2a_{\eta}a_{\zeta}\mu\sigma),$
$\displaystyle
f_{xv}=a_{\zeta}b_{\zeta}\mu^{2}+(a_{\zeta}b_{\eta}+a_{\eta}b_{\zeta})\mu\sigma+a_{\eta}b_{\eta}\
.$

Eqs (23)-(25) can be solved for the three scalar unknowns $\mu$, $\sigma$ and
$t_{f}$. The bounds of the variables are $\mu>0$, $-1<\sigma<1$ and $t_{f}>0$.
These three parameters can then be used to calculate $\eta,\xi$, from eqs (18)
and (19) which gives the full trajectory between two waypoints by substituting
them in eqs (16)

##  III Trajectory generation for surveying task

For the surveying task, we assume that the touring path is specified by the
user (zig-zag, spiral, circular, etc.). The zig-zag survey flight pattern is
used for our simulations and analysis. We consider coverage paths composed of
a sequence of waypoints and a prescribed upper bound on the drone speed at the
given waypoints. This is done to avoid motion blur in the images captured by
the drone when passing through each waypoint. The objective is to generate the
trajectory that utilizes maximum thrust from the quadrotor to minimize total
survey time. For this, we use our trajectory generator as the local planner
responsible for computing the intermediate keyframes (position, velocity,
acceleration and yaw).

Definition 3 For a RoI ‘$A$’, the optimal touring path for the survey is given
by $W(A)={w_{1},w_{2},...,w_{n}}$ where, $w_{k}\in\mathcal{R}^{3}\ \forall k$,
where $n$ is the total number of waypoints.

The time-optimal coverage trajectory generation problem can be formulated as:
Time optimal coverage trajectory generation for a region of interest ’A’ is

$\displaystyle{\min_{\Omega\in\mathcal{R}^{4}_{+}}\int_{0}^{t_{N}}dt}$
$\displaystyle{s.t.:}\ \ $ $\displaystyle{\dot{x}(t)=f(x,\Omega)}$ (26)
$\displaystyle{r(t_{i})=w_{i},\quad w_{i}\in}W(A)$ (27)

where ${x(t)=\mathrm{col}(r(t),v(t),q(t),\omega(t))}\in\mathcal{R}^{14}$ is
the state (column) vector, consisting of the position
${r(t)}\in\mathcal{R}^{3}$, velocity ${v(t)}\in\mathcal{R}^{3}$, the
orientation represented by the quaternion ${q(t)\in}\mathcal{R}^{4}$, and the
angular velocity ${\omega(t)\in}\mathcal{R}^{3}$. The input in such a system
is the vector ${\Omega(t)\in\mathcal{R}^{4}_{+}}$ where ${\Omega_{j}(t)}$ is
the speed of the ${j^{th}}$ rotor. The non-linear dynamics are expressed in
(26), where $f(x,\Omega)$ is the translational and rotational dynamics of a
quadrotor described in [17].

The constraint (27) corresponds to the imposed waypoint positions. Since the
quadrotor dynamics are differentially flat, [17], i.e. the states and the
inputs of the original system can be written as algebraic functions of flat
outputs and their derivatives, position ${r}\in\mathcal{R}^{3}$, and yaw angle
${\psi}$, we plan the trajectory in the flattened output space. We assume yaw
angle, ${\psi}$, is always 0. The limitations of the original system are
usually addressed by enforcing additional constraints such as maximum jerk, to
account for the thrust, angular velocity and angular acceleration limits of
the quadrotor. The optimization method presented here can be extended to
include such constraints. In this work, we consider only the particle dynamics
and consider a state vector, ${x(t)=\mathrm{col}(r(t),v(t))}$. Hence the
system dynamics ${f(x,\Omega)}$ simplifies to the system of eqs shown in (2).

###  III-A Waypoint Velocity Determination

The first step in computing the time-optimal trajectory in Problem III is to
find the optimal velocities the drone must have at each desired waypoint. We
then interpolate the optimal waypoint velocities to get a smooth trajectory.
To determine the optimal waypoint velocities, we discretize the domains of
Problem III. To approximate the optimal waypoint velocities, we assume the
input between two waypoints is of the form:

${u(t)=}\begin{cases}{\bar{u}\frac{\eta_{n,1}}{\left\lVert\eta_{n,1}\right\rVert}}&{t_{n}\leq
t\leq t_{n,s}}\\\
{\bar{u}\frac{\eta_{n,2}}{\left\lVert\eta_{n,2}\right\rVert}}&{t_{n,s}<t<t_{n+1}}\end{cases},$
(28)

where ${n\in\\{1,2,...,N\\}}$ is the waypoint index, ${t_{n}}$ denotes the
time at the ${n^{th}}$ waypoint and ${t_{n,s}}$ is the switching time where
the input switches from one constant vector to another. The control input can
be represented in a discrete form as follows:

$\displaystyle{u_{k}=\left[\bar{u}\frac{\eta_{1,1}}{\left\lVert\eta_{1,1}\right\rVert},\bar{u}\frac{\eta_{1,2}}{\left\lVert\eta_{1,2}\right\rVert},...,\bar{u}\frac{\eta_{N-1,1}}{\left\lVert\eta_{N-1,1}\right\rVert},\bar{u}\frac{\eta_{N-1,2}}{\left\lVert\eta_{N-1,2}\right\rVert}\right]^{T}}$
(29)

with variable discrete time steps:

$\displaystyle{dt_{k}=\left[t_{1,s}-t_{1},t_{2}-t_{1,s},...,t_{N-1,s}-t_{N-1},t_{N}-t_{N-1,s}\right]}$
(30)

where ${k\in\\{1,2,...,2(N-1)\\}}$. The velocities and positions at respective
time steps can be determined through integration:

$\displaystyle{r_{k+1}=r_{k}+v_{k}dt_{k}+u_{k}\frac{dt_{k}^{2}}{2}}$ (31)
$\displaystyle{v_{k+1}=v_{k}+u_{k}dt_{k}}$

Using discretized inputs and flight times, the time-optimal trajectory of
Problem III is cast as the NLP in Problem III-A. Discretized time-optimal
trajectory NLP

$\begin{gathered}\begin{aligned} &{\min\sum_{k=1}^{2N-1}dt_{k}}\\\ {s.t.}\ \ \
&{r_{k+1}=r_{k}+v_{k}dt_{k}+u_{k}\frac{dt_{k}^{2}}{2}}\\\
&{v_{k+1}=v_{k}+u_{k}dt_{k}}\\\ &{r_{j}=w_{i},\ \ i=1,...,N}\\\
&\quad\quad\quad\quad\ {j=1,3,...,2N-1}\\\ &{\left\lVert
u_{k}\right\rVert=\bar{u}}\\\ \end{aligned}\end{gathered}$

The method to solve such a switching point NLP is known in the literature. The
NLP problem here can be solved using an Interior Point Algorithm for Nonlinear
Optimization (IPOPT) [19] to find the required waypoint velocities.

###  III-B Optimal Trajectory Generation

After calculating the optimal velocities at each waypoint by solving Problem
III-A, we will use them to generate the continuous-time optimal input for the
drone. The continuous-time optimal input will generate the trajectory
(position, velocity, acceleration) passing through the prescribed waypoint
positions and with the calculated optimal waypoint velocities. The common
approach to generate the continuous trajectory is to directly integrate the
solution to Problem III-A between the time steps.

$\displaystyle{u(t)}$ $\displaystyle={u_{k}};\quad{v(t)}={v_{k}+u_{k}t}$ (32)
$\displaystyle{r(t)}$ $\displaystyle={r_{k}+v_{k}t+u_{k}\frac{t^{2}}{2}}$
$\displaystyle{k}$ $\displaystyle=\begin{cases}{0}&{0\leq
t\leq\sum_{k=1}^{1}dt_{k}}\\\ {1}&{\sum_{k=1}^{1}dt_{k}\leq
t\leq\sum_{k=1}^{2}dt_{k}}\\\ {\vdots}&{\vdots}\\\
{2N-2}&{\sum_{k=1}^{2N-2}dt_{k}\leq t\leq\sum_{k=1}^{2N-1}dt_{k}}\end{cases}$

While this approach is easier to solve, it does not give us any indication of
how far the solution is from the true continuous-time optimal solution. To
obtain the true time-optimal solution between waypoints, one needs to solve
the optimal control problem by utilizing PMP and computing the input that
minimizes the Hamiltonian at each instant along with the optimal state and co-
state trajectories. This can be done by using the NLP generated velocities at
each waypoint ($v_{0},v_{2},v_{4},...v_{2N-2}$). We use these velocities as
boundary conditions for the time-optimal steering problem whose solution was
presented in Section II

$\displaystyle{u^{*}(t)}$
$\displaystyle={\bar{u}\frac{-p_{i}t+q_{i}}{\left\lVert-
p_{i}t+q_{i}\right\rVert}}$ (33) $\displaystyle{v(t)}$
$\displaystyle={v_{i}+\int_{0}^{t}u^{*}(t)dt}$ $\displaystyle{r(t)}$
$\displaystyle={r_{i}+\int_{0}^{t}v(t)dt}$ $\displaystyle{i}$
$\displaystyle=\begin{cases}{0}&{0\leq t\leq\sum_{k=0}^{1}dt_{k}}\\\
{1}&{\sum_{k=0}^{1}dt_{k}\leq t\leq\sum_{k=0}^{2}dt_{k}}\\\ \vdots&\vdots\\\
{N-1}&{\sum_{k=0}^{N-1}dt_{k}\leq t\leq\sum_{i=0}^{N}dt_{k}}\end{cases}$

We solve the system of eqs in (33) in terms of the unknowns ${dt_{i},p_{i}}$
and ${q_{i}}$ for each set of waypoints. The optimal inputs obtained from
$p_{i}$, $q_{i}$’s steers the drone between two waypoints (${r_{i-1}}$ to
${r_{i}}$) with the start(${v_{i-1}}$) and end velocities(${v_{i}}$). The
${v_{i}}$’s come from solving the NLP in Problem III-A and ${r_{i}}$’s come
from waypoint locations ${w_{i}}$.

##  IV Simulation results

For our simulations, we assume the drone mass to be $1\ kg$, the maximum
magnitude of the thrust of the drone to be $10.5N$ and acceleration due to
gravity as $g=9.8m/s^{2}$. The set of optimal inputs and times ($dt_{i}$’s)
that solve eqs (33) were computed using the trust region dogleg algorithm in
MATLAB r2021b. Depending on the boundary conditions, the problem may be ill
conditioned for trust region based optimizers and multiple starting points are
needed for the optimizer. The search space was narrowed by using upper bounds
on time ($dt_{i}$’s) and bounds on $\mu$ and $\sigma$. Multiple solutions
existed even at residual tolerance values of order $10^{-6}$. To obtain a a
sufficiently accurate result, the residual was reduced to an order of
$10^{-11}$.

###  IV-A Comparison of interpolation methods

It is of interest to see how good is the solution obtained from discretizing
the search domain when solving for the optimal trajectory. The NLP solution of
Problem III-A produces the required waypoint velocities. To interpolate with
the optimal trajectory, we use those velocities as the boundary conditions for
the optimal control problem solved earlier using PMP(refer to Section II). The
entire domain was non-dimensionalized w.r.t. $\bar{u}$ and the minimum-time
solution based on PMP is computed. The resultant optimal positions, velocities
and accelerations are then scaled back to compare how good a direct
interpolation is.

Figures 2 and 3 show a comparison of the trajectory obtained from direct
interpolation v/s the trajectory obtained from solving the optimal control
problem. It can be seen that the trajectory interpolations are not very
different. Direct interpolation is easy to perform, incurs small computation
cost, and gives a solution which is very close to the optimal one. Hence a
discretized minimum-time control strategy produces a trajectory that is close
to the continuous minimum-time optimal solution. The time taken to go from a
waypoint $w_{i}\in W(A)$ to the next, $w_{i+1}\in W(A)$, is shown in in Figure
4. We observe that the times obtained from the NLP solution is not far from
the solution obtained from the PMP-based approach. While the time saved by
using the PMP approach is low, our method generates a much smoother trajectory
with a continuous input through all waypoints as shown in Figures 5 and 7.
This makes it easier to implement our method on real systems. Works such as
[14] propose using a time-optimal solution for the particle steering problem
as an initial guess for the time-optimal steering of a UAV with it’s rotary
dynamics included. The approach here provides a continuous-time optimal input
which potentially leads to better solutions after refinement for a UAV’s
underactuated dynamics.

The graph of the Hamiltonian v/s time shown in Figure 3 is a horizontal line
between two consecutive waypoints but is not constant throughout the
trajectory. This means that better waypoint velocities can exist that ensure
the Hamiltonian remains constant across all waypoints, which is a requirement
from PMP. That would make it the truly time-optimal solution. Refining the
waypoint velocity guesses is a possible direction for future work.

Figure 2: Paths obtained from direct interpolation vs continuous optimal input
solution (paths truncated for clarity) Figure 3: Velocity components from
direct interpolation and continuous optimal input solution Figure 4:
Comparison of time taken between consecutive waypoints using discretized and
continuous-time optimal solution,i.e. ${t(w_{i}\rightarrow w_{i+1})}$ vs
${w_{i}}$ Figure 5: Optimal input functions obtained from direct interpolation
and PMP approaches using 1 switching point

###  IV-B Influence of switching points

In our formulation, the input function is discretized by using switching
points. Only 1 intermediate switching point was used but it is possible to
have two or more switching points. The result of using more than one switching
point is shown in Figure 6. Using multiple switching points doesn’t greatly
affect the final objective function. The final time difference between using
one and multiple switching points is not significant. The inputs inferred from
those velocities are nearly identical and appear slightly smoother near the
waypoints. This can be seen when comparing Figure 5 and Figure 7. Clearly,
using more switching points increases the computation time of the solution to
the NLP. For a large number of switching points, some solvers had difficulty
converging to a solution. To reduce computation time, it is better to use one
switching point. If the requirement is to have a smoother time-optimal input,
then multiple switching points may be used.

Figure 6: Using more switching points doesn’t significantly impact the total
flight time Figure 7: Optimal input functions obtained from direct
interpolation and PMP approaches using 4 switching points

###  IV-C Comparison with alternative trajectory generation methods

For a comparison of our algorithm to existing methods, we use the trajectory
generated from the minimum snap problem described in [17, 15]. Most
conventional UAV planning methods use a polynomial trajectory that minimizes
jerk or snap. This is done to make the derivatives smooth and minimize body
torques, but such polynomial planners do not ensure that the drone uses
maximum actuator capacity/acceleration at all trajectory points, as seen in
Figures 8 and 9. Our approach produces a continuous input function while
ensuring that the boundary of the input value set is achieved at all points in
the trajectory, as seen in Figure 9. As a result, the total time taken by the
minimum snap planner to complete the trajectory is 20% slower than our method.

Figure 8: Comparison of minimum snap trajectory (solid) vs our trajectory
generator (dotted) Figure 9: Comparison of modified min. snap trajectory
(solid) vs the NLP trajectory generator (dotted)

##  V Conclusion

In this work we presented a time-optimal surveying planner that generates a
continuous trajectory which uses maximum thrust at all points and minimizes
the aggregate flight time. Solution optimality is verified by comparing the
generated trajectory to the solution obtained by solving the time-optimal
control problem between two waypoints. We did this by discretizing and
formulating a NLP problem for time-optimal surveying. We then use that to
generate position and velocity references that can be directly interpolated to
generate a continuous trajectory. We show that a direct interpolation is not
far off from the optimal control input trajectory obtained by solving for the
PMP solution. Next we study the effect of using multiple switching points and
how much improvement can be obtained from increasing them. No significant time
savings were observed on adding more intermediate switching points. However,
the continuous-time optimal input that was inferred using multiple switching
points was slightly smoother.

While this work focuses on aggregate flight time, total computation time is
another aspect that can be optimized for reduction in time. This can be
addressed in a future work. Possible approaches include the use of a better
initial guess(convexifying the NLP and using that solution) or using data-
driven approaches, [20], that learn from the planner described in this work.

##  VI Acknowledgement

This work was funded by the Rig Automation and Performance Improvement in
Drilling (RAPID) Consortium of University of Texas at Austin.

## References

  * [1] LD Akulenko and AP Koshelev  “Time-optimal steering of a point mass to a specified position with the required velocity”  In _Journal of Applied Mathematics and Mechanics_ 71.2 Elsevier, 2007, pp. 200–207
  * [2] LD Akulenko and AM Shmatkov  “Time-optimal steering of a point mass onto the surface of a sphere at zero velocity”  In _Journal of Applied Mathematics and Mechanics_ 66.1 Elsevier, 2002, pp. 9–21
  * [3] Aneesh Venkatraman and Sanjay P Bhat  “Optimal planar turns under acceleration constraints”  In _Proceedings of the 45th IEEE Conference on Decision and Control_ , 2006, pp. 235–240
  * [4] Efstathios Bakolas and Andrei Marchidan  “Time-optimal control of a self-propelled particle in a spatiotemporal flow field”  In _International Journal of Control_ 89.3 Taylor & Francis, 2016, pp. 623–634
  * [5] Efstathios Bakolas and Panagiotis Tsiotras  “Optimal synthesis of the Zermelo–Markov–Dubins problem in a constant drift field”  In _Journal of Optimization Theory and Applications_ 156.2 Springer, 2013, pp. 469–492
  * [6] Ernest Bruce Lee and Lawrence Markus  “Foundations of Optimal Control Theory”, 1967
  * [7] Zhong-Hui Du, Huai-Ning Wu and Shuang Feng  “Boundary Control of a Quadrotor UAV with a Payload Connected by a Flexible Cable”  In _2018 37th Chinese Control Conference (CCC)_ , 2018, pp. 1151–1156  IEEE
  * [8] Li-Chun Lai, Chi-Ching Yang and Chia-Ju Wu  “Time-optimal control of a hovering quad-rotor helicopter”  In _Journal of Intelligent and Robotic Systems_ 45.2 Springer, 2006, pp. 115–135
  * [9] Wannes Van Loock, Goele Pipeleers and Jan Swevers  “Time-optimal quadrotor flight”  In _2013 European Control Conference (ECC)_ , 2013  IEEE
  * [10] Robin Ritz, Markus Hehn, Sergei Lupashin and Raffaello D’Andrea  “Quadrocopter performance benchmarking using optimal control”  In _2011 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2011
  * [11] Sergei Lupashin and Raffaello D’Andrea  “Adaptive Open-Loop Aerobatic Maneuvers for Quadrocopters” 18th IFAC World Congress  In _IFAC Proceedings Volumes_ 44.1, 2011
  * [12] Xiaohui Zhou et al.  “Survey on path and view planning for UAVs”  In _Virtual Reality and Intelligent Hardware_ 2.1, 2020, pp. 56–69
  * [13] Giuseppe Loianno, Davide Scaramuzza and Vijay Kumar  “Special issue on high-speed vision-based autonomous navigation of uavs”  In _Journal of Field Robotics_ 1.1 Wiley-Blackwell Publishing, Inc., 2018, pp. 1–3
  * [14] Philipp Foehn, Angel Romero and Davide Scaramuzza  “Time-optimal planning for quadrotor waypoint flight”  In _Science Robotics_ 6.56 American Association for the Advancement of Science, 2021
  * [15] Gilhyun Ryou, Ezra Tal and Sertac Karaman  “Multi-fidelity black-box optimization for time-optimal quadrotor maneuvers”  In _arXiv preprint arXiv:2006.02513_ , 2020
  * [16] Tauã M. Cabreira, Lisane B. Brisolara and Paulo R. Ferreira Jr.  “Survey on Coverage Path Planning with Unmanned Aerial Vehicles”  In _Drones_ 3.1, 2019
  * [17] Daniel Mellinger and Vijay Kumar  “Minimum snap trajectory generation and control for quadrotors”  In _2011 ICRA_ , 2011, pp. 2520–2525  IEEE
  * [18] Michael Athans and Peter L Falb  “Optimal control: an introduction to the theory and its applications”  Courier Corporation, 2013
  * [19] Andreas Wächter and Lorenz T Biegler  “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming”  In _Mathematical programming_ 106.1 Springer, 2006, pp. 25–57
  * [20] Marcelino M Almeida, Rahul Moghe and Maruthi Akella  “Real-time minimum snap trajectory generation for quadcopters: Algorithm speed-up through machine learning”  In _2019 ICRA_ , 2019