# Nonlinear Controllability Assessment of Aerial Manipulator Systems using
Lagrangian Reduction

Skylar X. Wei  Matthew R. Burkhardt  Joel Burdick  California Institute of
Technology, Pasadena, CA 91125 USA (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>jwb@robby.caltech.edu).

###### Abstract

This paper analyzes the nonlinear Small-Time Local Controllability (STLC) of a
class of underatuated aerial manipulator robots. We apply methods of
Lagrangian reduction to obtain their lowest dimensional equations of motion
(EOM). The symmetry-breaking potential energy terms are resolved using
advected parameters, allowing full $SE(3)$ reduction at the cost of additional
advection equations. The reduced EOM highlights the shifting center of gravity
due to manipulation and is readily in control-affine form, simplifying the
nonlinear controllability analysis. Using Sussmann’s sufficient condition, we
conclude that the aerial manipulator robots are STLC near equilibrium
condition, requiring Lie bracket motions up to degree three.

###### keywords:

Controllability; Geometric Mechanics; Lagrangian Reduction; Flying robots;
UAVs.

††thanks: This work was funded in part by a contract from the NASA Jet
Propulsion Laboratory, and by a NASA NSTRF fellowship for the second author.

##  1 Introduction

Aerial manipulators, which combine a multi-rotor aerial platform with a multi-
jointed robot arm, have the mobility of multi-rotors and can interact with the
environment utilizing its manipulator. A variety of aerial manipulators have
been introduced, ranging from a fully-actuated aerial base with a three-link
manipulator (Sanchez-Cuevas et al.(2020)) to a helicopter with a seven degree-
of-freedom manipulator (Kondak et al. (2014)). As summarized by Ruggiero et
al. (2018), no matter the complexity of the fully-actuated appended arm, the
underactuation of the aerial base poses limitations on the overall system
ability.

Despite many successful designs, there is no formal guideline for designing an
aerial manipulator based on task objectives, nor has there been any analytical
assessment of the limitation inherent in an underactuated base. To the best of
the authors’ knowledge, there has not been any theoretical work studying the
underactuated aerial system’s controllability properties. Controllability,
which describes a system’s ability to drive its states to arbitrary values by
choices of feasible inputs, is widely used in the study of stabilization,
feedback controller design, and state-space reduction. Nonlinear
controllability analysis is more challenging than the linear case, but it must
be considered for aerial manipulators to ensure sufficient manipulability.

To tackle the algebraic complexity of nonlinear controllability analysis, we
establish a lower-dimensional EOM by utilizing Lagrangian reduction. This
reduction technique yields a set of the first-order EOM on a lower-dimensional
phase space that excludes the symmetry group, which is algebraically simpler
in the subsequent controllability assessment. Methods of Lagrangian reduction
have successfully been applied to evaluating controllability of complex
undulatory robots (Ostrowski (1996)) and the modeling and control of spherical
robotic vehicles (Burkhardt (2018)). We consult our earlier work (Burkhardt
and Burdick (2016)) to tackle the symmetry-breaking potential terms of aerial
manipulators via advected parameters.

We analyze nonlinear controllability under small-time, local, and near-
equilibrium conditions, using the sufficient conditions given by Sussmann
(1987), which applies to systems with a control-affine form, a non-zero drift,
and bilateral (both positive and negative) control inputs. However, aerial
manipulators have unilateral control inputs (non-negative thrust), making them
ill-suited for a direct application. Based on the rotor dynamics, we select
bilateral rotor RPM rate as the replacement control inputs. Our
controllability result can be viewed as a certification of an aerial
manipulator’s fundamental capability to control the multirotor platform during
aerial manipulation. The set of reduced EOM is mathematically compact and
reveals geometric properties that could serve in the motion planning and
controller synthesis process.

Paper Outline. Section 1 defines the class of aerial manipulator system
studied. Section 3 reviews the basics of reduction and reconstruction under
broken symmetries, as well as controllability of nonlinear systems. The EOM
for a general class of aerial manipulator system is developed in Section 4.
Section 5 presents the main contribution: proofs that our class of aerial
manipulators are Small-Time Locally Accessible and Small-Time Locally
Controllable. We conclude in Section 6 with future applications.

##  2 Systems Description

We consider the following class of Aerial Manipulators:

Definition: A multirotor aerial platform with the following characteristics is
an Aerial Manipulator (AM):

  * •

The multi-rotor includes $n$-pairs of identical rotors where $n\geq 2$. Each
rotor pair consists of one clockwise (CW) and one counterclockwise (CCW)
rotating rotor, distributed in a cross configuration.

  * •

All rotor thrusts point in the positive direction of the $z_{b}$ axis, and the
distances from their rotating axis to the $z_{b}$ axis are equal, see Fig. 1.

  * •

The 2-link manipulator is attached to the multi-rotor’s geometric center and
operates in $x_{b}-z_{b}$ plane. Each link is approximated by a uniform
cylinder.

  * •

All system components are rigid and complex fluid structure interactions are
ignored.

The manipulator is limited to planar operation only for simplicity: it is
sufficiently complicated to demonstrate the effectiveness of the reduction
process and the dynamic coupling between the multi-rotor and arm dynamics.

Our model is derived using the following reference frames:

  * •

The earth-fixed inertial frame $E=\\{O^{e},x_{e},y_{e},z_{e}\\}$.

  * •

The aerial-base body frame $B=\\{O^{b},x_{b},y_{b},z_{b}\\}$.

  * •

Manipulator $i^{th}$ link frame
$L_{i}=\\{O^{L_{i}},x_{L_{i}},y_{L_{i}},z_{L_{i}}\\}$.

  * •

The $j^{th}$ rotor frame
$T_{j}=\\{O^{T_{j}},x_{T_{j}},y_{T_{j}},z_{T_{j}}\\}$.

Hereafter, $s_{ab}\in\mathbb{R}^{3}$ and $R_{ab}\in SO(3)$ denote the position
and orientation of the origin of frame $B$ relative to the origin of frame
$A$, respectively. $\eta_{1},\eta_{2}\in\mathbb{S}$ be the relative joint
angles. The aerial-base linear velocity in $B$ frame is defined to as
$\dot{s}_{b}\triangleq R_{eb}^{T}\dot{s}_{eb}$.

Figure 1: Pictorial illustration of relevant reference frames for the defined
aerial manipulator system.

Using standard roll, pitch, and yaw (Euler) angles
$\mathbf{\Theta}=[\phi,\theta,\psi]^{T}$, the angular velocity of the aerial-
base in $B$ frames, $\omega_{b}$, is the following:

$\omega_{b}=\left[\begin{array}[]{ccc}1&0&-\sin(\theta)\\\
0&\cos(\phi)&\sin(\phi)\cos(\theta)\\\
0&-\sin(\phi)&\cos(\phi)\cos(\theta)\end{array}\right]\left[\begin{array}[]{c}\dot{\phi}\\\
\dot{\theta}\\\
\dot{\psi}\end{array}\right]=T(\mathbf{\Theta})\dot{\mathbf{\Theta}}.$ (1)

Lastly, $\dot{R}_{eb}\triangleq R_{eb}S(\omega_{b})$ where $S(\cdot)$ is the
$3\times 3$ skew-symmetric matrix such that $S(\omega_{b})x=\omega_{b}\times
x$.

##  3 Preliminaries

A mechanical system is characterized by the tuple $\Sigma=(Q,L,\mathcal{T})$,
where $Q$ is its finite-dimensional configuration space. $TQ$ denotes the
tangent bundle of $Q$, and $T_{q}Q$ is the tangent space to $Q$ at $q$.
$L:TQ\to\mathbb{R}$ is the system Lagrangian. $\mathcal{T}(q,\dot{q})\in
T^{*}Q$ represents the external forces field acting on the system, where
$T^{*}Q$ is the dual of $TQ$.

Broadly speaking, a Lagrangian possesses a symmetry if there is an action on
its arguments that renders the Lagrangian invariant. This symmetry allows the
reduction of the dynamical system to a lower dimensional phase space. This
paper is concerned with Lie group symmetries, which naturally arise in rigid
body systems. For the rest of the paper, we will assume that $Q$ is a trivial
principle bundle, i.e., $Q=G\times M$, where the fiber-space $G$ is a Lie
group. Smooth manifold $M$ denotes the internal shape space. Consequently,
coordinates on $q\in Q$ can be partitioned as $q=(\mathtt{g},r)$, where
$\mathtt{g}\in G$ and $r\in M$. Associated with a Lie group is its Lie
algebra, $\mathfrak{g}$, a vector space isomorphic to the tangent space at the
group identity, i.e. $\mathfrak{g}\simeq T_{q}G$. In the context of reduction,
we are interested in the left action of a Lie group $X$ on a smooth manifold
$P$, which is defined as the map $\Phi_{\mathtt{x}}:P\times P\to
P:p\to\mathtt{x}p$ for any $p\in P$. The Lie group of an aerial manipulator
robot’s configuration space is a semidirect product space, defined as
$G=H\otimes\mathcal{V}$ where $H$ is a Lie subgroup and $\mathcal{V}$ is a
vector space. We denote the Lie algebra of $H$ as $\mathfrak{h}$, and the
group $H$ acts on the space $\mathcal{V}$ from the left, and on $Q$ via a
left-action. The Lie algebra of the semidirect-product group can be written as
$\mathfrak{g}=\mathfrak{h}\otimes
T\mathcal{V}\simeq\mathfrak{h}\otimes\mathcal{V}$ with elements
$(\xi_{\mathtt{h}},\xi_{\mathcal{V}})\in\mathfrak{g}$. In local/body
coordinates, $\xi_{\mathtt{h}}=\mathtt{h}^{-1}\dot{\mathtt{h}}\in\mathfrak{h}$
and $\xi_{\mathcal{V}}=\mathtt{h}^{-1}\dot{v}\in T\mathcal{V}$ where
$(\dot{\mathtt{h}},\dot{v})\in TG$ is an arbitrary tangent vector. The lifted
action results from the left-action $\Phi_{\mathfrak{g}}(q)$ on tangent
vectors is defined as the map $T\Phi_{\mathfrak{g}}:TQ\to
TQ:\,(q,v)\to(\Phi_{\mathtt{g}}(q),T_{q}\Phi_{\mathtt{g}}(v))$. For
$\mathtt{g}=(\mathtt{h},v)$,
$T_{q}\Phi_{\mathtt{g}}(\dot{l},\dot{p},\dot{r})=(\mathtt{h}\dot{l},\mathtt{h}\dot{p},\dot{r})$.

A mechanical system possess a symmetry with respect to Lie group $G$ if its
Lagrangian $L:TQ\to\mathbb{R}$ and external forces $F(q,\dot{q})$ are Lie
Group Invariant. See Schwarz (1975). As the appended arm moves during its
manipulation tasks, the system’s overall center of mass displaces, breaking
symmetry in the Lagrangrian’s potential. Thus, we consider the mechanics under
a broken symmetry condition, which leads to a modified Euler-Poincaré set of
equations as will be discussed next.

###  3.1 Lagrangian Reduction and reconstruction

The reduction technique based on advected parameter, introduced by Holm et al.
(1998), is employed to reformulate the symmetry breaking potential
contributions. In the context of mechanical system, an advected parameter,
$\gamma(t)$, is a vector expressed in a body-fixed reference frame satisfying
the following differential equation:

$\left(\frac{d}{dt}+\mathtt{g}^{-1}(t)\dot{\mathtt{g}}(t)\right)\gamma(t)=0.$
(2)

where $\mathtt{g}\in G$. For AM, we choose an advected parameter,
$\gamma(t)\triangleq R_{eb}^{T}e_{3}\in\mathcal{V}$, to address the symmetry
breaking gravity. Another symmetry breaking term of the AM is its dependency
on the aerial base position. Let $\zeta\triangleq R_{eb}^{T}s_{eb}$, together
with $\gamma$, the AM potential energy can be expressed as
$V(r,\gamma,\zeta)$, where $r$ denotes the shape variables.

Define the Augmented Lagrangian
$\overline{L}:TG\times\mathcal{V}\to\mathbb{R}$ by augmenting the state with
advected parameter $\gamma$ and base position $\zeta$. If the Augmented
Lagrangian is $\Phi_{\mathtt{h}}$-invariant, then it can be reduced to
$\mathfrak{h}\times\mathcal{V}\times TM$ (Burkhardt (2018)). Applying the
Lagrange-d’Alembert principle to the Augmented Lagrangian yields a modified
form of the Euler-Poincaré equations. See Schneider (2002) for details. One
important takeaway from Schneider’s result is the following reduced
variational principle:

$\delta\left(\int_{t_{0}}^{t_{f}}K(\xi_{\mathtt{h}},\xi_{\mathtt{\mathcal{V}}},r,\dot{r})dt-\int_{t_{0}}^{t_{f}}V(r,\gamma,\zeta)\right)=0.$
(3)

Merging (3) with the Lagrange-d’Alembert principle yields:

$\delta\int_{t_{0}}^{t_{f}}K(\xi_{\mathtt{h}},\xi_{\mathtt{\mathcal{V}}},r,\dot{r})dt=\delta\int_{t_{0}}^{t_{f}}dV(r,\gamma,\zeta)\cdot\delta
qdt,$ (4)

where $-dV(r,\gamma,\zeta)$ represents conservative forces and moments arise
from the potential energy.

Based on Bullo and Lewis (2004), we can conclude the dynamical equivalence
between the mechanical systems $\Sigma=(Q,L,\mathcal{T})$ and the reduced
system $\Sigma_{r}=(Q,K,-dV+\mathcal{T})$. At the expanse of two additional
equations for $\dot{\gamma}$ and $\dot{\zeta}$, reduction by $G$ ($SE(3)$ for
AMs) is made possible on an augmented configuration space.

###### Theorem 1

(Burkhardt (2018)) Let a mechanical system with Lagrangian $L:=K-V$ satisfy
the following:

  1. 1.

The configuration space $Q$ has the product structure $Q=G\times
M=H\otimes\mathcal{V}\times M$, where $H$ is a Lie group, with a left action
on vector space $\mathcal{V}$, and $M$ is the shape space. For AM, $H=SO(3)$,
$\mathcal{V}=\mathbb{R}^{3}$, and $M$ is a smooth manifold. Hence,
$(q,\dot{q})=(h,\dot{h},\xi,\dot{\xi},r,\dot{r})$.

  2. 2.

The kinetic energy $K(q,\dot{q})$ and Augmented Lagrangian,
$\overline{L}=TG\times\mathcal{V}\to\mathbb{R}$, are invariant with respect to
a left $G$-action.
$L(\Phi_{\mathtt{g}}(q),T_{q}\Phi_{\mathtt{g}}\dot{q})=L(q,\dot{\mathtt{g}})$.

  3. 3.

The Reduced Augmented Lagrangian can be expressed as the following function
form:
$\overline{l}(\xi_{\mathtt{h}},\xi_{\mathcal{V}},\gamma,\zeta,r,\dot{r})$
where $r,\dot{r}\in M$, $\xi_{\mathtt{h}}\in\mathfrak{h}$,
$\xi_{\mathcal{V}}\in\mathcal{V}$. $\gamma$ and $\zeta$ are the advected
parameter and aerial base position.

The system’s reduced equations of motion takes the form:

$\dot{\varrho}=\dot{r}^{T}\alpha(r,\gamma)\dot{r}+\dot{r}^{T}\beta(r,\gamma)\varrho+\varrho^{T}\kappa(r,\gamma)\varrho+\tau_{\varrho}(q,\gamma),$
(5) $M(r)\ddot{r}=-C(r,\dot{r})+N(r,\dot{r},\varrho)+\tau_{r}(r,\gamma),$ (6)
$\dot{\gamma}=-\xi_{h}\gamma,\quad\quad\dot{\zeta}=-\xi_{h}(\zeta-A(r,\gamma)),$
(7)

where $\varrho$ are momenta (defined below), $\tau_{\varrho}(r,\gamma,\zeta)$
represents the conservative forces and moments resulting from gravity
projected along the momenta directions. The functions $\alpha(r,\gamma)$,
$\beta(r,\gamma)$, and $\kappa(r,\gamma)$ are smoothly dependent upon the
shape variables $r$ and advected parameter $\gamma$. $M(r)$ is the reduced
mass-matrix, $C(r,\dot{r})$ contains the Coriolis terms,
$N(r,\dot{r},\varrho)$ is the potential contribution, and $\tau_{r}(r,\gamma)$
represents the actuation forces associated with shape variables. $A(r,\gamma)$
is an invertible mass-like matrix which will arise in the reconstruction
equation, defined later (9). The structure of (5) and (6) follows from the
reduced variational principle with an extended base space consisting of the
generalized momenta $\varrho=(\rho_{1},\cdots,\rho_{n})$ and shape-space
variables $(r,\dot{r})$. (5) - (7) are the momentum equation, shape dynamics,
advection equation, and position dynamics of the system, respectively.
Together, they form a complete, reduced representation of the system dynamics.

The momenta are defined as follows. Recall the group orbit for the left action
$\Phi_{\mathtt{g}}$(q) on a manifold $Q$ is the set
$\mbox{Orb}(q)=\\{\Phi_{\mathtt{g}}(q)|\forall\mathtt{g}\in G\\}$. We denote
the tangent space to the orbit at $q\in Q$ as $T_{q}\mbox{Orb}(q)$. When $G$
is the Lagrangian symmetry group, the tangent space $T_{q}\mbox{Orb}(q)$
represents velocities along the symmetry directions. Let ${b_{i}}$ denote any
basis of $T_{q}\mbox{Orb}(q)$. For $SE(3)$, we choose the basis

$T_{q}\mbox{Orb}(q)=\mbox{span}\left(\frac{\partial}{\partial
s_{bx}},\frac{\partial}{\partial s_{by}},\frac{\partial}{\partial
s_{bz}},\frac{\partial}{\partial\omega_{b_{x}}},\frac{\partial}{\partial\omega_{b_{y}}},\frac{\partial}{\partial\omega_{b_{z}}}\right).$
(8)

The generalized momenta $\varrho_{i}$ are defined as
$\varrho_{i}(\xi)=\left<b_{i},\xi\right>$ along the symmetry directions. To
recover the spatial motion of the system, we employ a reconstruction equation
or connection, using trajectories in the reduced phase space. It defines a
horizontal space of $T_{q}Q$ as $H_{q}Q\triangleq\mbox{Ker}(\mathcal{A}(q))$,
where $\mathcal{A}$ is a principle connection form and describes motion along
the fiber as the flow of a left-invariant vector field. See Bloch et al.
(1996) and Ostrowski (1996). The general form of connection is:

$\xi_{\mathtt{h}}=\mathtt{h}^{-1}\dot{\mathtt{h}}=-A(r,\gamma)\dot{r}+B^{-1}(r,\gamma)\varrho,$
(9)

where $A(r,\gamma)$ and $B(r,\gamma)$ are mass and inertia tensors.

###  3.2 Nonlinear Control Theory Review

Consider steering problems in control-affine form:

$\Sigma:\dot{q}=f(q)+\sum_{i=1}^{m}g_{i}(q)u_{i}\quad
u\in\mathcal{U}\subset\mathbb{R}^{m},$ (10)

from an initial state $q_{0}\in Q$ to a final state $q_{f}\in Q$ by controls
$u:[0,T]\to\mathcal{U}$. $f$ is the drift vector field. $g_{i}$,
$i\in\\{1,...,m\\}$, are actuation vector fields which assumed to be real-
valued smooth functions. Define a distribution $\Delta$ of the system
$\Sigma$: $\Delta\triangleq\mbox{span}\\{g_{1},\cdots,g_{m}\\}$. The
followings are fundamentals for analyzing $\Sigma$’s controllability.

Definition: The reachable set, denoted as $\mathcal{R}^{V}(q_{0},\geq T)$, is
the set of all points $q_{f}$ such that there exists an input
$u(t)\in\mathcal{U}$, $0\leq t\leq T$, that steers the system (10) from
$q(0)=q_{0}$ to $q(T)=q_{f}$, where the trajectories $q(t)$, $0\leq t\leq T$,
remain inside a neighborhood $V$ of $q_{0}\in Q$.

Definition: $\Sigma$ is Small-Time Locally Accessible (STLA) from $q_{0}$ if
the reachable set contains a full n-dimensional subset of $Q$ $\forall$
neighborhoods $V$ and all $T>0$.

Definition: $\Sigma$ is Small-Time Locally Controllable (STLC) from $q_{0}$ if
the reachable set contains a neighborhood of $q_{0}$ $\forall$ neighborhoods
$V$ and all $T>0$.

The Lie bracket (product) between vector fields $g_{1}(q)$ and $g_{2}(q)$ is
$[g_{1}(q),g_{2}(q)]\triangleq\frac{\partial g_{1}(q)}{\partial
q}g_{2}(q)-\frac{\partial g_{2}(q)}{\partial q}g_{1}(q)$ , which quantifies
how the derivative of vector field $g_{2}(q)$ varies along the flow of
$g_{1}(q)$. The Lie bracket $[g_{1}(q),g_{2}(q)]$ (and higher order brackets)
may enable infinitesimal movements locally in the system tangent space along
directions that are not in $\Delta$. In summary, STLA is to access the Lie
brackets’ ability to generate independent actuation vector fields. STLC is to
further examined the brackets’ ability to overcome the drift $f$, leading to
the following definition.

Definition: A Lie product is considered to be a bad bracket if the drift term
$f$ appears an odd number of times in the product and each control vector
field $g_{i},i\in\\{1,...,m\\}$, appears an even number of times (including
zero). If a Lie product is not bad, it is a good bracket.

###### Theorem 2

(Hermann (1963)) The system (10) is STLA from $q$ if it satisfies the Lie
algebra rank condition (LARC):
$\overline{\mbox{Lie}}({f,g_{1}...,g_{m}})(q)=T_{q}Q$, where
$\overline{\mbox{Lie}}(\cdot)$ is the closure of Lie algebra over the span of
all input vector fields and their iterated Lie brackets,
$\\{f(q),g_{1}(q),\cdots,g_{m}(q)\\}$.

For driftless systems, i.e. $f(q)=0,\forall q\in Q$, STLA and STLC properties
are equivalent, and LARC alone is sufficient to establish both properties.
However, for systems with drift like an aerial manipulator robot, the
following condition, proposed by Sussmann (1987), is instead sufficient to
conclude STLC.

###### Theorem 3

(Sussmann (1987)) System (10) is STLC from $q^{*}$ if (1) $f(q^{*})=0$, (2)
$\mathcal{U}$ is bilateral, (3) the LARC is satisfied by good Lie bracket
terms up to degree $k$, and (4) any bad bracket of degree $j\leq k$ can be
expressed as a linear combination of good brackets of degree less than $j$.

##  4 System Dynamics

Suppose $m_{b}$ and $I_{b}$ are the mass and inertia tensor of the bare multi-
rotor structure, expressed in frame $B$. The kinetic energy of the this rigid
body is
$K_{b}=\frac{1}{2}m_{b}\dot{s}_{eb}^{T}\dot{s}_{eb}+\frac{1}{2}\omega_{b}^{T}I_{b}\omega_{b}$
, and potential energy is $V_{s}=m_{b}ge_{3}^{T}s_{eb}$, where vectors
$e_{1}$, $e_{2}$, and $e_{3}$ be standard Cartesian basis vectors.

To properly capture the yaw dynamics and control, we separately model the
rotors’ dynamics. Explicitly, the translational kinetic energy of the $j^{th}$
rotor $K_{t,t_{j}}=\frac{1}{2}m_{r}\dot{s}_{et_{j}}^{T}\dot{s}_{et_{j}}$ where
$m_{r}$ denotes the mass of an individual rotor, and its rotational kinetic
energy is

$K_{r,t_{j}}=\frac{1}{2}\left(J\Omega_{j}^{2}+\omega_{b}^{T}I_{r}\omega_{b}+2J\Omega_{j}e_{3}^{T}\omega_{b}\right),$
(11)

where $I_{r}\triangleq diag(a,a,J)$ is the moment of inertia of the $j^{th}$
rotor w.r.t. frame $T_{j}$ , and rotors are modeled as cylinders.

Suppose the mass, length, and inertia tuple of the link 1 and 2 w.r.t. frame
$L_{1}$ and $L_{2}$ are $\\{m_{1},d_{1},I_{1}\\}$ and
$\\{m_{2},d_{2},I_{2}\\}$, respectively. The 2-link manipulator dynamics can
be conveniently expressed using the manipulator Jacobian matrix:

$\left[\begin{array}[]{cc}\dot{s}_{bl_{1}}&\omega_{l_{1}}^{b}\end{array}\right]^{T}=\left[\begin{array}[]{cc}J_{l,bl_{1}}^{b}(\eta_{1})&J_{r,bl_{1}}^{b}(\eta_{1})\end{array}\right]^{T}\dot{\eta}_{1},$
(12)
$\left[\begin{array}[]{cc}\dot{s}_{bl_{2}}&\omega_{l_{2}}^{b}\end{array}\right]^{T}=\left[\begin{array}[]{cc}J_{l,bl_{2}}^{b}(\eta_{1},\eta_{2})&J_{r,bl_{2}}^{b}(\eta_{1},\eta_{2})\end{array}\right]^{T}\left[\begin{array}[]{cc}\dot{\eta_{1}}&\dot{\eta_{2}}\end{array}\right]^{T},$
(13)

where $J_{bl_{i}}^{b}(\eta_{1},\cdots,\eta_{i})\in\mathbb{R}^{6\times i}$ is
the manipulator spatial Jacobian matrices defined in Murray et al. (1994).
Similar to the multi-rotor case, the kinetic and potential energy of each link
can be calculated by summing the translational and rotational contributions.
Details will be omitted.

###  4.1 Overall System Dynamics

Let
$\dot{q}=[\dot{s}_{b}^{T},\omega_{b}^{T},\Omega_{j},\dot{\eta}_{1},\dot{\eta}_{2}]^{T}\in\mathbb{R}^{8+2n}$.
The kinetic energy of the AM system can be rewritten as the following:

$K_{AM}=\frac{1}{2}\dot{q}^{T}\mathcal{M}(\eta_{1},\eta_{2})\dot{q}=\frac{1}{2}\dot{q}^{T}\left[\begin{array}[]{c|c|c|c}\mathcal{M}_{p}&\mathcal{M}_{p\omega}&0_{3\times
2n}&\mathcal{M}_{pl}\\\
\hline\cr\mathcal{M}_{p\omega}^{T}&\mathcal{M}_{\omega}&\mathcal{M}_{\omega
r}&\mathcal{M}_{\omega l}\\\ \hline\cr 0_{2n\times 3}&\mathcal{M}_{\omega
r}^{T}&\mathcal{M}_{r}&\mathbf{0}_{2n\times 2}\\\
\hline\cr\mathcal{M}_{pl}^{T}&\mathcal{M}_{\omega l}^{T}&\mathbf{0}_{2\times
2n}&\mathcal{M}_{l}\\\ \end{array}\right]\dot{q},$

where $\mathcal{M}(\eta_{1},\eta_{2})\in\mathbb{R}^{(8+2n)\times(8+2n)}$ is
the overall system mass matrix, which is block partitioned.
$\mathcal{M}_{p}\in\mathbb{R}^{3\times 3}$,
$\mathcal{M}_{\omega}\in\mathbb{R}^{3\times 3}$,
$\mathcal{M}_{r}\in\mathbb{R}^{2n\times 2n}$, and
$\mathcal{M}_{l}\in\mathbb{R}^{2\times 2}$ are symmetric mass and inertia
matrices of the multirotor structure, the rotors, and the manipulator with
respect to $B$ frame. Matrices $\mathcal{M}_{p\omega}\in\mathbb{R}^{3\times
3}$, $\mathcal{M}_{pl}\in\mathbb{R}^{3\times 2}$, $\mathcal{M}_{\omega
r}\in\mathbb{R}^{3\times 2n}$, and $\mathcal{M}_{\omega
l}\in\mathbb{R}^{3\times 2}$ describe coupling effects. The total potential
energy of the AM system with respect to inertial frame is
$V_{AM}=ge_{3}^{T}\left((m_{b}+2nm_{r}){s}_{eb}+m_{1}s_{el_{1}}+m_{2}s_{el_{2}}\right).$
Notice $V_{AM}$ is a function of the direction of gravity, the position of the
vehicle, and manipulator joint angles, thus making the AM systems suitable to
apply theorem 1.

###  4.2 Gravitational Potential Forces and Moments

The previously defined advected parameter $\gamma$ and $\zeta$ satisfies the
following advection equations:

$\dot{\gamma}=-S(\omega_{b})\gamma\quad\quad\dot{\zeta}=-S(\omega_{b})\zeta+\dot{s}_{b}$
(14)

The Augmented Lagrangian, parametrized by $\gamma$ and $\zeta$, is
$SE(3)$-invariant by Theorem 1. The Augmented Lagrangian
$\overline{L}(q,\gamma,\zeta)$ for the system is:

$\begin{array}[]{ll}\overline{L}=\frac{1}{2}\dot{q}^{T}\mathcal{M}(q)\dot{q}-g\gamma^{T}(m_{t}\zeta-
m_{1}s_{bl_{1}}(\eta_{1})-m_{2}s_{bl_{2}}(\eta_{1},\eta_{2}))\end{array}$ (15)

where the total vehicle mass is lumped as $m_{t}=m_{b}+m_{1}+m_{2}+2nm_{r}$.
The conservative forces and momenta due to gravity can be derived from (4):

$\begin{array}[]{ll}-dV&=-\frac{\partial
V}{\partial\gamma}\dot{\gamma}dt-\frac{\partial
V}{\partial\zeta}\dot{\zeta}dt-\frac{\partial
V}{\partial\eta_{1}}\dot{\eta_{1}}dt-\frac{\partial
V}{\partial\eta_{2}}\dot{\eta_{2}}dt\\\ \end{array}$ (16)

Explicitly, the force and torque of gravity $\tau_{\mathbf{p}}$ and
$\tau_{\mathbf{l}}$ are:

$\displaystyle\tau_{\mathbf{p}}$ $\displaystyle=-gm_{t}\gamma^{T},$ (17)
$\displaystyle\tau_{\mathbf{l}}$
$\displaystyle=g\left(m_{1}s_{bl_{1}}^{T}+m_{2}s_{bl_{2}}^{T}\right)S(\gamma)$
(18) $\displaystyle-g\gamma^{T}\left(m_{1}\frac{\partial
s_{bl_{1}}}{\partial\eta_{1}}+m_{2}\left[\begin{array}[]{cc}\frac{\partial
s_{bl_{2}}}{\partial\eta_{1}}&\frac{\partial
s_{bl_{2}}}{\partial\eta_{2}}\end{array}\right]\mathbf{1}_{2\times
1}\right)e_{2}.$ (20)

###  4.3 Non-Conservative Forces and Moments

This section describes the non-conservative forces arising from aerodynamic
effects and manipulator actuation, and we only address aerodynamic effects
near hovering conditions. The thrust force produced by the $j^{th}$ rotor is
modeled as the product of a lumped thrust coefficient $c_{T}$ with the
rotating velocity squared. Letting
$\mathbf{\Omega}=\left[\begin{array}[]{ccc}\Omega_{1}&\cdots&\Omega_{2n}\end{array}\right]^{T}$,
be the rotor RPMs, and the total thrust is
$\tau_{f}=c_{T}\mathbf{\Omega}^{T}\mathbf{\Omega}$. The torques imparted on
the vehicle in the frame $B$ are:

$\displaystyle\left[\begin{array}[]{c}\tau_{r}\\\ \tau_{p}\\\
\tau_{y}\end{array}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cccc}0&rc_{T}&0&\cdots\\\
-rc_{T}&0&rc_{T}&\cdots\\\
-c_{Q}&c_{Q}&-c_{Q}&\cdots\end{array}\right]\left[\begin{array}[]{c}\Omega_{1}^{2}\\\
\vdots\\\ \Omega_{2n}^{2}\end{array}\right]\triangleq
W\left[\begin{array}[]{c}\Omega_{1}^{2}\\\ \vdots\\\
\Omega_{2n}^{2}\end{array}\right]$ (33)

where $c_{Q}$ is the lumped resistive torque coefficient that depends on rotor
geometry. Together, they form the non-conservative forces and moments vector
$\tau_{nc}\triangleq\left[\begin{array}[]{ccc}\tau_{r}&\tau_{p}&\tau_{y}\end{array}\right]^{T}$.
The manipulator joint actuators are modeled as a bilateral torque input at the
joints. Explicitly, the torque input at the revolute joint connecting aerial
base and link 1 is $\tau_{l_{1}}$ and connecting link 1 and link 2 is
$\tau_{l_{2}}$.

Variable RPM rotor thrusts are unilateral, making it incompatible with the
second requisite for Sussmann’s theorem. Alternatively, we adopt the time rate
of change of rotor RPMs as control inputs. It is physically reasonable to
assume that $\\{\dot{\Omega}_{1},\cdots,\dot{\Omega}_{n}\\}$ can be driven
bilaterally since the rotor dynamics can be modeled as
$\dot{\Omega}_{i}=\frac{1}{J}(\tau_{d}+\tau_{r})$, where $\tau_{d}>0$ is the
motor drive torque and $\tau_{r}<0$ is a passive resistive torque due to air
damping. For $\dot{\Omega}_{i}=0$, $\tau_{d}+\tau_{r}$ implies balance between
the drive and resistive torque, and $\tau_{r}<0$ when the rotor is spinning.
From this equilibrium position, one can increase $\tau_{d}$ to produce
positive $\dot{\Omega}_{i}$. Vice versa, if $\tau_{d}$ decreases yielding
$\tau_{d}+\tau_{r}<0$, a negative rate of change in RPM is possible,
$\dot{\Omega}_{i}<0$.

###  4.4 Reduced Dynamics:

In the standard basis for $\mathfrak{se}(3)$, the generalized rotation
momentum and translational momentum in $B$ frame are defined as
$\mathbf{l}\triangleq\frac{\partial\overline{l}}{\partial\omega_{b}}\in
T_{q}^{*}SO(3)$ and
$\mathbf{p}\triangleq\frac{\partial\overline{l}}{\partial\dot{s}_{b}}\in
T_{q}^{*}\mathbb{R}^{3}$:

$\left[\begin{array}[]{c}\mathbf{p}\\\
\mathbf{l}\end{array}\right]=\underbrace{\left[\begin{array}[]{cc}\mathcal{M}_{p}&\mathcal{M}_{p\omega}\\\
\mathcal{M}_{p\omega}^{T}&\mathcal{M}_{\omega}\end{array}\right]}_{\text{$\triangleq\mathcal{M}_{s}$}}\left[\begin{array}[]{c}\dot{s}_{b}\\\
\omega_{b}\end{array}\right]+\left[\begin{array}[]{cc}\mathcal{M}_{sr}&\mathcal{M}_{sl}\end{array}\right]\left[\begin{array}[]{c}\mbox{diag}(\mathbf{\Omega})\mathbf{\Omega}\\\
\dot{\eta}_{1}\\\ \dot{\eta}_{2}\end{array}\right].$

where $\mathcal{M}_{sr}\triangleq[0_{3\times 2n},\mathcal{M}_{\omega r}]^{T}$
and $\mathcal{M}_{sl}\triangleq[\mathcal{M}_{pl},\mathcal{M}_{\omega l}]^{T}.$
Since the AM system is constraint free, the connection is simply:

$\left[\begin{array}[]{c}\dot{s}_{b}\\\
\omega_{b}\end{array}\right]=\mathcal{M}_{s}^{-1}\left(\left[\begin{array}[]{c}\mathbf{p}\\\
\mathbf{l}\end{array}\right]-\left[\begin{array}[]{cc}\mathcal{M}_{sr}\mbox{diag}(\mathbf{\Omega})&\mathcal{M}_{sl}\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{\Omega}\\\ \dot{\eta}_{1}\\\
\dot{\eta}_{2}\end{array}\right]\right).$ (34)

Following Ostrowski (1998)’s work, the non-conservative/ dissipative forces
$\tau_{f}$ and torques (33) can be projected along the acting momentum
directions. Since the symmetry directions are a simple basis of
$\mathfrak{se}$(3), the momentum equation including conservative forces $-dV$
and nonconservative forces and moments $\tau_{nc}$ becomes:

$\dot{\varrho}_{i}=\left<\frac{\partial\overline{l}}{\partial\xi},[\xi,b_{i}]+\left<\frac{\partial
b_{i}}{\partial
r},\dot{r}\right>\right>+\left<-dV,b_{i}\right>+\left<\tau_{nc},b_{i}\right>$
(35)

Substituting $\tau_{f}$, (16), and (33), we can obtain the translational and
rotational momenta differential equations:

$\begin{array}[]{ll}\dot{\mathbf{p}}&=\mathbf{p}\times\omega_{b}+\tau_{\mathbf{p}}+c_{T}\mathbf{\Omega}^{T}\mathbf{\Omega}e_{3}\\\
\dot{\mathbf{l}}&=\mathbf{p}\times\dot{s}_{b}+\mathbf{l}\times\omega_{b}+\tau_{\mathbf{l}}+W\mbox{diag}(\mathbf{\Omega})\mathbf{\Omega}\end{array}$
(36)

Next, we derive the dynamics of the manipulator. Let
$\eta=[\eta_{1},\eta_{2}]^{T}$. With system Lagrangian (15), one can apply the
Euler-Lagrange equation for the shape variable $\eta_{i}$ to extract the
following the second-order ODE:

$\ddot{\eta}=\mathcal{M}_{\eta}^{-1}\left(-f_{\eta}+g_{\eta}\dot{\mathbf{\Omega}}+\left[\begin{array}[]{cc}\tau_{l_{1}}&\tau_{l_{2}}\end{array}\right]^{T}\right).$
(37)

where
$\mathcal{M}_{\eta}=\mathcal{M}_{l}+\mathcal{M}_{sl}^{T}\mathcal{M}_{s}^{-1}\mathcal{M}_{sl}$.
Further, $f_{\eta}$ and $g_{\eta}$ are:

$\displaystyle f_{\eta}$
$\displaystyle=\dot{\mathcal{M}}_{l}\dot{\eta}+\frac{d}{dt}\left(\left(\mathcal{M}_{sl}^{T}\mathcal{M}_{s}^{-1}\right)\left(\left[\begin{array}[]{c}\mathbf{p}\\\
\mathbf{l}\end{array}\right]-\left[\begin{array}[]{cc}\mathcal{M}_{sr}\mbox{diag}(\mathbf{\Omega})&\mathcal{M}_{sl}\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbf{\Omega}\\\
\dot{\eta}\end{array}\right]\right)\right)$
$\displaystyle+\frac{1}{2}\frac{\partial(\dot{q}^{T}\mathcal{M}(q))}{\partial\eta}\dot{q}-g\left(m_{1}\frac{\partial
s_{bl_{1}}}{\partial\eta}+m_{2}\frac{\partial
s_{bl_{2}}}{\partial\eta}\right)^{T}\gamma,$ $\displaystyle g_{\eta}$
$\displaystyle=-\mathcal{M}_{sl}^{T}\mathcal{M}_{s}^{-1}\mathcal{M}_{sr}\mbox{diag}(\mathbf{\Omega}).$

In summary, the reduced EOM for AM system are: the momenta equations (36), the
shape dynamics (37), and the advection equations (14). The Connection (34) can
be incorporated to reconstruct the system’s motion in $SE(3)$.

##  5 Nonlinear Controllability Assessment

To prepare for nonlinear controllability analysis, the EOM need to be
organized into control affine form. We defined the overall system state
variable to be:

$\mathbf{q}\triangleq[\mathbf{p}^{T},\mathbf{l}^{T},\phi,\theta,\psi,\eta_{1},\eta_{2},\dot{\eta}_{1},\dot{\eta}_{2},\Omega_{1},\cdots,\Omega_{2n}]^{T}\in\mathbb{R}^{2n+13}$

where we use the Euler angle instead of the advected parameter $\gamma$ to
track the the multirotor orientation. Further, the control affine form of the
full AM system is:

$\dot{\mathbf{q}}=f(\mathbf{q})+G(\mathbf{q})\mathbf{u},$ (38)

where the control input of the system is

$\mathbf{u}=[\tau_{l_{1}},\tau_{l_{2}},\dot{\Omega}_{1},\dot{\Omega}_{2},\cdots,\dot{\Omega}_{2n}]^{T}\in\mathbb{R}^{2n+2}.$
(39)

Further, the drift term and actuation vector fields are

$f(\mathbf{q})=\left[\begin{array}[]{c}\mathbf{p}\times\omega_{b}-m_{t}g\gamma+e_{3}c_{T}\mathbf{\Omega}^{T}\mathbf{\Omega}\\\
\mathbf{p}\times\dot{s}_{b}+\mathbf{l}\times\omega_{b}+\tau_{\mathbf{l}}+W\mbox{diag}(\mathbf{\Omega})\mathbf{\Omega}\\\
T(\mathbf{\Theta})^{-1}\omega_{b}\\\ \dot{\eta}\\\
-\mathcal{M}_{\eta}^{-1}f_{\eta}\\\ 0_{(2n)\times(2n+2})\end{array}\right],$
$G(\mathbf{q})=\left[\begin{array}[]{ccc}g_{1}(\mathbf{q})&\cdots&g_{2n+2}(\mathbf{q})\end{array}\right]=\left[\begin{array}[]{cc}0_{11\times
2}&0_{11\times 2n}\\\ \mathcal{M}_{\eta}^{-1}&g_{\eta}\\\ 0_{2n\times
2}&I_{2n\times 2n}\\\ \end{array}\right],$

respectively, where unknown states $\omega_{b}$ and $\dot{s}_{b}$ can be
expressed in terms of $\mathbf{q}$ using the connection (34) and
$\tau_{\mathbf{l}}$ is a function of $\gamma$ and $\eta$ given in (18) - (20).

###  5.1 Equilibrium Condition

Nonlinear controllability can only be assessed in the neighborhood of
equilibrium conditions. Setting $\dot{\mathbf{q}}=0$, one can verify the
following relationship must be true:

$\displaystyle 0_{1\times 9}$
$\displaystyle=\left[\begin{array}[]{ccccccccc}\mathbf{p}_{x}^{*}&\mathbf{p}_{y}^{*}&\mathbf{p}_{z}^{*}&\mathbf{l}_{x}^{*}&\mathbf{l}_{y}^{*}&\phi^{*}&\theta^{*}&\dot{\eta}_{1}^{*}&\dot{\eta}_{2}^{*}\end{array}\right],$
(41) $\displaystyle\mathbf{l}^{*}_{z}$
$\displaystyle=J(\mathbf{\Omega}^{*})^{T}\mathbf{\Omega}^{*},\quad\quad
c_{T}(\mathbf{\Omega}^{*})^{T}\mathbf{\Omega}^{*}=m_{t}g,$ (42) $\displaystyle
W\mbox{diag}(\mathbf{\Omega}^{*})\mathbf{\Omega}^{*}$
$\displaystyle=-ge_{3}\times\left(m_{1}s_{bl_{1}}(\eta_{1})+m_{2}s_{bl_{2}}(\eta_{1},\eta_{2})\right),$
(43) $\displaystyle+ge_{3}^{T}\left(m_{1}\frac{\partial
s_{bl_{1}}}{\partial\eta_{1}}+m_{2}\left[\begin{array}[]{cc}\frac{\partial
s_{bl_{2}}}{\partial\eta_{1}}&\frac{\partial
s_{bl_{2}}}{\partial\eta_{2}}\end{array}\right]\mathbf{1}_{2\times
1}\right)e_{2}.$ (45)

Further, $\eta_{1}$, $\eta_{2}$, and $\psi$ are free, meaning any combinations
of the free states can be an equilibrium condition of the overall system.
Physically speaking, these conditions correspond to a hovering AM at arbitrary
vehicle position and manipulator configurations with zero roll and pitch.

###  5.2 Small-Time Locally Accessible and Controllability

###### Theorem 4

All AMs are Small-Time Locally Accessible (STLA) evaluated at the equilibrium
conditions (41) - (43), except at manipulator singularity
$\eta=[\pi/2,0]^{T}$. By Theorem 2, we must check that system (38) satisfied
the LARC when evaluated at the equilibrium condition, equivalently
$\overline{\mbox{Lie}}({f,g_{1}...,g_{m}})(\mathbf{q}^{*})=T_{\mathbf{q}^{*}}Q$.

First note that $g_{1}(\mathbf{q}^{*})$ and $g_{2}(\mathbf{q}^{*})$ are
functions of the free state $\eta_{2}$, and $\forall\eta_{2}$,
$\mbox{rank}([g_{1}(\mathbf{q}^{*}),g_{2}(\mathbf{q}^{*})])=2$, spanning the
states corresponds to $\ddot{\eta}_{1}$ and $\ddot{\eta}_{2}$. Despite
coupling with the $\ddot{\eta}_{1}$ and $\ddot{\eta}_{2}$, the remaining first
degree actuation fields have $\mbox{rank}(g_{3},\cdots,g_{2n+2})=2n$,
annihilating any nonzero $\dot{\mathbf{\Omega}}$ with appropriate combinations
of $\mathbf{u}$. Since $g_{1}$ and $g_{2}$ corresponds to manipulator
maneuvers in the $x_{b}-z_{b}$ plane, the Lie brackets between $g_{1}$,
$g_{2}$ with the drift $f$ are:

$\begin{array}[]{ll}\left[f,g_{i}\right]\big{|}_{\mathbf{q}^{*}}=\left[\begin{array}[]{ccccccccc}0^{\dot{\mathbf{p}}}_{1\times
3}&f^{[\dot{\mathbf{l}}_{x},\dot{\mathbf{l}}_{y}]}_{1\times
2}&0^{\dot{\mathbf{l}}_{z}}&f^{\dot{\mathbf{\Theta}}}_{1\times
3}&f^{[\dot{\eta},\ddot{\eta}]}_{1\times 4}&0^{\mathbf{\dot{\Omega}}}_{1\times
2n}\end{array}\right]^{T},\end{array}$

which are two independent actuation fields that primarily influences
$\dot{\eta}_{1}$ and $\dot{\eta}_{2}$ and also perturbs the angular velocities
and rotational momentum because of the shifting in center of gravity. The
remaining $2n$ actuation vector fields together span $\dot{\mathbf{l}}$ and
$\dot{\mathbf{p}}_{z}$. For $i\in\\{3,2n+2\\}$, we have

$\left[f,g_{i}\right]\big{|}_{\mathbf{q}^{*}}=\left[\begin{array}[]{ccccccccc}0^{[\dot{\mathbf{p}}_{x},\dot{\mathbf{p}}_{y}]}_{1\times
2}&f^{\dot{\mathbf{p}}_{z}}&f^{\dot{\mathbf{l}}}_{1\times
3}&f^{\dot{\mathbf{\Theta}}}_{1\times 3}&f^{[\dot{\eta},\ddot{\eta}]}_{1\times
4}&0^{\mathbf{\dot{\Omega}}}_{1\times 2n}\end{array}\right]^{T},$

where $f^{\dot{\mathbf{\Theta}}}_{1\times 3}$ is non-zero because of the
unbalance in rotor thrust leads to a non-zero resistive torque resulting in
non-zero body rates. The manipulator joints are also coupled with the rotor
thrust rate via $g_{\eta}$ in (37). Therefore, any perturbation in rotor
thrust will leads to a non-zero rotation. Aided by the asymmetry, all degree
two brackets together add six more independent control distributions. Looking
at the third degree brackets and for $i\in\\{1,2\\}$,

$\left[[f,g_{i}],f\right]\big{|}_{\mathbf{q}^{*}}=\left[\begin{array}[]{ccccccc}f^{[\dot{\mathbf{p}}_{x},\dot{\mathbf{p}}_{y}]}_{1\times
2}&0^{\dot{\mathbf{p}}_{z}}&f^{\dot{\mathbf{l}}}_{1\times
3}&f^{\dot{\mathbf{\Theta}}}_{1\times 3}&f^{[\dot{\eta},\ddot{\eta}]}_{1\times
4}&0^{\mathbf{\dot{\Omega}}}_{1\times 2n}\end{array}\right]^{T},$

which provides two more independent bases for $T_{\mathbf{q^{*}}}Q$.
Physically speaking, the third degree brackets with $g_{1}$ and $g_{2}$ are
similar to the second degree brackets where shifting the vehicle center of
gravity perturbs the hovering motion. It is important to note that the
perturbation will come in the form of forces which leads to nonzero
translational accelerations. Further, the remaining third order brackets
provided 3 more independent bases. For $i\in\\{3,2n+2\\}$,

$\left[[f,g_{i}],f\right]\big{|}_{\mathbf{q}^{*}}=\left[\begin{array}[]{ccccccccc}f^{[\dot{\mathbf{p}}_{x},\dot{\mathbf{p}}_{y}]}&0^{\dot{\mathbf{p}}_{z}}&f^{\dot{\mathbf{l}}}_{1\times
3}&f^{\dot{\mathbf{\Theta}}}_{1\times
3}&&f^{[\dot{\eta},\ddot{\eta}]}_{1\times
4}&0^{[\mathbf{\dot{\Omega}}]}_{1\times 2n}\end{array}\right]^{T}.$

Together, there are $13+2n$ good brackets that can be used in spanning the
tangent space $T_{\mathbf{q}^{*}Q}$ except at the manipulator singularity
($\eta=[\pi/2,0]$). Therefore, by theorem 2, we conclude AM are STLA evaluated
at the equilibrium condition (41)-(43) except at system singularity.
$\blacksquare$

###### Theorem 5

All AMs are Small-Time Locally Controllable (STLC) evaluated at the
equilibrium conditions (41) - (43), except at manipulator singularity
$\eta=[\pi/2,0]^{T}$. By definition of good and bad brackets, all Lie brackets
used to establish STLA are good brackets, and the highest degree bracket used
in spanning $T_{\mathbf{q}^{*}}Q$ is three. To apply theorem 3, we need to
show all bad brackets of degree $j\leq 3$ can be expressed as a linear
combination of good brackets of degree less than $j$. By its definition, there
will be no bad brackets associated with an even degree. The third degree bad
brackets can be expressed generally as $[[f,g_{i}],g_{i}],\forall
i\in\\{1,\cdots,2n+2\\}$. More specifically, for $i\in\\{1,2\\}$ the bad
brackets evaluated at the equilibrium condition takes the following form:

$[[f,g_{i}],g_{i}]\big{|}_{\mathbf{q}^{*}}=\left[\begin{array}[]{cccc}0^{[\dot{\mathbf{p}},\dot{\mathbf{l}},\dot{\mathbf{\Theta}}]}_{1\times
12}&0^{\dot{\eta}}_{1\times 2}&f^{\ddot{\eta}}_{1\times
2}&0^{\mathbf{\dot{\Omega}}}_{1\times 2n}\end{array}\right],$

which can be annihilated with good, degree one bracket $g_{1}$ and $g_{2}$.
The remaining $2n$ third degree bad brackets take the following form where
$i\in\\{3,\cdots,2n\\}$:

$\left[\begin{array}[]{cccccc}[[f,g_{i}],g_{i}]\big{|}_{\mathbf{q}^{*}}=0^{[\dot{\mathbf{p}}_{x},\dot{\mathbf{p}}_{y}]}_{1\times
2}&f^{\dot{\mathbf{p}}_{z}}&f^{\dot{\mathbf{l}}}_{3\times
1}&f^{\dot{\mathbf{\Theta}}}_{3\times 1}&f^{[\dot{\eta},\ddot{\eta}]}_{4\times
1}&0^{\dot{\mathbf{\Omega}}}_{2n\times 1}\end{array}\right],$

which can be annihilated by good degree two brackets $[f,g_{i}]$ and degree
one brackets $g_{1}$ and $g_{2}$. The final bad bracket is the drift term
evaluated at equilibrium conditions. Since by the definition of equilibrium
condition, we have
$\dot{\mathbf{q}}^{*}=f(\mathbf{q}^{*})+G(\mathbf{q}^{*})\mathbf{u^{*}}=0$.
Thus, all bad brackets of degree $k\leq 3$ can be written as combinations of
good Lie brackets of degree $j\leq k$, establishing STLC.$\blacksquare$

It is important to note that when interpreting these results, the net motions
generate from higher degree Lie brackets ($d>1$) are “slower” than motions
driven directly by actuation vectors fields $g_{1},\cdots,g_{m}$. In fact, the
net motions are $O(\epsilon^{d})$ for time $O(\epsilon)$, where $\epsilon\ll
1$. See Choset et al. (2005). Aided by net external motions resulting from the
coupling of internal shape changes, the AM system can achieve local
accessibility and controllability. However, we suggest the more valuable
takeaway from our SLTA and SLTC analysis should be the compositions of good
and bad lie brackets since they contain information about how well the system
can maneuver.

In summary, Lagrangian reduction process, directly produces the minimum set of
first-order dynamical equations, substantially simplifies AM systems’
controllability analysis. Using the connection, the “reduced” EOM can be
reformulated into the control-affine form effortlessly. In comparison, to
assess controllability without the reduced dynamics, using the EOM given by
Kim et al. (2013) requires a costly symbolic inversion of a
$(12+2n)\times(12+2n)$ sized mass matrix to reformulate. Further, the
reduction and reconstruction process can be generalized to analyze other
flying vehicles that include conventional quadcopter, tilt-rotors, multi-
rotors with asymmetric rotor placements, tethers, and pendulum appendages.
Concurrently, we are extending the geometric reduction and controllability
analysis to more general flying vehicles.

##  6 Conclusion

This paper analyzed two aspects of a class of aerial-manipulator robots.
First, aerial manipulator robots’ EOM was developed using Lagrangian reduction
and reconstruction despite having broken symmetry. Second, the AM system
properties, STLA and STLC, are formally analyzed. We have concluded that
underactuated multi-rotor with planar two-link manipulators are STLA and STLC
near equilibrium condition. A physical interpretation of this result is if the
AM is hovering, there exists control actions to maintain hover while the
manipulator joints track a non-singular and smooth trajectory. However, our
STLC result should be applied with caution since a controllability certificate
does not assess the relative control effort needed to realize the trajectory
or hovering state.

## References

  * Bloch et al. (1996) Bloch, A.M., Krishnaprasad, P.S., Marsden, J.E., and Murray, R.M. (1996).  Nonholonomic Mechanical Systems with Symmetry.  _Archive for Rational Mechanics and Analysis_ , 136(1), 21–99.
  * Bullo and Lewis (2004) Bullo, F. and Lewis, A. (2004).  _Geometric Control of Mechanical Systems. Modeling, Analysis, and Design for Simple Mechanical Control Systems_.
  * Burkhardt and Burdick (2016) Burkhardt, M. and Burdick, J.W. (2016).  Reduced dynamical equations for barycentric spherical robots.  In _IEEE Int. Conf. Robot. and Autom._ , 2725–2732.
  * Burkhardt (2018) Burkhardt, M.R. (2018).  _Dynamic Modeling and Control of Spherical Robots_.  Ph.D. thesis, Caltech.
  * Choset et al. (2005) Choset, H., Lynch, K.M., Hutchinson, S., Kantor, G.A., Burgard, W., Kavraki, L.E., and Thrun, S. (2005).  _Principles of Robot Motion: Theory, Algorithms, and Implementations_.
  * Hermann (1963) Hermann, R. (1963).  On the accessibility problem in control theory.  In _Int. Symp. Nonlinear Diffe. Eqns. and Nonlinear Mechanics_ , 325–332. Academic Press.
  * Holm et al. (1998) Holm, D.D., Marsden, J.E., and Ratiu, T.S. (1998).  The euler–poincaré equations and semidirect products with applications to continuum theories.  _Advances in Mathematics_ , 137(1), 1–81.
  * Kim et al. (2013) Kim, S., Choi, S., and Kim, H.J. (2013).  Aerial manipulation using a quadrotor with a two dof robotic arm.  In _IEEE/RSJ Int. Conf. Intel. Robots Systems_ , 4990–4995.
  * Murray et al. (1994) Murray, R.M., Li, Z., and Sastry, S. (1994).  _A mathematical introduction to robotic manipulation_.
  * Ostrowski (1996) Ostrowski, J.P. (1996).  _The mechanics and control of undulatory robotic locomotion_.  Ph.D. thesis, Caltech.
  * Ostrowski (1998) Ostrowski, J. (1998).  Reduced equations for nonholonomic mechanical systems with dissipative forces.  _Reports on Mathematical Physics_ , 42(1), 185–209.
  * Schneider (2002) Schneider, D. (2002).  Nonholonomic euler-poincarè equations and stability in chaplygin’s sphere.  _Dynamical System_ , 17.
  * Schwarz (1975) Schwarz, G.W. (1975).  Smooth functions invariant under the action of a compact lie group.  _Topology_ , 14, 63–68.
  * Sussmann (1987) Sussmann, H. (1987).  A general theorem on local controllability.  _SIAM J. Control and Optimization_ , 25, 158–194.


# 6G and Beyond: Joint Waveform and Wavefront Engineering To Realize Terahertz
Communications

Arjun Singh, Vitaly Petrov, Hichem Guerboukha, Innem V.A.K. Reddy,
Edward W. Knightly, Daniel M. Mittleman, and Josep M. Jornet Arjun Singh is
with the SUNY Polytechnic Institute, USA. Vitaly Petrov and Josep M. Jornet
are with Northeastern University, USA. I.V.A.K. Reddy is with the University
at Buffalo, USA. Hichem Guerboukha and Daniel M. Mittleman are with Brown
University, USA. Edward W. Knightly is with Rice University, USA.The work has
been supported in part by the National Science Foundation CNS-1955004,
1954780, 2148132, 2211616, 1955075, 1923782, 1824529, and 2148132, and the US
Air Force Office of Scientific Research Award FA9550-23-1-0254 and
FA9550-22-1-0412, and the Army Research Laboratory grant W911NF-19-2-0269. H.
Guerboukha also gratefully acknowledges support from Fonds de recherche du
Quebec – Nature et technologies.

# 6G and Beyond: Wavefront Engineering to Realize Efficient Terahertz Band
Communications in 6G and Beyond

Arjun Singh, Vitaly Petrov, Hichem Guerboukha, Innem V.A.K. Reddy,
Edward W. Knightly, Daniel M. Mittleman, and Josep M. Jornet Arjun Singh is
with the SUNY Polytechnic Institute, USA. Vitaly Petrov and Josep M. Jornet
are with Northeastern University, USA. I.V.A.K. Reddy is with the University
at Buffalo, USA. Hichem Guerboukha and Daniel M. Mittleman are with Brown
University, USA. Edward W. Knightly is with Rice University, USA.The work has
been supported in part by the National Science Foundation CNS-1955004,
1954780, 2148132, 2211616, 1955075, 1923782, 1824529, and 2148132, and the US
Air Force Office of Scientific Research Award FA9550-23-1-0254 and
FA9550-22-1-0412, and the Army Research Laboratory grant W911NF-19-2-0269. H.
Guerboukha also gratefully acknowledges support from Fonds de recherche du
Quebec – Nature et technologies.

# Wavefront Engineering to Realize Efficient Terahertz Band Communications in
6G and Beyond

Arjun Singh, Vitaly Petrov, Hichem Guerboukha, Innem V.A.K. Reddy,
Edward W. Knightly, Daniel M. Mittleman, and Josep M. Jornet Arjun Singh is
with the SUNY Polytechnic Institute, USA. Vitaly Petrov and Josep M. Jornet
are with Northeastern University, USA. I.V.A.K. Reddy is with the University
at Buffalo, USA. Hichem Guerboukha and Daniel M. Mittleman are with Brown
University, USA. Edward W. Knightly is with Rice University, USA.The work has
been supported in part by the National Science Foundation CNS-1955004,
1954780, 2148132, 2211616, 1955075, 1923782, 1824529, and 2148132, and the US
Air Force Office of Scientific Research Award FA9550-23-1-0254 and
FA9550-22-1-0412, and the Army Research Laboratory grant W911NF-19-2-0269. H.
Guerboukha also gratefully acknowledges support from Fonds de recherche du
Quebec – Nature et technologies.

# Wavefront Engineering: Realizing Efficient Terahertz Band Communications in
6G and Beyond

Arjun Singh, Vitaly Petrov, Hichem Guerboukha, Innem V.A.K. Reddy,
Edward W. Knightly, Daniel M. Mittleman, and Josep M. Jornet Arjun Singh is
with the SUNY Polytechnic Institute, USA. Vitaly Petrov and Josep M. Jornet
are with Northeastern University, USA. I.V.A.K. Reddy is with the University
at Buffalo, USA. Hichem Guerboukha and Daniel M. Mittleman are with Brown
University, USA. Edward W. Knightly is with Rice University, USA.The work has
been supported in part by the National Science Foundation CNS-1955004,
1954780, 2148132, 2211616, 1955075, 1923782, 1824529, and 2148132, and the US
Air Force Office of Scientific Research Award FA9550-23-1-0254 and
FA9550-22-1-0412, and the Army Research Laboratory grant W911NF-19-2-0269. H.
Guerboukha also gratefully acknowledges support from Fonds de recherche du
Quebec – Nature et technologies.

###### Abstract

Terahertz (THz) band communications is envisioned as a key technology for
future wireless standards. Substantial progress has been made in this field,
with advances in hardware design, channel models, and signal processing. High-
rate backhaul links operating at sub-THz frequencies have been experimentally
demonstrated. However, there are inherent challenges in making the next great
leap for adopting the THz band in widespread communication systems, such as
cellular access and wireless local area networks. Primarily, such systems have
to be both: (i) wideband, to maintain desired data rate and sensing
resolution; and, more importantly, (ii) operate in the massive near field of
the high-gain devices required to overcome the propagation losses. In this
article, it is first explained why the state-of-the-art techniques from lower
frequencies, including millimeter-wave, cannot be simply repurposed to realize
THz band communication systems. Then, a vision of wavefront engineering is
presented to address these shortfalls. Further, it is illustrated how novel
implementations of specific wavefronts, such as Bessel beams and Airy beams,
offer attractive advantages in creating THz links over state-of-the-art far-
field beamforming and near-field beamfocusing techniques. The paper ends by
discussing novel problems and challenges in this new and exciting research
area.

###### Index Terms:

Terahertz Communications; 6G; Wavefront Engineering; Bessel beams; Near field;
Orbital Angular Momentum

##  I Introduction

The terahertz (THz) band (0.3–3 THz) has been heavily investigated by the
scientific community over the past decades for wireless sensing and
communications. With developments relating to device technology, new waveform
modulations, signal processing techniques, and robust digital back-ends, the
so-called “THz technology gap” is slowly closing [1]. However, the open
question that remains is how to design a practical THz band communication
system that delivers on the promises of the next generations of wireless
networks, which are tentatively marked as 1 terabit-per-second (Tbps) data
rate, 0.1 ms latency, and “$10^{-9}$” reliability [1].

Here, Shannon’s capacity gives us major metrics to work with. More
specifically, the channel capacity depends on: (i) the available bandwidth,
(ii) the received signal-to-noise ratio (SNR), and (iii) the spatial reuse
factor. While the THz band provides massive available bandwidth for ultra-
broadband modulations, the corresponding noise power, crippling path losses,
and the fact that THz signals are easily obstructed by everyday objects can
have devastating consequences on the SNR, leading to link breakage. The THz
channel is also low rank, making many un-correlated multiple-input-multiple-
output (MIMO) links unlikely [2]. These problems have thus spurred significant
research efforts [3, 1, 4].

To this end, in addition to the utilization of very high gain antenna arrays
[5], intelligent reflecting surfaces (IRSs) are expected to support THz links
[3]. The main goal with these devices is to utilize beamforming theory to
increase the radiation gain, spatial multiplexing capability, and robustness
of reconfigurable links including non-line-of-sight (NLoS) capability around
obstacles. The very small wavelengths at THz frequencies result in such high-
gain devices being relatively compact but having a massive near-field region.
Thus, in most practical setups, the conventional far-field beamforming
assumptions of the microwave, and to a large extent, the millimeter-wave
(mmWave) bands, do not hold [6, 5].

Here, we observe that the properties of a beam generated from an aperture are
completely defined by the phase and, in some cases, amplitude distribution of
the electric field at the aperture itself. Thus, the knowledge of the
wavefront, i.e., the imaginary line that connects all the points of a wave
with the same phase, is enough to completely characterize the beam [7].
However, exploring the concept of wavefront engineering at THz frequencies to
generate beams that have well-defined propagation characteristics in the near-
field is a significant task, much more complicated than canonical phased array
beamforming. Aside from _near-field beamfocusing_ which utilizes a lens-like
approach from the optical domain, other wavefronts and corresponding beams
have also been explored, primarily for imaging and sensing purposes [8, 9].
The THz band however provides a Goldilocks zone for a novel implementation of
such wavefronts in realizing efficient communication links as well.

In this paper, we present _a vision for THz wavefront engineering to realize
practical THz band communications_ by utilizing wavefronts that can operate in
the near field, mitigate the blockage problem, and guarantee multiple un-
correlated channels matching perfect MIMO expectations, without the need of
advanced signal processing algorithms or active relays. The present tutorial
should boost further research in the growing area of combining
electromagnetics and communications for physical layer design in 6G and beyond
wireless systems.

Figure 1: Near-field and far-field 6G-grade sub-THz prospective use cases (300
GHz carrier and a 10-15 cm aperture at the larger unit or the base station)
and beyond-6G true THz setups (1 THz carrier and a 25-cm aperture at the base
station for increased coverage).

##  II The Inherent Problem of THz Communications

Modern radio communications have been, for the most part, explored and
developed under two major postulations: (i) the far-field assumption; and (ii)
the narrow-band assumption. _Prospective ultra-high data rate THz systems
contradict both these assumptions._

###  II-A Near field vs. Far field

The small wavelengths of THz signals result in small radiating elements. Thus,
even moderately sized physical devices (up to a few cm) can pack many
radiating elements, providing a high gain for overcoming large path losses. At
the same time, this results in large electrical apertures (ratio of physical
size to signal wavelength). The relative phase delays and signal strengths of
the radiation from the elements of such apertures are configured such that
their electromagnetic (EM) superposition creates the desired beam profile,
such as concentrating the power in the desired spot/direction.

Under the far-field assumption, the resultant superimposed beam has a uniform
phase across its cross-section, effectively like a plane wave. However, as the
near-field region of an array increases exponentially with its physical size
and linearly with the frequency, the far field Fraunhofer distance for
practical THz communication links can easily be _hundreds of meters away_ from
the transmitter [5]. Hence, many prospective THz use cases will operate in the
near field, as illustrated in Fig. 1.

###  II-B Wideband and MIMO issues

The second fundamental assumption that does not always hold at THz frequencies
is the narrow-band assumption. Following it, the bandwidth of the signal is
considered negligibly small compared to the central frequency; the entire
signal can be approximated as a single “wave" at this design frequency.

However, one of the key motivations in exploring new spectrum is to achieve
extremely high data rates and sensing resolution by harnessing tens or even
hundreds of gigahertz (GHz) of the spectrum. Thus, the narrow-band assumption
may not hold, and the same phase delay across the array elements no longer
provides a uniform time delay across the signal bandwidth when moving away
from a central frequency. This effect is also present in mmWave communications
but becomes more profound in wider-band THz systems [1].

Exploitation of MIMO techniques such as spatial multiplexing is one way to
address the wideband-related difficulties without compromising on the desired
capacity. However, the sub 6 GHz MIMO point-to-point spatial multiplexing
methods that rely on rich multipath may not be viable in THz channels.
Specifically, preliminary THz MIMO studies indicate that only a few dominant
paths of communication exist in the THz channel, and are strongly correlated,
thus making the THz channel strongly Rician, low rank, and sparse [2]. Further
still, communication devices would still require an advanced THz MIMO
architecture with powerful back-end processing for processing the independent
streams accurately, naturally increasing the system complexity and costs even
further [1].

###  II-C Blockage of the THz signal

Another inherent limitation of THz propagation is that THz signals are easily
blocked by many objects in the environment, including building walls,
furniture, and even parts of the human body. Today, the default approaches to
mitigate blockage obstructions are: (i) applying IRSs to extend the coverage
and engineer NLoS paths [3]; and (ii) exploiting THz multi-connectivity
solutions [1]. Both approaches are feasible but come with an increased system
design complexity.

##  III Wavefront Engineering
to Conquer The THz Channel

TABLE I: Comparison between different wavefront candidates for THz communications Type |  Phase profile | Beam profile |  Self-healing |  OAM |  Self-accelerating
---|---|---|---|---|---
|  |  Near field |  Far field |  |  |
Beamforming |  Planar |  Undefined |  Gaussian |  X |  ✓ |  X
Beamfocusing |  Quadratic |  Point |  Undefined |  X |  X |  X
Bessel beams |  Conical |  Bessel |  Annular |  ✓ |  ✓ |  X
Airy beams |  Exponential |  Airy |  Undefined |  ✓ |  ✓ |  ✓

Here, we introduce the vision of utilizing novel wavefronts to realize robust
THz communications. The discussed THz wavefronts are illustrated in Fig. 2 and
their key features are summarized in Table I.

###  III-A Wavefront Generation

The Huygens-Fresnel principle shows that the characteristics of a generated
beam are completely defined by the phase and amplitude distribution of the
electric field at the radiating aperture [7]. Therefore, to generate any of
the beams, all that is required is the corresponding wavefront which can be
engineered by a specific phase profile. In addition to implementing this
through an array, the same is also possible through custom-designed lenses, or
in reflection through the utilization of reflectarrays and metasurfaces. A
more in-depth discussion of these different design implementations is
presented in Sec V.

###  III-B Beamfocusing and beamforming

Beamfocusing and beamforming are two canonical, interrelated wave types in the
near field and far field, respectively. With _beamfocusing_ , a transmitting
array is configured with a lens-like quadratic phase profile to focus the
signal at a particular focal point, at which point the wavefront converges to
a singularity [6]. Beam intensity is maximum at this focal point, and the Abbe
limit, which dictates the resolution capability of the array, gives the
resultant beam spot size.

_Beamforming_ is the far-field extension of beamfocusing, such that the focal
point is, effectively, moved to infinity – the quadratic phase becomes a
constant, uniform phase and the wavefront becomes planar, resulting in a beam
that diverges in the far-field due to limited aperture size. Both beams have
their limitations. Beamfocusing is focused at a given point, thus requiring
perfect CSI and mobility tracking, while beamforming is valid only in the far
field. Neither of these has any strong resilience to blockage.

###  III-C Bessel beams

When discussing the candidate solutions for near-field THz systems, one of the
solutions constitutes a beam intensity profile given by a zeroth-order Bessel
function, giving these their name – _“Bessel beams"_ [8]. A radially symmetric
linear phase creates this Beam profile, resulting in a conical wavefront.

As illustrated in Fig. 2, the zero-order Bessel beam profile has a central
bright spot along the cone axis, with multiple concentric rings around it,
which are interference patterns created from the interference of plane waves
from the opposite sides of this central axis. Bessel beams are exact solutions
to Maxwell’s equations with a constant beam intensity. Thus, as the wave
propagates, the beam profile remains invariant (no spreading losses). The
phase and total power within the central spot and all the concentric rings are
the same, with much greater intensity in the central spot. However, for a
finite aperture and power, the propagation distance over which a practical
Bessel beam retains its properties is limited. Some of the major design
considerations for Bessel beams are that, for a given aperture size, the angle
of the conical wavefront directly affects the resultant diameter of the
central spot, the number of rings, and the maximum propagation distance.
Further, the electrical size of the aperture also factors into these
considerations and all can be customized. For example, a particular design
constraint is that Bessel beams generated from the same aperture but with a
larger number of rings theoretically lead to better self-healing capabilities,
while a smaller number of rings and a wider central bright spot leads to
greater propagation distance [8].

Bessel beams propagate in the near field, a feature absent in both beamforming
(not valid in the near field) and beamfocusing (does not propagate but
converges at a focal point). In addition, as Bessel beams can be understood to
be the interference pattern of plane waves traveling inwards on a cone, they
are resilient to blockage; even if some of the waves are blocked by an
obstruction, the remaining waves still reconstruct the interference pattern
and help regenerate the Bessel beam after the obstruction; hence these are
also referred to as _self-healing beams_ [10]. These properties make Bessel
beams a promising wavefront candidate for THz communications.

Figure 2: Wavefront engineering for creating wireless THz links, in which the
beam profiles can significantly deviate from the Gaussian.

###  III-D Curved and Airy beams

While Bessel beams can “plow” through an obstacle, there is another class of
beams that can completely circumvent blockage by following a curved
trajectory. When viewed in the transverse direction, the beams appear to have
acceleration without any external energy – giving them the name of _self-
accelerating beams_. Such beams have been recently demonstrated in EM
applications [9], where the beam profile can be described with the Airy
function. The resultant Airy-like beam is a subset of self-accelerating beams,
with the cross-sectional intensity illustrated in Fig. 2.

Using the concept of caustics from ray optics, an arbitrary curve can be
generated from a large input aperture where the local required phase, and
subsequent wavefront, are found using tangents of the curve extended to the
aperture. A particular design consideration of these beams is that greater
beam curvature requires both a steeper phase progression across the array
aperture as well as a larger aperture size. In THz wireless communications,
their ability to follow almost any curved, or convex, path makes them very
promising to avoid signal blockage. In addition, accelerating beams also have
the potential to change the angle of arrival at the receiver side, making them
useful during beam alignment.

###  III-E OAM – A New Feature

Certain wavefronts possess new properties which open the door to new
communication strategies. The most widely researched and explored among these
is orbital angular momentum (OAM), a property of light first demonstrated in
[11]. A beam that has OAM manifests a spiral phase in the transverse
direction, resulting in a helical wavefront with a zero-intensity vortex in
the center, as shown in Fig. 2.

Beams with distinct OAM modes are always orthogonal to each other, and the
resulting cross-correlation is guaranteed to be zero. Therefore, OAM mode
multiplexing can be utilized to create simultaneous uncorrelated data streams
[4]. Gaussian, Bessel, and Airy beams can all be customized to carry OAM by
adding a spiral phase on top of the existing beam profile. Thus, in the near
field, Bessel beams with OAM further augment the attractive features of Bessel
beams [10].

###  III-F Receiver Performance and Operation

Employing wavefront engineering at the transmitter and exploiting the
properties of the discussed beams can truncate some of the more significant
drawbacks of the THz channel, such as the near-field effect, blockage, and
wideband-related difficulties. Ultimately, the objective is to provide greater
signal power available at the receiver, thereby improving the SNR and
performance metrics.

Here, we note that during near-field operation, both the magnitude of the
signal reaching the receiver and its coherent reception at the receiver
becomes important. For example, with beamfocusing, if the receiver is not at
the focal spot but slightly closer to the transmitter, the incident power is
spread over an arc. Hence, the phase is non-uniform across the receiver cross-
section. The receiver must then be able to accommodate the required delay and
sum method for coherently adding such a non-uniform power profile [6].

With THz Bessel beams, the entire Bessel intensity pattern has the same phase,
thereby making the receiver operation similar to that within far-field
beamforming [7]. With THz Airy beams, the receiver should be tuned to be
coherent towards the incident angle defined by the curve of the Airy beam
(similar to the far-field beamsteering direction of arrival).

When required to demultiplex multiple OAM modes, the receiver must employ a
large enough aperture with a specific phase-based diffraction grating that
exploits the principle of zero cross-correlation in OAM modes for ease of
demultiplexing [12]. Nonetheless, the fundamental principles remain the same;
more power received coherently allows for greater SNR and performance metrics.

##  IV THz Wavefronts: A Benchmark Study

In this section, we present a numerical study illustrating how the described
wavefronts and beams can quantitatively assist in mitigating the challenges of
THz band communications. We consider a planar array of a side length of 25 cm
operating at 1 THz and designed with a typical half-wavelength spacing of
elements. The nominal gain of such an array is 70 dBi, which can help sustain
local area, multi-GHz wide links.

###  IV-A Improved Radiation Gain

Figure 3: Improved radiation gain with wavefront engineering.

The effective normalized radiation gain is a key measure of how accurately an
applied phase profile at the array constructively superimposes the waves from
each of the array elements at a given evaluation point, and linearly affects
the derived directivity of the array, with a maximum value of $1$. Both the
received power and the SNR at the receiver are linearly proportional to the
effective normalized radiation gain.

In Fig. 3, we present this normalized radiation gain under the application of
Bessel beams, beamfocusing, and beamforming, all designed to steer in the
broadside direction of the array. The peculiar beam profile of Airy beams
cannot be grouped with the other beam shapes, since the Airy beam follows a
curved trajectory. Nonetheless, for the path specified during wave generation,
the Airy beam will have an intensity maximum along that curve. The diameter of
the central spot of the Bessel is set to 20 mm, which, given the aperture
size, results in a maximum propagation distance of 23 m.

It is observed from Fig. 3 that under beamforming the normalized gain is
notably crippled at distances of up to 20 m, which are typical ranges for
prospective indoor THz WLANs. In contrast, the near field validity of Bessel
beams results in a significantly improved radiation gain for the same
settings. Arbitrarily, the focal point of beamfocusing is set to the distance
where the normalized gain with the Bessel beam is maximized (at 13.5 m). The
corresponding curves intersect at that point. However, even small deviations
from this focal point (i.e., due to mobility) result in the efficiency of
beamfocusing reducing rapidly and much faster than the Bessel curve. Hence, in
the absence of a perfectly updating CSI (impractical in many setups),
sustaining acceptable performance for mobile THz systems with beamfocusing is
unlikely. The improvement in beamforming will be slightly faster for smaller
apertures or lower frequencies since, in both cases, the electrical size of
the aperture decreases. However, these changes simultaneously lead to a
reduced aperture gain that challenges coverage and capacity. The qualitative
system behavior remains the same: _THz Bessel beams are more efficient for
shorter distances (up to several tens of meters), while THz beamforming is
more efficient for longer links closer to the far-field region._

###  IV-B Blockage Mitigation

Figure 4: Blockage mitigation with different wavefronts. From top to bottom:
beamforming, beamfocusing, Bessel, Airy-like curved beam. Full wave finite-
element-method (FEM) with COMSOL Multiphysics is utilized.

The potential of wavefront engineering in mitigating the impact of obstacles
to THz signal propagation is illustrated in Fig. 4. The wave is steered in the
broadside and encounters an obstacle that is 10% the size of the aperture,
except in the case of the Airy-like beam which curves around a much larger,
indefinitely extending obstacle. This also highlights additional possibilities
in blockage mitigation with THz Airy beams.

With beamforming, the radiation intensity behind the obstacle is significantly
diluted and the wave is blocked. However, other cases show certain resilience
to blockage. With beamfocusing, the radiation from the portion of the aperture
not blocked by the obstacle can still converge to focus at the desired focal
point. However, this is not technically “self-healing”, as the focal spot had
not yet been formed; were the obstacle closer to the focal point, the blockage
would be severe. In contrast, the self-healing nature of Bessel beams _reforms
the beam_ beyond the obstacle and the radiation intensity is almost
unaffected. The performance of the Airy-like beam also remains impervious to
the blockage since the beam can be designed to curve around this obstruction.

###  IV-C Increased System Capacity

The exploitation of wavefronts that carry OAM, such as higher-order Bessel
beams in the near field, allows the creation of perfect parallel orthogonal
channels. In Fig. 5, we show the required bandwidth to achieve a point-to-
point 1 Tbps wireless link for different OAM configurations and various orders
of quadrature amplitude modulation (QAM), assuming no distortions on the
spiral phase for illustrative purposes.

Figure 5: Required bandwidth for a 1 Tbps link capacity with OAM-Multiplexing
enabled.

It is observed from Fig. 5 that the demand on the bandwidth can be
significantly reduced when both OAM and higher-order modulations are used. For
example, utilizing 32 beams, each with a distinct OAM mode, the required
bandwidth is less than 10 GHz, even with 16-QAM. This can be further reduced
down to 3 GHz with a more sophisticated 1024-QAM modulation. The latter could
not only facilitate the spectrum utilization and the co-existence of
prospective THz links with other services (i.e., the Earth exploration
satellites) but also allows to avoid certain wideband-specific challenges
discussed earlier.

##  V Research Challenges and Opportunities

The efficient harnessing of novel wavefronts and their attractive properties
within widespread 6G and beyond scenarios is a wide research area with many
open challenges:

####  V-1 Near-field THz channel

The work should start with a better understanding of the near-field THz
channel itself. Solving Maxwell’s equations numerically for all the possible
THz channel configurations is computationally extensive in many setups. Hence,
there is a need to develop THz-specific channel models or approximations for
near-field propagation in target setups, such as indoor and outdoor wireless
access. These could not only account for the equivalent path loss in the near
field but also scattering and reflection effects from possible complex
propagation environments in the near field.

####  V-2 Near-field THz sensing

Adding sensing capabilities to wireless access systems is today considered one
of the major innovations toward 6G. Nonetheless, the peculiarities of near-
field propagation demand updating the relevant metrics and models. In
addition, the discussed wavefronts (and other counterparts) have been
previously utilized in sensing and imaging applications in the optical domain.
When it comes to THz imaging and sensing, as well as the promising research
topic of joint communications and sensing (JCS), the applicability and impact
of wavefront engineering must be explored with the specifics of the underlying
beam profile taken into account. Here, an interesting research problem is the
fact that one of the important metrics in sensing, i.e., the radar cross-
section (RCS) is, strictly speaking, defined in a far-field sense. While
recent findings on near-field RCS utilizing the spherical wave model are
summarized in [13], evaluating the effectiveness of the developed near-field
model for the other wavefronts is a pressing task. Further, in the advent that
the model is not valid ubiquitously, developing a general solution that
captures the significant complexity of these beams for near-field RCS could be
a promising research avenue.

####  V-3 Wavefront Distortions in Ultra-broadband THz Links

One of the fundamental expectations in practical THz systems is enabling
_broadband_ communications. Wavefront engineering is a narrow-band technique:
a targeted feature of a wavefront does not extend uniformly across a large
bandwidth, as the same phase delay does not generate the same time delay that
truly governs wave propagation, superposition, and final beam shape. True time
delay lines can be utilized to counter this problem through a frequency-
dependent phase profile that makes the time delay uniform. However, the
required architecture is more complex [7].

Here, OAM multiplexing can help boost the system capacity, matching optimal
(uncorrelated) MIMO. In addition, OAM multiplexing/demultiplexing is less
demanding from a signal processing perspective compared to THz MIMO, since, as
mentioned in Sec. III-F, the receiver only needs the correct phase profile
across its aperture. Nonetheless, if the signal’s bandwidth within the OAM
mode is large enough, the generated OAM modes are no longer pure and can lead
to crosstalk or interference between multiplexed OAM streams. Therefore,
further advances are needed to properly characterize these new performance-
complexity trade-offs, such as by characterizing and quantifying this
spillover effect.

####  V-4 Mobile THz systems

Mobility is one of the principal challenges when it comes to operating in
practical, everyday scenarios. Both Bessel and Airy beams propagate in the
near field, and thus the application of beamsteering techniques from the far
field is also valid here. However, it has been shown that OAM-carrying beams,
when steered, have a spillover effect that results in crosstalk (interference)
with neighboring OAM modes. Additionally, to apply both steering and wave-
generating phase profiles, the phase shifters need to be significantly
complex; 4-bit (16-level) phase shifters are usually sufficient for almost all
applications [7], specifying the phase resolution threshold. On the positive
side, the blockage mitigation capabilities may relax certain mobility-related
requirements such as the degree of THz multi-connectivity [1].

####  V-5 Interference Modeling

Most, if not all, studies on interference modeling for mmWave and THz
communications assume plane waves [14]. Thus, the results and observations
from prior studies on interference in THz networks must be revisited to
account for the wavefront specific properties.These findings would be of
paramount importance for efficient frequency and space reuse, medium access
control (MAC), and resource allocation in prospective multi-user THz networks.

####  V-6 Mechanisms to Generate Wavefronts

The generation and detection of custom wavefronts can be achieved with a lower
complexity of the digital backends preceding the array apertures.
Particularly, holographic beamforming architectures may significantly reduce
the power consumption and cost in developing large-scale arrays [15].
Moreover, the extreme scalability of IRSs can constitute large intelligent
surfaces (LISs), which can extend the applicability of these wavefronts to
large-scale cellular deployments.

This setup, however, depends on the accuracy of the phase-shifting elements
(as discussed above) and the spatial resolution of the array [3]. Utilizing a
metasurface approach with sub-wavelength radiating elements could increase the
spatial resolution [3]. Nonetheless, the radiation response of such
metasurfaces is dependent on the mutual coupling between these densely packed
elements; accurately applying the phase profile for a desired wavefront is
challenging.

Clearly, for dynamic wavefront customization, the utilization of lenses and
non-reconfigurable dielectrics will prove insufficient. However, the
significant recent advances in reconfigurable arrays for THz communications
indicate that dynamic and advanced wavefront engineering will be possible in
the 6G and beyond landscape for improved mobility.

####  V-7 Energy Efficiency

The energy cost of THz communications is still a significant challenge,
primarily due to the lack of efficient THz radiation generation in comparison
to power generation at lower frequencies. In addition to significant advances
in device technology [1], wavefront engineering can also be a potential
solution since the beam energy can be better focused towards the intended
directions. The substantial energy challenge of blockages can be addressed by
leveraging the self-healing and self-accelerating properties of Bessel and
Airy beams. Still, comparative research is required to identify which of the
beams to use to address this problem.

##  VI Conclusion

With continuous growth in user demands, the THz band is envisioned to be
explored for both high-rate data exchange and high-precision sensing in next-
generation wireless systems. However, the dogmas of conventional RF systems
from the microwave and mmWave bands cannot hold infinitely for higher
frequencies. Particularly, the THz physical layer must have a symbiosis of
communication and wave theory. Here, the design of novel custom THz wavefronts
discussed in this article is one of the key enablers for efficient practical
THz communication systems that will constitute an inherent part of the 6G and
beyond-6G landscapes.

## References

  * [1] I. F. Akyildiz _et al._ , “Terahertz band communication: An old problem revisited and research directions for the next decade,” _IEEE Trans. on Commun._ , vol. 70, no. 6, pp. 4250–4285, June 2022.
  * [2] C. Han, L. Yan, and J. Yuan, “Hybrid beamforming for terahertz wireless communications: Challenges, architectures, and open problems,” _IEEE Wireless Commun._ , vol. 28, no. 4, pp. 198–204, Aug. 2021.
  * [3] M. Di Renzo _et al._ , “Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come,” _EURASIP J. on Wireless Commun. and Netw._ , vol. 2019, no. 1, pp. 1–20, May 2019.
  * [4] H. Zhou, X. Su, A. Minoofar, R. Zhang, K. Zou, H. Song, K. Pang, H. Song, N. Hu, Z. Zhao _et al._ , “Utilizing multiplexing of structured THz beams carrying orbital-angular-momentum for high-capacity communications,” _Optics Express_ , vol. 30, no. 14, pp. 25 418–25 432, June 2022.
  * [5] C. A. Balanis, _Antenna theory: Analysis and design_. John Wiley & Sons, 2016.
  * [6] M. Cui, Z. Wu, Y. Lu, X. Wei, and L. Dai, “Near-Field MIMO Communications for 6G: Fundamentals, Challenges, Potentials, and Future Directions,” _IEEE Communications Magazine_ , vol. 61, no. 1, pp. 40–46, Jan. 2023.
  * [7] D. Headland _et al._ , “Tutorial: Terahertz beamforming, from concepts to realizations,” _Apl Photonics_ , vol. 3, no. 5, Feb. 2018.
  * [8] J. Durnin, J. Miceli, and J. H. Eberly, “Comparison of Bessel and Gaussian beams,” _Optics letters_ , vol. 13, no. 2, pp. 79–80, Feb. 1988.
  * [9] G. Siviloglou _et al._ , “Observation of accelerating Airy beams,” _Phys. Review Letters_ , vol. 99, no. 21, p. 213901, Nov. 2007.
  * [10] S. Li and J. Wang, “Adaptive free-space optical communications through turbulence using self-healing Bessel beams,” _Scientific reports_ , vol. 7, no. 1, pp. 1–8, Feb. 2017.
  * [11] L. Allen _et al._ , “Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes,” _Physical review A_ , vol. 45, no. 11, p. 8185, June 1992.
  * [12] V. Petrov, D. Bodet, and A. Singh, “Mobile near-field terahertz communications for 6G and 7G networks: Research challenges,” _Frontiers in Communications and Networks_ , vol. 4, p. 6, Mar. 2023.
  * [13] T. Dogaru, “Near-field target scattering characterization and radar modeling,” US Army Combat Capabilities Development Command, Army Research Laboratory, Tech. Rep., May 2021.
  * [14] H. Shokri-Ghadikolaei, C. Fischione, and E. Modiano, “Interference model similarity index and its applications to millimeter-wave networks,” _IEEE Trans. on Wireless Commun._ , vol. 17, no. 1, pp. 71–85, Jan. 2018\.
  * [15] R. Deng _et al._ , “Reconfigurable holographic surface: Holographic beamforming for metasurface-aided wireless communications,” _IEEE Trans. on Vehicular Tech._ , vol. 70, no. 6, pp. 6255–6259, May 2021.

## Authors’ Biographies

Arjun Singh<EMAIL_ADDRESS>is an Assistant Professor in the Department
of Engineering at the State University of New York Polytechnic Institute,
Utica, NY, USA. His research interests include realizing THz band wireless
communications for 6G and beyond.

Vitaly Petrov<EMAIL_ADDRESS>is an Associate Research Scientist,
Northeastern University, Boston, MA, USA. His research interests include THz
band communications for 6G and beyond.

Hichem Guerboukha<EMAIL_ADDRESS>is a Postdoctoral Research
Fellow with the School of Engineering, Brown University, Providence, RI, USA.
His research interests include THz instrumentation and waveguides, THz
computational imaging, and THz communications.

Innem V.A.K. Reddy<EMAIL_ADDRESS>is a dual-Ph.D. graduate from the
University at Buffalo, SUNY, NY, USA, and King Abdullah University of Science
and Technology, Saudi Arabia. His current research interests are in the field
of micro-optics, beam shaping, 3D nanofabrication, optical tweezers,
plasmonics, and THz communications.

Edward W. Knightly<EMAIL_ADDRESS>is the Sheafor-Lindsay Professor of
Electrical and Computer Engineering at Rice University, Houston, TX, USA. His
research interests are in experimental wireless networking and wireless
network security.

Daniel M. Mittleman<EMAIL_ADDRESS>is a Professor in the School
of Engineering at Brown University. His research interests involve the science
and technology of THz radiation.

Josep Miquel Jornet<EMAIL_ADDRESS>is an Associate Professor in
the Department of Electrical and Computer Engineering at Northeastern
University. His research interests are in THz communications and wireless
nano-bio-communication networks.