# Automated Layout Design and Control of Robust Cooperative Grasped-Load
Aerial Transportation Systems

Carlo Bosio, Jerry Tang, Ting-Hao Wang, and Mark W. Mueller The authors are
with the Dept. of Mechanical Engineering, University of California, Berkeley
{c.bosio, jerrytang, wtyngh<EMAIL_ADDRESS>
###### Abstract

We present a novel approach to cooperative aerial transportation through a
team of drones, using optimal control theory and a hierarchical control
strategy. We assume the drones are connected to the payload through rigid
attachments, essentially transforming the whole system into a larger flying
object with “thrust modules” at the attachment locations of the drones. We
investigate the optimal arrangement of the thrust modules around the payload,
so that the resulting system is robust to disturbances. We choose the
$\mathcal{H}_{2}$ norm as a measure of robustness, and propose an iterative
optimization routine to compute the optimal layout of the vehicles around the
object. We experimentally validate our approach using four drones and
comparing the disturbance rejection performances achieved by two different
layouts (the optimal one and a sub-optimal one), and observe that the results
match our predictions.

##  I INTRODUCTION

The increasing demand for effective and autonomous transportation in
applications such as logistics, search and rescue operations, and load lifting
has increased interests in efficient aerial transportation solutions [1].
Being able to apply aerial robots to logistics and transportation would not
only reduce costs and enhance efficiency, but also allow to reduce the burden
on ground infrastructures. The application of UAVs to these fields is
challenging, mainly due to their payload and battery life limitations. We
address the problem of using multiple drones for payload lifting and
transportation. A lot of work has been carried out on design, control and path
planning for single drone transportation systems, some examples of which are
[2, 3, 4, 5, 6, 7]. However, practical limitations on dimensions, vehicle
complexity, load capacity, and costs limit the application of such
technologies.

The widespread use and cost-effectiveness of smaller vehicles, like
quadcopters, have made them the preferred choice for practical applications
[8]. However, using multiple smaller vehicles in a cooperative manner, along
with potential payload interactions (e.g. aerodynamic interactions,
vibrations), also introduce greater complexity in terms of control and
trajectory planning for cooperative aerial transportation systems [9, 10].
Connecting the vehicles to the payload through tethers is a popular choice,
and allows to keep distance between drones and payload. The dynamics and
control of tethered systems have been extensively investigated and yielded
interesting results, e.g. [11]. Other examples are [12], which employed an
interconnected structure to enhance the stability of the tethered system, or
[13], which implemented a decentralized control scheme not relying on
assumptions about the payload’s state. [14] presented a cooperative control
approach which showed significant disturbance rejection performances.

Figure 1: Picture of four drones cooperatively carrying a single panel
payload.

Compared to the challenges of managing tethered or moving payloads, which can
be complex due to their internal dynamics and limited maneuverability, rigid
attachments are often the preferred option for transportation purposes. Some
work has been carried out also along this direction. [15] studied the control
for cooperative payload stabilization with a team of drones equipped with
grippers. [16] focused on cooperative grasped transportation control using
onboard visual and inertial sensing, while [17] implemented an adaptive
control framework based on a Kalman Filter for systems with uncertainty about
drone location and payload parameters. [18] explores the possibility of
transporting flexible payloads, taking into account internal deformations and
dynamics and aiming to minimize deformations.

Our research addresses the synthesis of rigidly attached transportation
systems. The core concept is to address both the design layout and the optimal
control of the combined drone-payload system as one unified problem. We do
this by formulating and solving an optimization problem derived from robust
control theory for a first order model of the system. In particular,
$\mathcal{H}_{2}$ control proved to be a promising approach to maximize
disturbance rejection performances [19]. The main contribution of this
research is an automated design tool to determine the optimal thrust module
arrangement for payload stability and robustness. Experimentally validated
with four drones, our method confirms the effectiveness of the optimal layout
in rejecting disturbances.

The paper is structured as follows. In Section II we describe the mathematical
framework employed to compute the optimal placement locations and the linear
feedback law to control the system. We then describe the control
infrastructure adopted in our hardware setup. In Section III we show some
example outputs of our optimization routine applied to different payload
shapes, masses and number of drones used for transportation. In Section IV we
introduce the experimental setup and the flight tests which were carried out
to validate the theoretical predictions. Finally, in Section V we discuss the
results, relevance of the work and future developments.

##  II MATERIALS AND METHODS

The problem we solve is twofold. We i) optimize the placement of the drones on
the object to be transported, so as to maximize the system’s disturbance
rejection performance. This involves solving an optimization problem where the
cost function is derived from $\mathcal{H}_{2}$ control theory, and the
decision variables are the drones’ attachment locations. We then ii) address
the constrained drones dynamics through a hierarchical control strategy of the
integrated drone-payload system. The motion of the system is controlled using
the gain obtained from the same $\mathcal{H}_{2}$ control used for layout
optimization. Each drone (thrust module) is sent thrust commands. An example
picture of a team of drones carrying a payload with rigid attachments is shown
in Fig. 1.

###  II-A Layout Optimization

We consider the transportation of a planar object, i.e. with a geometry which
can be seen as the extrusion of a planar shape along the direction orthogonal
to the plane. This assumption simplifies the mathematical derivation but does
not limit significantly the space of possibilities: many objects, as delivery
packages, boxes, and pallets are all included in this category. This
assumption, since the drones can all lie on the same plane, is useful because
it eventually simplifies both the control of the system and the layout
optimization problem. In fact, each drone attachment location is described by
one variable, i.e. its position along the perimeter of the payload shape. We
also assume that the inertial properties of the payload are known. With
reference to Fig. 2, the placement variables, are defined as
$\boldsymbol{\theta}=[\theta_{i},...,\theta_{N}]^{T}$. The reference frame is
defined at the geometric center of the in-plane object geometry. The drone-
payload system is itself a rigid body, and can therefore be described by a
state vector
$\mathbf{x}=[\mathbf{p}^{T},\dot{\mathbf{p}}^{T},\boldsymbol{\gamma}^{T},\boldsymbol{\omega}^{T}]^{T}\in\mathbb{R}^{12}$,
where $\mathbf{p}$ and $\dot{\mathbf{p}}$ are the position and velocity of the
center of the object with respect to a global world reference frame,
$\boldsymbol{\gamma}$ is the Roll-Pitch-Yaw triplet defining the orientation
of the local frame, and $\boldsymbol{\omega}$ the angular velocity with
respect to the local body frame. Furthermore, each drone is seen as a thrust
module providing four inputs to the system. Therefore, if $N$ drones are
deployed, we have a thrust input vector $\mathbf{u}\in\mathbb{R}^{4N}$. Due to
unavoidable hardware limitations, we assume that each thrust component is
limited to the interval $[u_{l},u_{h}]$. It is possible to linearize the rigid
body dynamic model about hover configuration, and decompose the inputs in
$\mathbf{u}=\bar{\mathbf{u}}+\mathbf{u^{\prime}}$, in which $\bar{\mathbf{u}}$
are the feedforward thrusts and $\mathbf{u}^{\prime}$ are the first order
components of the linearized model. We do the same for the state vector
$\mathbf{x}=\bar{\mathbf{x}}+\mathbf{x}^{\prime}$. We model the disturbances
applied to the system as a random variable
$\mathbf{d}=[\mathbf{F}_{d}^{T},\mathbf{T}_{d}^{T}]^{T}$, in which
$\mathbf{F}_{d}$ and $\mathbf{T}_{d}$ represent white noise force and a torque
applied to the origin of the local body reference frame. The problem is
formulated as the minimization of the probability of infeasible inputs when
the system is subject to the defined random disturbances.

Figure 2: Schematics of drones attached to a planar payload. Representation of
the local body reference frame, placement variables $\boldsymbol{\theta}$, and
disturbance actions $F_{d}$ and $T_{d}$.

We analyze the integrated drone-payload system as a whole. We can write its
continuous-time linearized dynamics as:

$\dot{\mathbf{x}}^{\prime}=A\,\mathbf{x}^{\prime}+B(\boldsymbol{\theta})\,\mathbf{u}^{\prime}+B_{d}(\boldsymbol{\theta})\,\mathbf{d}.$
(1)

The dynamic matrix $A$ is obtined by linearization of the rigid body dynamics
equations about hovering and does not depend on the placement variables
$\boldsymbol{\theta}$. The matrices $B(\boldsymbol{\theta})$ and
$B_{d}(\boldsymbol{\theta})$ depend on $\boldsymbol{\theta}$ either directly
or through the inertia of the system (which also depends on the placement
variables). With reference to the $\mathcal{H}_{2}$ optimal control literature
[20], we define an additional auxiliary variable:

$\mathbf{z}=C\mathbf{x}+D\mathbf{u}.$ (2)

The goal is to find the feedback gain $K^{*}$, such that

$\mathbf{u}=\bar{\mathbf{u}}-K^{*}(\boldsymbol{\theta})\mathbf{x}^{\prime}$
(3)

and minimize the integral $\mathcal{H}_{2}$ control cost
$J=\sum_{k=1}^{6}\left[\int_{0}^{\infty}\mathbf{z}^{T}\mathbf{z}\,dt\,\,\,\text{s.t.}\,\,\,d=\mathbf{e}_{k}\delta(t)\right]$,
where $\mathbf{e}_{k}$ is the $k$-th element of the euclidean base and
$\delta(t)$ is Dirac’s delta function. It is a known result that the optimal
feedback controller is obtained as

$K^{*}(\boldsymbol{\theta})=(D^{T}D)^{-1}B(\boldsymbol{\theta})^{T}S_{1}(\boldsymbol{\theta})$
(4)

where $S$ is the solution of the algebraic Riccati equation:

$A^{T}S_{1}+S_{1}A-S_{1}B(D^{T}D)^{-1}B^{T}S_{1}+C^{T}C=0$ (5)

in which we omitted the dependency on $\boldsymbol{\theta}$ for clarity of
writing. The feedback gain does not depend on the matrix $B_{d}$, which
represents how the disturbance comes into the dynamics. The optimal cost value
turns out to be:

$J^{*}=\text{tr}(B_{d}^{T}S_{1}B_{d}).$ (6)

The optimal cost value does indeed depend on the matrix $B_{d}$. This is a
valid cost function which can be used to optimize the drones’ attachment
locations. However, it does not guarantee that the inputs generated by eq. 3
are feasible. Therefore, we make an additional step: we compute the covariance
$S_{2}(\boldsymbol{\theta})$ of the state vector in a condition in which the
optimal feedback is active and the system is subject to white noise
disturbance. This can be computed using the optimal continuous-time state
observer theory, which leads to the equation:

$A_{f}S_{2}+S_{2}A_{f}^{T}+B_{d}B_{d}^{T}=0$ (7)

where we defined
$A_{f}(\boldsymbol{\theta})=A-B(\boldsymbol{\theta})K^{*}(\boldsymbol{\theta})$
as the dynamics matrix with linear feedback. Note that the assumption of white
noise disturbance makes the process ergodic, and therefore the state sample
mean is equivalent to the temporal mean. The matrix
$S_{2}(\boldsymbol{\theta})$ is then an evaluation of the covariance of the
first order state linearization $\mathbf{x}^{\prime}$ when subject to white
noise disturbance. At this point we have the feed-forward thrust values, which
are computed from the hovering equilibrium condition, and the covariance
matrix of the inputs, obtained through the feedback law
$K^{*}(\boldsymbol{\theta})^{T}S_{2}(\boldsymbol{\theta})K^{*}(\boldsymbol{\theta})$.
The goal is to maximize the probability of feasible inputs. If we set the
element-wise thrust lower bound to be $\mathbf{u}_{l}=u_{l}\mathbf{1}_{4N}$
and the upper bound to be $\mathbf{u}_{h}=u_{h}\mathbf{1}_{4N}$, where
$\mathbf{1}_{4N}$ indicates the $4N$-dimensional vector of ones. We can then
write the optimization problem we are trying to solve as:

$\max_{\boldsymbol{\theta}}F(\mathbf{u}_{h})-F(\mathbf{u}_{l})$ (8)

where $F(\cdot)$ is the cumulative distribution over thrust inputs, and its
mean and covariance depend on $\boldsymbol{\theta}$. It is possible to compute
this cost function through sampling, but the sampling approach would lead to a
noisy cost function whose variance decreases linearly with the number of
samples used. An additional limitation of sampling is that it needs a
distribution to draw samples from, and this would have to be assumed. In
addition, sampling methods are computationally expensive and lead to long
solution times. Therefore, we choose to leverage the concept of minimum
Mahalanobis distance [21] of the mean (feedforward thrusts) to the thrust
bounds. The Mahalanobis distance of a point $\mathbf{u}$ to a distribution
with mean $\bar{\mathbf{u}}$ and covariance $\Sigma_{\mathbf{u}}$ is defined
as:

$d_{M}(\mathbf{u};\bar{\mathbf{u}},\Sigma_{\mathbf{u}})=\sqrt{(\mathbf{u}-\bar{\mathbf{u}})\Sigma_{\mathbf{u}}^{-1}(\mathbf{u}-\bar{\mathbf{u}})}.$
(9)

We can extend the concept of Mahalanobis distance to a hyperplane as the
minimum Mahalanobis distance of a point of the hyperplane from the
distribution. A two dimensional example is shown in Fig. 3. Since in our case
we have a lower and an upper bound for each input, there are $8N$ different
hyperplanes to check in the input space. Therefore, going back to the goal of
maximizing the probability of feasible inputs, we aim to maximize the
Mahalanobis distance of the mean (feedforward thrusts) to the thrust bounds.
The optimization problem can be formulated as:

$\displaystyle\max_{\boldsymbol{\theta}}\biggl{\\{}\,\,\min_{k}$
$\displaystyle\left[\min_{u_{a}\in\\{u_{l},u_{h}\\}}\biggl{(}\,\,\min_{\mathbf{u}}d_{M}\Bigl{(}\mathbf{u};\bar{\mathbf{u}}(\boldsymbol{\theta}),\Sigma_{\mathbf{u}}(\boldsymbol{\theta})\Bigl{)}\,\,\biggl{)}\,\,\right]\,\,\biggl{\\}}$
(10)
$\displaystyle\text{s.t.}\,\,\mathbf{u}_{l}\leq\mathbf{u}\leq\mathbf{u}_{h}\,\,\,\text{(element-
wise)}$ $\displaystyle\,\,\,\,\,\,\,\,\,u^{(k)}=u_{a}$

In which with the superscript $u^{(k)}$ we refer to the $k$-th component of
the vector $\mathbf{u}$. We are then dealing with a nested optimization
problem. Trying to unpack what the routine is about: we are looking for the
value $\boldsymbol{\theta}^{*}$ which maximizes the minimum Mahalanobis
distance of the thrust bounds from the feedforward thrust. This minimum
distance is computed by evaluating it on each bound hyperplane separately,
hence the minimization over $k$ (component-wise) and over
$u_{a}\in\\{u_{l},u_{h}\\}$, where the subscript $a$ stands for “active”. The
inner minimization over $\mathbf{u}$ is necessary because we are considering
all the points on the active bound hyperplane, and are interested in the one
with minimum distance from $\bar{\mathbf{u}}$. The optimization routine is
also described in Algorithm 1. On the computational cost level, there are no
particular expensive operations. Solving Riccati equations is relatively fast,
and minimizing the Mahalanobis distance from a hyperplane is a Quadratic
Program (QP). Therefore, for each update of $\boldsymbol{\theta}$ in the
optimization routine, $8N$ QPs are solved.

Figure 3: Two dimensional example of minimum Mahalanobis distance to input
bounds. The point of minimum distance is the point of tangency between the
largest internal ellipsoid and the closest bounding hyperplane. Algorithm 1
$\mathcal{H}_{2}$ based drone placement optimization

$N,\text{drones inertia},\text{object inertia},\text{object shape}$

$\text{weight matrices}\,\,C\,\,\text{and}\,\,D$

Initial guess for $\boldsymbol{\theta}$

while Optimization not completed do

Compute overall system inertia

Compute $A$, $B(\boldsymbol{\theta})$, $B_{d}(\boldsymbol{\theta})$ and
feedforward $\bar{\mathbf{u}}$

Compute $S_{1}(\boldsymbol{\theta})$ and $K^{*}(\boldsymbol{\theta})$ from eq.
5

Compute $S_{2}(\boldsymbol{\theta})$ from eq. 7

Compute input covariance as $\Sigma_{\mathbf{u}}\leftarrow K^{*T}S_{2}K^{*}$

Initialize $d_{M}^{*}$

for each thrust input component $k$ do

for $u_{a}$ in $\\{u_{l},u_{h}\\}$ do

Solve QP
$d^{\prime}\leftarrow\min_{\mathbf{u}}d_{M}(\mathbf{u};\bar{\mathbf{u}},\Sigma_{\mathbf{u}})$
as in eq.

10 s.t. $\mathbf{u}_{l}\leq\mathbf{u}\leq\mathbf{u}_{h}$, $u^{(k)}=u_{a}$

if $d^{\prime}<d_{M}^{*}$ then

$d_{M}^{*}\leftarrow d^{\prime}$

end if

end for

end for

Update $\boldsymbol{\theta}$

end while

Figure 4: Control infrastructure diagram

###  II-B Control Infrastructure

The dynamics of motion of each drone is highly constrained due to the rigid
attachments. Therefore, using a traditional drone controller would be
challenging. However the integrated system corresponds itself to a rigid body
(neglecting effects due to internal elasticity and vibrations). Therefore,
approaching the overall dynamics and control as a whole (as was done for the
placement optimization) is preferable. This strategy leads to a hierarchical
control infrastructure, in which a single estimator is used for the overall
motion and the drones receive thrust commands individually through $N$
parallel communication channels. Furthermore, this control strategy comes
naturally from the layout optimization computation, during which the optimal
LQR gains have been computed. A block diagram of the suggested control
infrastructure is shown in Fig. 4. The overall behavior of the controller is
dominated by the choice of the matrices $C$ and $D$ (eq. 2), which directly
impact the computation of the $\mathcal{H}_{2}$ cost. For simplicity, the $C$
matrix was chosen to just generate a contribution to the $\mathcal{H}_{2}$
cost from the position of the system’s center of mass and yaw angle (i.e.
$C\in\mathbb{R}^{4\times 12}$). After evaluating the system’s behavior in
simulation, we set $C(1,1)=C(2,2)=0.5$ (corresponding to the position on a
horizontal plane $x-y$), $C(3,3)=10$ (corresponding to vertical position $z$),
$C(4,9)=50$ (corresponding to yaw angle), and zero for all the other
components. The input contribution was set to the $2N$-dimensional identity
matrix $D=I_{2N}$.

###  II-C Hardware Setup

The vehicles we used as a reference for layout design and eventually for
experiments are quadcopters equipped with a PX-4 flight controller. The
relevant vehicle parameters can be found in Tab. I.

drone size (motor to motor) | $330\,\mathrm{mm}$ | propeller diameter | $114.5\,\mathrm{mm}$
---|---|---|---
drone frame mass | $525\,\mathrm{g}$ | battery mass | $437\,\mathrm{g}$

TABLE I: Table of parameters

The battery capacity ($5300\,\mathrm{mAh}$) has been chosen to guarantee $10$
minutes of flight time at maximum payload. The system theoretically is able to
lift an external payload of $\sim 0.8\,\mathrm{kg}$ per drone. It is
interesting to note that the main contribution to the overall inertia of the
system is due to the drones, since they are attached to the perimeter. It is
interesting to note that contributing a greater inertia than the actual
payload is in line with all the ground and aerial transportation systems in
use nowadays.

##  III RESULTS

The optimization framework which we developed allows us to determine the
optimal placement of drones for transportation across many different shapes,
masses and number of drones used. In Fig. 5 we show some examples of optimal
placements using the vehicles described in Tab. I as a reference. It is
interesting to observe that the drones are placed in a way which seems to
trade off two different contributions. On the one hand, the area of the
polygon defined by the drones’ attachment points tends to be maximized in
order to maximize control authority. On the other hand the distances of the
attachment points to the system’s center of mass tend to be equalized. This is
an effect of employing the Mahalanobis distance as a metric. In fact, even if
the control authority is maximized (i.e. the input covariance
$\Sigma_{\mathbf{u}}$ has smaller eigenvalues), we also need the drones to
have a comparable margin between feedforward thrusts $\bar{\mathbf{u}}$ and
bounds, otherwise the Mahalanobis distance objective is decreased.

Figure 5: Optimal layouts for various object shapes and masses using different
number of drones.

##  IV EXPERIMENTAL VALIDATION

For the flight experiments the drones were attached to two different square
panels of different masses and sizes. A heavier panel (A) of side length
$0.61\,\mathrm{m}$ and mass $2.54\,\mathrm{kg}$ was used to test the payload
capacity, and a lighter one (B) of side length $0.45\,\mathrm{m}$ and mass
$1.088\,\mathrm{kg}$ to test the disturbance rejection performances. The rigid
attachments have been obtained by screwing directly the drones to the panels
through wood dowels. The experiments are conducted in an indoor flight space
equipped with a Motion Capture system, and markers for motion tracking are
placed directly on the panels.

###  IV-A Payload lifting

Panel A was used to test the payload capacity of the system. The four drones
were attached to the corners of the panel and a hovering position was
commanded. The system was able to successfully take off, stabilize, hover and
land safely. An important observation to highlight is the close proximity of
the feedforward thrusts $\bar{\mathbf{u}}$ to the upper thrust limit
$\mathbf{u}_{h}$. Given the significant mass of Panel A ($\sim
2.5\,\mathrm{kg}$) compared to the drone masses ($\sim 0.95\,\mathrm{kg}$
including battery), the narrow margin causes the system to have limited
disturbance response capabilities, especially within the spatial limits of the
flight space available for experiments.

###  IV-B Disturbance Rejection

The four drones were attached to panel B using two distinct configurations: an
optimal one, placed at the corners, and a sub-optimal one. These
configurations can be seen in Fig. 6. As a disturbance, a weight of
$0.5,\mathrm{kg}$, or approximately $10\%$ of the system’s total mass, was
attached to a corner of the panel during the flight, specifically under Drone
1 (Fig. 6). The mass attachment is done through a magnet to guarantee the
repeatability of the experiment. Both the state and input data were collected
and compared to validate their consistency with the layout optimization
framework. Representative response plots over time are presented in Fig. 7.

From the data, it is clear that the choice of drone layout has a significant
impact on disturbance rejection performance. The optimal layout results in
lower oscillations in the attitude of the vehicle, specifically in the pitch
and roll angles, and also recovers more quickly from disturbances. In
comparison, the sub-optimal layout sees its peak deviation from hovering
conditions nearly double the one of the optimal setup. Another observation is
related to thrust behavior after adding the payload. In the sub-optimal
configuration, the increase in thrust following the addition of the payload is
more noticeable. Additionally, the steady-state thrust from the drone
positioned closest to the added payload is also higher in the sub-optimal
setup. This suggests that, despite facing disturbances of the same magnitude,
the thrust distribution becomes more uneven in the sub-optimal case,
emphasizing the benefits of an optimal drone layout. It is important to note
that in the sub-optimal configuration, the model and controller gains were
recomputed from scratch to maximize performances (i.e. we did not use the same
controller as in the optimal configuration).

Figure 6: Attachment layouts tested on panel B for disturbance rejection
comparison. Left: optimal configuration, right: sub-optimal configuration.

##  V CONCLUSIONS

In this study, we introduced a novel approach to cooperative aerial
transportation with a team of drones, employing optimal control theory and a
hierarchical control strategy. This methodology offers an original approach in
which the layout design and control problems are solved in a coupled way
through an optimization routine derived from $\mathcal{H}_{2}$ control theory.

Figure 7: Response to step disturbance over time. The mass is attached
instantly (through a magnet) underneath Drone 1. The attachment time is
indicated with a vertical dashed line. On the top row, attitude response. On
the bottom row, total quadcopter thrusts. Left: optimal configuration, right:
sub-optimal configuration

Our optimization framework allows to cater our system to different payload
shapes, inertias, and drone counts. The hierarchical control structure offers
a systematic and simple way to deal with the dynamics viewing the system as a
single, integrated rigid body. The choice of the weight matrices $C$ and $D$
represents the freedom given to the control designer to obtained the desired
performances. Overall, the main contributions of this paper are:

  * •

An automated tool to optimally arrange thrust modules on a payload to maximize
robustness of the integrated system;

  * •

A control framework general enough to be able to fly different payload shapes
using masses the preferred number of drones.

We also emphasize the idea that dealing with system design and control in a
coupled manner simplifies the control complexity and increases the overall
system performances compared to a sequential siloed approach (control design
following system design).

Experimental results validate the effectiveness of our approach. Flight tests
under disturbance highlighted the stability and robustness of the system, and,
most importantly, the agreement between the performance predicted by our
optimization tool and experimental evidence.

Future developments for this work include:

  * •

Exploring more advanced control strategies and infrastructures. At this stage,
for example, the vehicles’ onboard sensor data were not used.

  * •

Removing the assumption of known payload inertia and developing a control
scheme under this uncertainty component.

  * •

Using the system as a moving tray that can be used to transport and balance a
payload placed on its top.

In conclusion, this study highlights the coupling between design variables and
control. The physical placement of drones significantly affect the
performances of the controller. Optimizing the layout allows to take more
advantage of optimal control tools rather than applying them at a second
stage, once the design process is finished.

##  VI ACKNOWLEDGMENTS

This work was funded by the Tsinghua-Berkeley-Shenzhen Institute (TBSI). The
authors would also like to thank Roman Ibrahimov and Ruiqi Zhang for their
precious help.

## References

  * [1] D. K. Villa, A. S. Brandao, and M. Sarcinelli-Filho, “A survey on load transportation using multirotor uavs,” _Journal of Intelligent & Robotic Systems_, vol. 98, pp. 267–296, 2020.
  * [2] B. Arbanas, A. Ivanovic, M. Car, T. Haus, M. Orsag, T. Petrovic, and S. Bogdan, “Aerial-ground robotic system for autonomous delivery tasks,” in _2016 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2016, pp. 5463–5468.
  * [3] J. Zeng, P. Kotaru, M. W. Mueller, and K. Sreenath, “Differential flatness based path planning with direct collocation on hybrid modes for a quadrotor with a cable-suspended payload,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 3074–3081, 2020.
  * [4] F. Schiano, P. M. Kornatowski, L. Cencetti, and D. Floreano, “Reconfigurable drone system for transportation of parcels with variable mass and size,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 12 150–12 157, 2022.
  * [5] S. Lee and H. Lee, “Trajectory generation of a quadrotor transporting a bulky payload in the cluttered environments,” _IEEE Access_ , vol. 10, pp. 31 586–31 594, 2022.
  * [6] M. Moshref-Javadi and M. Winkenbach, “Applications and research avenues for drone-based models in logistics: A classification and review,” _Expert Systems with Applications_ , vol. 177, p. 114854, 2021.
  * [7] G. Macrina, L. D. P. Pugliese, F. Guerriero, and G. Laporte, “Drone-aided routing: A literature review,” _Transportation Research Part C: Emerging Technologies_ , vol. 120, p. 102762, 2020.
  * [8] M. W. Mueller, S. J. Lee, and R. D’Andrea, “Design and control of drones,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 5, no. 1, pp. 161–177, 2022.
  * [9] H. Lee, H. Kim, and H. J. Kim, “Planning and control for collision-free cooperative aerial transportation,” _IEEE Transactions on Automation Science and Engineering_ , vol. 15, no. 1, pp. 189–201, 2016.
  * [10] H. Lee, H. Kim, W. Kim, and H. J. Kim, “An integrated framework for cooperative aerial manipulators in unknown environments,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 2307–2314, 2018.
  * [11] C. Meissen, K. Klausen, M. Arcak, T. I. Fossen, and A. Packard, “Passivity-based formation control for uavs with a suspended load,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 13 150–13 155, 2017.
  * [12] H. Rastgoftar and E. M. Atkins, “Cooperative aerial lift and manipulation (calm),” _Aerospace Science and Technology_ , vol. 82, pp. 105–118, 2018\.
  * [13] K. Mohammadi, M. Jafarinasab, S. Sirouspour, and E. Dyer, “Decentralized motion control in a cabled-based multi-drone load transport system,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 4198–4203.
  * [14] J. Geng and J. W. Langelaan, “Implementation and demonstration of coordinated transport of a slung load by a team of rotorcraft,” in _AIAA Scitech 2019 Forum_ , 2019, p. 0913.
  * [15] D. Mellinger, M. Shomin, N. Michael, and V. Kumar, “Cooperative grasping and transport using multiple quadrotors,” in _Distributed Autonomous Robotic Systems: The 10th International Symposium_. Springer, 2013, pp. 545–558.
  * [16] G. Loianno and V. Kumar, “Cooperative transportation using small quadrotors using monocular vision and inertial sensing,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 680–687, 2017.
  * [17] K. Webb and J. Rogers, “Adaptive control design for multi-uav cooperative lift systems,” _Journal of Aircraft_ , vol. 58, no. 6, pp. 1302–1322, 2021.
  * [18] R. Ritz and R. D’Andrea, “Carrying a flexible payload with multiple flying vehicles,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 3465–3471.
  * [19] N. Bucki and M. W. Mueller, “A novel multicopter with improved torque disturbance rejection through added angular momentum,” _International Journal of Intelligent Robotics and Applications_ , vol. 3, no. 2, pp. 131–143, 2019.
  * [20] G. E. Dullerud and F. Paganini, _A course in robust control theory: a convex approach_. Springer Science & Business Media, 2013, vol. 36.
  * [21] P. C. Mahalanobis, “On the generalized distance in statistics,” _Sankhyā: The Indian Journal of Statistics, Series A (2008-)_ , vol. 80, pp. S1–S7, 2018.


# Axis-symmetric Onsager Clustered States of Point Vortices in a Bounded
Domain

Yanqi Xiong  Graduate School of China Academy of Engineering Physics, Beijing
100193, China  Jiawen Chen  Graduate School of China Academy of Engineering
Physics, Beijing 100193, China  Xiaoquan Yu<EMAIL_ADDRESS>Graduate School
of China Academy of Engineering Physics, Beijing 100193, China  Department of
Physics, Centre for Quantum Science, and Dodd-Walls Centre for Photonic and
Quantum Technologies, University of Otago, Dunedin, New Zealand.

###### Abstract

We study axis-symmetric Onsager clustered states of a neutral point vortex
system confined to a two-dimensional disc. Our analysis is based on the mean
field of bounded point vortices in the microcanonical ensemble. The clustered
vortex states are specified by the inverse temperature $\beta$ and the
rotation frequency $\omega$, which are the conjugate variables of energy $E$
and angular momentum $L$. The formation of the axis-symmetric clustered vortex
states (azimuthal angle independent) involves the separating of vortices with
opposite circulation and the clustering of vortices with same circulation
around origin and edge. The state preserves $\rm SO(2)$ symmetry and breaks
$\mathbb{Z}_{2}$ symmetry. We find that, near the uniform state, the rotation
free state ($\omega=0$) emerges at particular values of $L^{2}/E$ and $\beta$.
At large energies, we obtain asymptotically exact vortex density
distributions, whose validity condition gives rise the lower bound of $\beta$
for the rotation free states. Noticeably, the obtained vortex density
distribution near the edge at large energies provides a novel exact vortex
density distribution for the corresponding chiral vortex system.

Keywords: Vortex clusters, Negative temperature, Exact solutions, Quantum
vortices

##  I Introduction

In two-dimensional (2D) fluid turbulence, energy at small scales can transport
to large scales known as inverse energy cascade Eyink and Sreenivasan (2006);
Kraichnan and Montgomery (1980); Tabeling (2002); Boffetta and Ecke (2012).
This process involves formations of large scale vortex patterns. Onsager
explained the formation of large scale structures by studying equilibrium
statistical mechanics of point vortices in a bounded domain. The macroscopic
vortex structure is associated with clustering of like-sign point vortices at
negative temperature Onsager (1949); Eyink and Sreenivasan (2006). These
coherent large structures occur in various systems. Examples are Great Red
Spot in Jupiter’s atmosphere Vasavada  _et al._ (1998), giant vortex clusters
in atomic Bose-Einstein condensates (BECs) Gauthier  _et al._ (2019);
Johnstone  _et al._ (2019); Reeves  _et al._ (2022), and vortex clustering in
quantum fluids of exciton–polaritons Panico  _et al._ (2023).

Clustering phenomena of vortices have attracted much attention Montgomery and
Joyce (1974); Edwards and Taylor (1974); Williamson (1977); Pointin and
Lundgren (1976); Campbell and O’Neil (1991); Smith (1989); Smith and O’Neil
(1990); Yatsuyanagi  _et al._ (2005); Esler and Ashbee (2015); Billam  _et
al._ (2014); Siggia and Aref (1981); Reeves  _et al._ (2013, 2012); Billam
_et al._ (2015); Simula  _et al._ (2014); Gurarie and Chow (2004); Yu  _et
al._ (2016); Salman and Maestrini (2016); Valani  _et al._ (2018); Han and
Tsubota (2018, 2019); Skaugen and Angheluta (2016, 2017); Kanai and Guo
(2021); Patil and Dunkel (2021); Pandit  _et al._ (2017). For a given 2D
domain, searching for the maximum entropy clustered vortex state is at the
center of investigations. For circularly symmetric domains, previous studies
on neutral vortex systems focus on zero angular momentum case Yatsuyanagi  _et
al._ (2005); Yu  _et al._ (2016); Simula  _et al._ (2014); Skaugen and
Angheluta (2016, 2017). The role of finite angular momentum in formation of
clustered states in a neutral vortex system remains less well-explored.

In this paper, we study axis-symmetric clustered vortex states through the
mean field approach. The mean field theory to describe formation of negative
temperature clustered vortex states was formulated systematically by Joyce and
Montgomery Joyce and Montgomery (1973). The mean field equations, which were
obtained via maximizing the entropy of the vortex system, are essential for
analyzing possible clustered states. We consider a neutral vortex system
consisting of equal number of positive and negative vortices confined to a
disc. For given positive vortex number $N_{+}$ and negative vortex number
$N_{-}$, clustered vortex states are specified by energy $E$ and angular
momentum $L$ or their conjugate variables inverse temperature $\beta(E,L)$ and
rotation frequency $\omega(E,L)$. We find that in the limit $\beta\rightarrow
0$, $\omega\rightarrow\infty$ while keeping $\beta\omega$ finite, positive and
negative vortex density distributions are Gaussian distributions centered at
origin and edge, respectively. For rotation free sates ($\omega=0$), we find
asymptotically exact positive and negative vortex density distributions at
large energies. In particular, the one maximized on the edge provides a new
exact solution to the mean field equations for the corresponding chiral vortex
system. The lower bound of $\beta$ is obtained from the validity condition of
the asymptotically exact solutions at high energies, above which rotation free
clustered states exist. To analyze clustered states closed to the uniform
state at low energies, we generalized the perturbation theory, which was
initially developed for chiral systems Smith and O’Neil (1990), to the neutral
case. Using this perturbation theory we find the critical value of $\beta$ for
the onset of the rotation free clustered vortex state, providing an upper
bound of $\beta$.

##  II Model

The point-vortex model describes dynamics of well-separated quantum vortices
in a superfluid at low temperature Ambegaokar  _et al._ (1980), 2D classical
inviscid, incompressible fluids Aref  _et al._ (1999); Siggia and Aref (1981)
and guiding-center plasma Joyce and Montgomery (1973). Negative temperature
states occur due to the bounded phase space of a 2D confined point vortex
system. Above a certain energy the number of available states decreases as the
function of energy and consequently the system becomes more ordered as energy
increases Eyink and Sreenivasan (2006) .

We consider a system consisting a large number of point vortices confined to a
uniform disc with radius $R$. The system is neutral and contains $N_{+}$
positive vortices and $N_{-}=N_{+}$ negative vortices . The Hamiltonian is
Newton (2001)

$\displaystyle H$ $\displaystyle=-\sum_{i\neq
j}\kappa_{i}\kappa_{j}\log|{\boldsymbol{r}}_{i}-{\boldsymbol{r}}_{j}|+\sum_{i,j}\kappa_{i}\kappa_{j}\log\bigg{|}({\boldsymbol{r}}_{i}-\bar{\boldsymbol{r}}_{j})\frac{|{\boldsymbol{r}}_{j}|}{R}\bigg{|}.$
(1)

In a BEC, the Hamiltonian Eq. (1) is measured in unit $E_{0}=\rho
m_{a}\kappa^{2}/4\pi$, where $\rho$ is the superfluid density,
$\kappa=h/m_{a}$ is the circulation quantum and $m_{a}$ is the atomic mass. In
this unit, $\kappa_{i}=\pm 1/N_{\pm}$ and the $1/N$ scaling gives a well-
defined mean field limit Eyink and Spohn (1993); Caglioti  _et al._ (1995).
For a vortex at position ${\boldsymbol{r}}_{j}$, its image locates at
$\bar{{\boldsymbol{r}}}_{j}=R^{2}{\boldsymbol{r}}_{j}/|{\boldsymbol{r}}_{j}|^{2}$
to ensure that the fluid velocity normal to the boundary vanishes. The
Hamiltonian (1) has rotational $\rm SO(2)$ symmetry due to the disc geometry
and $\mathbb{Z}_{2}$ symmetry (invariant under
$\kappa_{i}\rightarrow-\kappa_{i}$). Hereafter we set $R=1$ without loosing
generality.

To investigate formations of large-scale clustered patterns, it is necessary
to consider the continuous effective Hamiltonian in the large $N$ limit Joyce
and Montgomery (1973):

$\displaystyle\ H_{\rm eff}=\frac{1}{2}\int{\rm d}^{2}{\boldsymbol{r}}{\rm
d}^{2}{\boldsymbol{r}}^{\prime}\,\sigma({\boldsymbol{r}})\phi({{\boldsymbol{r}}-{\boldsymbol{r}}^{\prime}})\sigma({\boldsymbol{r}}^{\prime}).$
(2)

Here $\sigma(\boldsymbol{r})\equiv
n_{+}(\boldsymbol{r})-n_{-}(\boldsymbol{r})$ is the vorticity field,

$\displaystyle
n_{\pm}(\boldsymbol{r})\equiv\frac{1}{N_{\pm}}\sum_{i}\delta(\boldsymbol{r}-\boldsymbol{r}^{\pm}_{i})$
(3)

is the local density of positive (negative) vortices, and
$\boldsymbol{r}^{\pm}_{i}$ is the position of the vortex $i$ with circulation
$\pm 1/N$. The vortex densities $n_{\pm}$ satisfy the normalization condition

$\displaystyle\int{\rm d}^{2}\boldsymbol{r}\,n_{\pm}=1.$ (4)

The Green’s function $\phi(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ satisfies
$\nabla^{2}\phi(\boldsymbol{r}-\boldsymbol{r}^{\prime})=-4\pi\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})$.
Here $\phi(\boldsymbol{r}-\boldsymbol{r}^{\prime})=0$ on the boundary
($|\boldsymbol{r}|=1$), and
$\phi(\boldsymbol{r}-\boldsymbol{r}^{\prime})\sim-2\mathop{log}\nolimits|\boldsymbol{r}-\boldsymbol{r}^{\prime}|$
as $|\boldsymbol{r}-\boldsymbol{r}^{\prime}|\rightarrow 0$ Lin (1941). The
stream function

$\displaystyle\psi(\boldsymbol{r})\equiv\int{\rm
d}^{2}\boldsymbol{r}^{\prime}\phi(\boldsymbol{r},\boldsymbol{r}^{\prime})\sigma(\boldsymbol{r}^{\prime}),$
(5)

satisfies the Poisson equation

$\displaystyle\nabla^{2}\psi=-4\pi\sigma(\boldsymbol{r})$ (6)

with the boundary condition $\psi(r=1,\theta)=C$. Here $C$ is a constant.
Recall that the radial velocity

$\displaystyle u_{r}=\frac{1}{r}\frac{\partial\psi}{\partial\theta}.$ (7)

This boundary condition ensures that there is no flow across the boundary of
the domain. Without losing generality, we choose $C=0$, which is equivalent to
including image terms in Eq. (1).

For a rotationally symmetric domain, energy

$\displaystyle E=\frac{1}{2}\int{\rm d}^{2}{\boldsymbol{r}}\,\sigma\psi$ (8)

and angular momentum

$\displaystyle L=\int{\rm d}^{2}\boldsymbol{r}\,r^{2}\sigma$ (9)

are conserved quantities.

The most probable density distribution is given by maximizing the entropy
function

$\displaystyle S=-\int{\rm d}^{2}\boldsymbol{r}\,n_{+}\log n_{+}-\int{\rm
d}^{2}\boldsymbol{r}\,n_{-}\log n_{-},$ (10)

at fixed values of $N_{+}$, $N_{-}$, $E$ and $L$ Joyce and Montgomery (1973).
From the variational equation

$\displaystyle\delta S-\beta\delta E-\alpha\delta L-\mu_{+}\delta
N_{+}/N_{+}-\mu_{-}\delta N_{-}/N_{-}=0,$ (11)

we obtain

$\displaystyle
n_{\pm}(\boldsymbol{r})=\exp\left[\mp\beta\psi(\boldsymbol{r})\mp\alpha
r^{2}+\gamma_{\pm}\right],$ (12)

where $\beta$, $\alpha$ and $\mu_{\pm}$ are Lagrange multipliers and
$\gamma_{\pm}=-\mu_{\pm}-1$. The parameters $\beta$,
$\omega\equiv\alpha/\beta$ and $\mu_{\pm}$ have the interpretation of inverse
temperature, rotation frequency and chemical potentials, respectively.

##  III Onset of clustering

In this section we analyze the possible stable large scale coherent structures
described by Eq. (6) and Eq. (12) near the uniform state. Here we generalized
the method which was developed for analyzing chiral vortex matter Smith and
O’Neil (1990), to the neutral case.

Let us start at a solution $n_{\pm}$ of Eq. (12) at energy $E$ and angular
momentum $L$, and consider a nearby solution $n_{\pm}+\delta n_{\pm}$ at
$E+\delta E$ and $L+\delta L$. The corresponding changes of the constraints
are

$\displaystyle 0$ $\displaystyle=\int{\rm d}^{2}{\boldsymbol{r}}\,\delta
n_{+},$ (13) $\displaystyle 0$ $\displaystyle=\int{\rm
d}^{2}{\boldsymbol{r}}\,\delta n_{-},$ (14) $\displaystyle\delta E$
$\displaystyle=\int{\rm
d}^{2}\boldsymbol{r}\,\psi\delta\sigma+\frac{1}{2}\int{\rm
d}^{2}\boldsymbol{r}\,\delta\psi\delta\sigma,$ (15) $\displaystyle\delta L$
$\displaystyle=\int{\rm d}^{2}\boldsymbol{r}\,r^{2}\delta\sigma.$ (16)

To leading order, we obtain

$\displaystyle\delta n_{+}\simeq$ $\displaystyle
n_{+}(-\psi\delta\beta-\beta\delta\psi+\delta\gamma_{+}-r^{2}\delta\alpha),$
$\displaystyle\delta n_{-}\simeq$ $\displaystyle
n_{-}(\psi\delta\beta+\beta\delta\psi+\delta\gamma_{-}+r^{2}\delta\alpha),$
(17)

where $\delta\gamma_{-},\delta\gamma_{+},\delta\beta$ and $\delta\alpha$ are
obtained from the variations of the constraints. Rewriting them explicitly, we
have

${\cal
Q}\delta\boldsymbol{\mu}=-\delta\boldsymbol{T}+\beta{\boldsymbol{V}}\delta\psi,$
(18)

where
$\delta\boldsymbol{\mu}=(\delta\beta,\delta\gamma_{+},\delta\gamma_{-},\delta\alpha)^{\rm
T}$, $\delta\boldsymbol{T}=(0,0,\delta L,\delta E)^{\rm T}$,
$\boldsymbol{V}\delta\psi=-\int{\rm
d}^{2}\boldsymbol{r}\,(n_{+},n_{-},nr^{2},n\psi)^{\rm T}\delta\psi$,

${\cal Q}\equiv\int{\rm
d}^{2}\boldsymbol{r}\left(\begin{array}[]{cccc}n_{+}\psi&-n_{+}&0&n_{+}r^{2}\\\
n_{-}\psi&0&n_{-}&n_{-}r^{2}\\\ n\psi r^{2}&-n_{+}r^{2}&n_{-}r^{2}&nr^{4}\\\
n\psi^{2}&-n_{+}\psi&n_{-}\psi&n\psi r^{2}\end{array}\right),$ (19)

and $n=n_{+}+n_{-}$ is the total density. Variation of the Eq. (6) gives us

$\displaystyle\nabla^{2}\delta\psi$ $\displaystyle=-4\pi\left(\delta
n_{+}-\delta n_{-}\right)$ (20) $\displaystyle=4\pi(\psi n\delta\beta+\beta
n\delta\psi-n_{+}\delta\gamma_{+}+n_{-}\delta\gamma_{-}+nr^{2}\delta\alpha).$

Our aim is to find stable clustered states which emerge from the homogeneous
state $n_{-}=n_{+}=n_{0}=1/\pi$. For the homogeneous state, $\sigma=0$,
$\psi=0$, $\alpha=0$, $L=0$ and $E=0$. Now let us assume that $\delta\alpha$
is the same order as $\delta\psi$ and from Eq. (18) we obtain

$\displaystyle 0$ $\displaystyle=\delta\gamma_{+}+\delta\gamma_{-},$ (21)
$\displaystyle 0$
$\displaystyle=\delta\gamma_{+}-\delta\gamma_{-}-\delta\alpha-2\beta
n_{0}\int{\rm d}^{2}\boldsymbol{r}\,\delta\psi,$ (22) $\displaystyle\delta L$
$\displaystyle=\beta n_{0}\int{\rm
d}^{2}\boldsymbol{r}(1-2r^{2})\delta\psi-\frac{1}{6}\delta\alpha,$ (23)
$\displaystyle\delta E$ $\displaystyle=\frac{1}{2}\int{\rm
d}^{2}\boldsymbol{r}\,\delta\psi\delta\sigma.$ (24)

Let us introduce operator ${\cal L}$:

$\displaystyle{\cal L}\delta\psi$
$\displaystyle\equiv\nabla^{2}\delta\psi-8\pi
n_{0}\bigg{[}\beta\delta\psi-\beta n_{0}\int{\rm
d}^{2}\boldsymbol{r}\delta\psi$ (25) $\displaystyle-3\beta
n_{0}(1-2r^{2})\int{\rm
d}^{2}\boldsymbol{r}(1-2r^{2})\delta\psi+3(1-2r^{2})\delta L\bigg{]}=0.$

Then Eq. (20) becomes a zero mode equation of the operator ${\cal L}$. The
onset of large scale vortex clusters occurs if Eq. (25) has non-zero
solutions. The value of $\beta$ is undefined in the homogeneous phase within
our mean field approach and depends on the mode developing from the uniform
state. Since the operator ${\cal L}$ is defined on a disc with the Dirichlet
boundary condition, it is natural to decompose Eq. (25) in azimuthal Fourier
harmonics ${\psi_{s}}$ which is characterized by the mode number $s$ and
satisfies $\partial^{2}{\psi_{s}}/\partial{\theta^{2}}=-s^{2}\psi_{s}$:

$\displaystyle\delta\psi=\sum_{s}\epsilon f_{s}\psi_{s}(r,\theta),$ (26)

where $\epsilon\ll 1$ is a small amplitude and $f_{s}$ is the mode
coefficient. We denote $\delta L=L_{0}\epsilon$, $\delta E=E_{0}\epsilon^{2}$
and $\delta\alpha=\epsilon\beta\omega$. Then each mode satisfies

$\displaystyle{\cal L}\psi_{s}(r,\theta)=0,$ (27)

where $\psi_{s}(r,\theta)$ satisfies the boundary condition
$\psi_{s}(r=1,\theta)=0$.

We find that

$\displaystyle\psi_{s}(r,\theta)$ $\displaystyle=c_{s}{\rm
J}_{s}(kr)\cos(s\theta)+b_{s}+ar^{2}$ (28)

solves Eq. (27) with

$\displaystyle a=-\omega=-2c_{s}\beta n_{0}\int{\rm d}^{2}\boldsymbol{r}\,{\rm
J}_{s}(kr)\cos(s\theta)$ (29)

and

$\displaystyle\beta=-\frac{k^{2}}{8\pi n_{0}}.$ (30)

Here ${\rm J}_{s}(r)$ is the Bessel function of the first kind. Consistently,

$\displaystyle L_{0}$ $\displaystyle=\beta n_{0}\int{\rm
d}^{2}\boldsymbol{r}\,(1-2r^{2})\left[\psi_{s}(r,\theta)+\omega r^{2}\right],$
(31) $\displaystyle E_{0}$ $\displaystyle=-\frac{1}{8\pi}\int{\rm
d}^{2}\boldsymbol{r}\,\psi_{s}(r,\theta)\nabla^{2}\psi_{s}(r,\theta).$ (32)

For given $L_{0}$ and $E_{0}$, the parameters $c_{s}$, $k$, and $b_{s}$ are
determined by Eqs. (31) and (32) combined with the Dirichlet boundary
condition

$\displaystyle\psi_{s}(r=1,\theta)=c_{s}{\rm
J}_{s}(k)\cos(s\theta)+b_{s}+a=0.$ (33)

The single-valueless of the stream function requires that $s$ has to be an
integer, namely, $s\in\mathbb{Z}$.

For $s\neq 0$, $L_{0}=0$, $a=-\omega=0$, $b_{s}=0$,

$\displaystyle c_{s}^{2}=\frac{16E_{0}}{k\left[k{\rm J}_{s-1}(k){}^{2}-2s{\rm
J}_{s}(k){\rm J}_{s-1}(k)+k{\rm J}_{s}(k){}^{2}\right]},$ (34)

and $k=j_{s,m}$, where $j_{s,m}$ is the $m$th zero of the Bessel function of
the first kind ${\rm J}_{s}(r)$.

For $s=0$,

$\displaystyle a$ $\displaystyle=-2c_{0}\beta n_{0}\int{\rm
d}^{2}\boldsymbol{r}\,{\rm J}_{0}(kr)=\frac{1}{2}c_{0}k{\rm J}_{1}(k),$ (35)
$\displaystyle b_{0}$ $\displaystyle=-a-c_{0}{\rm J}_{0}(k).$ (36)

For given $E_{0}$ and $L_{0}$, $c_{0}$ and $k$ are determined by

$\displaystyle E_{0}$ $\displaystyle=\frac{1}{8}c_{0}^{2}k^{2}\left[{\rm
J}_{0}(k){}^{2}+{\rm J}_{1}(k){}^{2}\right],$ (37) $\displaystyle L_{0}$
$\displaystyle=-\frac{1}{4}c_{0}k{\rm J}_{3}(k).$ (38)

It is useful to introduce

$\displaystyle\Gamma(k)\equiv\frac{(\delta L)^{2}}{\delta
E}=\frac{L^{2}_{0}}{E_{0}}=\frac{{\rm J}^{2}_{3}(k)}{2\left[{\rm
J}_{0}(k){}^{2}+{\rm J}_{1}(k){}^{2}\right]}$ (39)

as a control parameter.

The ratio $\Gamma(k)$ reaches its maximum value at $k=k^{*}$ with
$j_{1,1}<k^{*}<j_{2,1}$ (see Fig.1). For a given $\Gamma_{0}<\Gamma(k^{*})$,
there are infinitely many values of $k_{c}$ such that
$\Gamma(k_{c})=\Gamma_{0}$. Guided by the maximum entropy principle, the
minimal value of $k_{c}$ corresponds to the equilibrium state. For
$k\rightarrow 0$, $E_{0}\rightarrow 0$, $L_{0}\rightarrow 0$ and this mode
describes the uniform state.

Figure 1: $\Gamma(k)$ as a function of $k$. The maximum value of $\Gamma(k)$
is reached at $k=k^{*}$ and $j_{1,1}<k^{*}<j_{2,1}$.

The modes $s\neq 0$ break $\rm SO(2)$ symmetry and the maximum entropy state
for given energy is the clustered vortex dipole state which corresponds to the
$s=1$ mode Yu  _et al._ (2016). This clustered vortex dipole state has been
recently realized in BEC experiments Gauthier  _et al._ (2019). In this paper,
we focus on states related to the $s=0$ mode.

##  IV Axis-symmetric Clustered States

In this section, we present some (asymptotic) exact results on axis-symmetric
neutral vortex clusters. For axis-symmetric states, the boundary condition Eq.
(7) which is imposed by the most relevant physical condition is fulfilled
automatically.

###  IV.1 Gaussian vortex states

Let us firstly consider $\beta\rightarrow 0$. For finite $\omega$, vortex
distributions $n_{\pm}$ must be uniform. However, when
$\omega\rightarrow\infty$ simultaneously such that $\alpha=\omega\beta$ is
finite, non-trivial distributions can occur. In this special limit, the vortex
densities have the profile of Gaussian distribution:

$\displaystyle n_{+}(r)$
$\displaystyle=-\frac{\pi\left[\exp(-\alpha)-1\right]}{\alpha}\exp\left(-\alpha
r^{2}\right),$ (40) $\displaystyle n_{-}(r)$
$\displaystyle=\frac{\pi\left[\exp(\alpha)-1\right]}{\alpha}\exp\left(\alpha
r^{2}\right),$ (41)

where $\alpha\in(-\infty,\infty)$.

The corresponding stream function reads

$\displaystyle\psi(r)$
$\displaystyle=\frac{\exp(\alpha)\text{Ei}\left(-r^{2}\alpha\right)+\text{Ei}\left(r^{2}\alpha\right)-2\left[\exp(\alpha)+1\right]\log
r}{\exp(\alpha)-1}$
$\displaystyle-\frac{\exp(\alpha)\text{Ei}(-\alpha)+\text{Ei}(\alpha)}{\exp(\alpha)-1},$
(42)

where

$\displaystyle\text{Ei}(x)=-\int_{-x}^{\infty}\frac{\exp(-t)}{t}\,dt$ (43)

is the exponential integral function. The stream function satisfies
$\psi(r=1)=0$ and $d\psi/dr|_{r=1}=0$.

The angular momentum is

$\displaystyle
L=\alpha^{-1}\left[2-\alpha\coth\left(\frac{\alpha}{2}\right)\right].$ (44)

It is easy to see that $L\leq 1$. Figure 2 shows typical vortex densities for
different values of $\alpha$, and angular momentum and energy as functions of
$\alpha$. Note that the Gaussian state is available in the chiral vortex
system as well Smith and O’Neil (1990).

Figure 2: Vortex densities for $\alpha=1$ (a) and $\alpha=6$ (b). The angular
momentum and the energy as functions of $\alpha$ are shown in (c) and (d),
respectively.

###  IV.2 Rotation free vortex states

In this section, we consider clustered vortex states for $\omega=0$ and finite
$\beta<0$.

####  IV.2.1 Onset of axis-symmetric clustered states

Closed to the uniform state, the clustered states can be analyzed using the
formalism developed in Sec.III. The polar angle $\theta$-independent zero
modes $s=0$ carry non-zero angular momentum. For the $s=0$ mode, the rotation
free condition $a=-\omega=c_{0}k{\rm J}_{1}(k)/2=0$ [see Eq. (35)] requires
that $k=j_{1,m}$, where $j_{1,m}$ is the $m$th zero of the Bessel function of
the first kind ${\rm J}_{1}(r)$. These modes occur at
$\beta=\beta_{1,m}=-j_{1,m}^{2}/8\pi n_{0}$ and break $\mathbb{Z}_{2}$
symmetry. The $m=1$ mode starts to emerge at
$\beta=\beta_{t}=\beta_{1,1}\simeq-1.835$ and has the highest statistical
weight among the rotation free modes ($\omega=0$):

$\displaystyle\psi_{0}(r)=c_{0}{\rm J}_{0}(j_{1,1}r)+b_{0},$ (45)

where $c_{0}=\pm 2\sqrt{2}E_{0}/|{\rm J}_{0}(j_{1,1})|j_{1,1}$ and
$b_{0}=-c_{0}{\rm J}_{0}(j_{1,1})$. For this mode, $\Gamma(j_{1,1})\sim
0.545$. Since $\Gamma(j_{1,1})>\Gamma(j_{1,m>1})$, $\Gamma(j_{1,1})$ gives the
upper bond of $\Gamma$ for the rotation free modes. The rotation free axis-
symmetric phase emerges from the uniform phase by varying angular momentum and
energy such that $\Gamma=\Gamma(j_{11})$.

####  IV.2.2 High energy configuration

All the rotation free and axis-symmetry states satisfy

$\frac{1}{r}\frac{d}{dr}r\frac{d}{dr}\psi(r)=-4\pi\left[\exp(-\beta\psi(r)+\gamma_{+})-\exp(\beta\psi(r)+\gamma_{-})\right].$
(46)

The most relevant solution of Eq. (46) should be the nonlinear continuation of
the zero mode $\psi_{0}$ and describes the axis-symmetry equilibrium state
with zero rotation frequency.

At large energies, vorticies with opposite sign are well-separated and the
overlap between $n_{+}$ and $n_{-}$ can be neglected. In this limit, exact
results are available. Let us assume that positive vorticies are concentrated
in the center of the disc and negative vorticies are distributed along the
edge of the disc. The density distribution of positive vortices near $r=0$ can
be obtained analytically by neglecting the influence of negative vortices:

$\displaystyle n_{+}(r)=\frac{4A}{(2-\pi\beta Ar^{2})^{2}},$ (47)

with the boundary conditions $\psi(0)=0$ and $\psi^{\prime}(0)=0$ Smith and
O’Neil (1990). Here $A=[\pi(1-\beta/\beta_{\ast})]^{-1}$ is fixed by the
normalization condition of $n_{+}$ and $\beta_{\ast}=-2$. The
supercondensation occurs at $\beta=\beta_{\ast}$, involving point-like
concentration of the positive vortices and the divergence of energy Kraichnan
(1975); Smith and O’Neil (1990).

Near $r=1$, we can neglect the influence of positive vortices and find the
density distribution of negative vortices

$\displaystyle
n_{-}(r)=\frac{2(2/\beta-1)(1-\beta)^{2}r^{-2\beta}}{\pi\beta(r^{-2\beta+2}+1-2/\beta)^{2}},$
(48)

where the boundary conditions are $\psi(1)=0$ and $\psi^{\prime}(1)=0$.

Note that $\psi(0)$ and $\psi(1)$ can be chosen as arbitrary constants and
here we choose them to be zero for convenience. The boundary condition
$\psi^{\prime}(0)=0$ ensures that $n^{\prime}_{+}(0)=0$ and $n_{+}$ has no
singular behavior near $r=0$. Similarly, the boundary condition
$\psi^{\prime}(1)=0$ implies that $n^{\prime}_{-}(1)=0$ and the absence of
singular behavior of $n_{-}$ near $r=1$. As approximations of densities at
large energies, Eqs. (47) and (48) should be evaluated for
$\beta_{\ast}<\beta$. Combining the critical value of $\beta$ at which the
onset of clustering occurs, we obtain the parameter regime for the rotation
free clustered vortex state:

$\displaystyle\beta_{\ast}<\beta<\beta_{t}.$ (49)

Figure. 3 shows the vortex density distributions at high energies.

In the deep clustered sate, positive vortices are concentrated in a small
region and the total energy are contributed dominantly from positive vortices.
So as $\beta\rightarrow\beta_{\ast}$,

$\displaystyle
E\backsimeq-\frac{2}{\beta^{2}}\left[\log\left(1-\frac{\beta}{\beta_{\ast}}\right)-\frac{\beta}{2}\right].$
(50)

At large energies, the angular momentum is

$\displaystyle L$ $\displaystyle=\bigg{|}\int{\rm
d}^{2}\boldsymbol{r}r^{2}[n_{+}(\boldsymbol{r})-n_{-}(\boldsymbol{r})]\bigg{|}$
$\displaystyle\simeq\bigg{|}\int{\rm
d}^{2}\boldsymbol{r}r^{2}\bigg{[}\frac{4A}{\left(2-\pi\beta
Ar^{2}\right)^{2}}-\frac{2(2/\beta-1)(1-\beta)^{2}r^{-2\beta}}{\pi\beta(r^{-2\beta+2}+1-2/\beta)^{2}}\bigg{]}\bigg{|},$
(51)

and as $\beta\rightarrow\beta_{\ast}$, $L\rightarrow L_{\rm{max}}$ with

$\displaystyle L_{\rm max}=\left|\int{\rm
d}^{2}\boldsymbol{r}r^{2}\left\\{\delta(\boldsymbol{r})-18r^{4}[\pi(2+r^{6})^{2}]^{-1}\right\\}\right|\simeq
0.705.$ (52)

Figure 3: Vortex density distributions at high energies. The densities of
positive vortices (a) and negative vortices (b) are evaluated via Eq. (47) and
Eq. (48), respectively.

##  V Exact results for Chiral vortex clusters

As stated in the previous section, Eq. (48) is the exact solution to Eq. (46),
provided that $n_{+}$ is neglected. Hence Eq. (48) provides an exact vortex
density distribution for a chiral system, which is distinct from the well-
known exact distribution. In this section, we make a summary of relevant exact
results and make a comparison between our finding and the known distribution.

For a rotation free ($\omega=0$) and axis-symmetric chiral system, Eq. (46)
becomes

$\displaystyle\frac{1}{r}\frac{d}{dr}r\frac{d}{dr}\psi(r)$
$\displaystyle=-4\pi n(r),$ (53) $\displaystyle n(r)$
$\displaystyle=\exp(-\beta\psi+\gamma).$ (54)

There is a known exact solution to Eqs. (53) and (54), which is Eq. (47):

$\displaystyle n(r)$ $\displaystyle=\frac{2(\beta+2)}{\pi\left(\beta-\beta
r^{2}+2\right)^{2}},$ (55) $\displaystyle\psi(r)$
$\displaystyle=-\frac{2}{\beta}\log\frac{\beta+2}{\beta(1-r^{2})+2},\quad\psi(0)=0,\,\psi^{\prime}(0)=0,$
(56) $\displaystyle\psi(r)$
$\displaystyle=\frac{2}{\beta}\log\left(1+\frac{\beta(1-r^{2})}{2}\right),\quad\psi(1)=0,\,\psi^{\prime}(0)=0.$
(57)

This solution is valid for $\beta>-2$. The corresponding stream function could
be different depending on the boundary conditions. The vortex density Eq. (55)
exhibits distinct behaviors in different parameter regimes. The vortices
accumulate around the edge for $0<\beta$ while for $-2<\beta<0$ the vortices
are center-concentrated (see Fig.4). Note that in some literature, Eq. (53)
does not have the prefactor $4\pi$ and hence the solution looks slightly
different Menon (2008); Patil and Dunkel (2021).

Figure 4: Typical profiles of the vortex density distribution described by Eq.
(55) in two distinct parameter regimes: $0<\beta$ and $-2<\beta<0$.

Distinct from Eq. (55), our finding is Eq. (48):

$\displaystyle n(r)$
$\displaystyle=\frac{2(2/\beta-1)(1-\beta)^{2}r^{-2\beta}}{\pi\beta(r^{-2\beta+2}+1-2/\beta)^{2}},$
(58) $\displaystyle\psi(r)$
$\displaystyle=-\frac{2}{\beta}\log\left[\frac{2(\beta-1)r^{\beta}}{(\beta-2)r^{2\beta}+\beta
r^{2}}\right],$ (59)

with boundary conditions

$\displaystyle\psi(1)=0,\quad\psi^{\prime}(1)=0.$ (60)

The solution Eq. (58) holds for $\beta<1$ and $\beta\neq 0$. If requiring that
$n^{\prime}(r=0)$ is finite, $\beta<-1/2$. For $0<\beta<1$, Eq. (58) shows
center-concentrated distribution and $n(r\rightarrow 0)\rightarrow\infty$. For
$-1/2<\beta<0$, the vortex density is singular at origin, namely
$n^{\prime}(r\rightarrow 0)\rightarrow\infty$. Vortices distribute around the
edge for $\beta<-1/2$. In contrast to the known exact solution Eq. (55), the
distribution Eq. (58) is peaked on the boundary at negative temperature and is
maximized at origin at positive temperatures. Figure 5 shows typical behaviors
of the vortex density in these parameter regimes.

Figure 5: Typical profiles of the vortex density distribution described by Eq.
(58) in three distinct parameter regimes: $0<\beta<1$, $-1/2<\beta<0$ and
$\beta<-1/2$.

##  VI  Conclusions

Axis-symmetric clustered vortex states for a neutral vortex system confined to
a disc are investigated. Combining the perturbation theory near the uniform
state and asymptotic analysis at high energies, we find the parameter regime
for which the rotation free states are supported. At large energies, the
distributions of positive vortices and negative vortices are well-separated
and the edge-concentrated part provides a new exact vortex density
distribution for the corresponding chiral vortex system.

The onset of a non-axisymmetric vortex cluster appears to proceed via a
second-order phase transition Smith (1989); Smith and O’Neil (1990). It would
be interesting to investigate possible non-axisymmetric states for neutral
systems carrying finite angular momentum. Thanks to the recent experimental
advances Gauthier  _et al._ (2019); Johnstone  _et al._ (2019); Reeves  _et
al._ (2022), our work would motivate experimentally investigating axis-
symmetric clustered phases in a homogeneous Bose-Einstein condensate trapped
in cylindrically symmetric potentials. Due to the presence of conservation of
angular momentum, axis-symmetric clustered phases are expected to have longer
life time than the giant vortex dipole state Gauthier  _et al._ (2019).

## Acknowledgement

We acknowledge J. Nian, T. P. Billam, M. T. Reeves and A. S. Bradley for
useful discussions. X.Y. acknowledges support from the National Natural
Science Foundation of China (Grant No. 12175215), the National Key Research
and Development Program of China (Grant No. 2022YFA 1405300) and NSAF (Grant
No. U1930403).

## References

  * Eyink and Sreenivasan (2006) Gregory L. Eyink and Katepalli R. Sreenivasan, “Onsager and the theory of hydrodynamic turbulence,” Rev. Mod. Phys. 78, 87–135 (2006).
  * Kraichnan and Montgomery (1980) Robert H Kraichnan and David Montgomery, “Two-dimensional turbulence,” Reports on Progress in Physics 43, 547 (1980).
  * Tabeling (2002) Patrick Tabeling, “Two-dimensional turbulence: a physicist approach,” Physics Reports 362, 1–62 (2002).
  * Boffetta and Ecke (2012) Guido Boffetta and Robert E. Ecke, “Two-dimensional turbulence,” Annual Review of Fluid Mechanics 44, 427–451 (2012).
  * Onsager (1949) L. Onsager, “Statistical hydrodynamics,” Il Nuovo Cimento (1943-1954) 6, 279–287 (1949).
  * Vasavada  _et al._ (1998) Ashwin R. Vasavada, Andrew P. Ingersoll, Don Banfield, Maureen Bell, Peter J. Gierasch, Michael J.S. Belton, Glenn S. Orton, Kenneth P. Klaasen, Eric DeJong, H.Herbert Breneman, Todd J. Jones, James M. Kaufman, Kari P. Magee, and David A. Senske, “Galileo imaging of jupiter’s atmosphere: The great red spot, equatorial region, and white ovals,” Icarus 135, 265–275 (1998).
  * Gauthier  _et al._ (2019) Guillaume Gauthier, Matthew T. Reeves, Xiaoquan Yu, Ashton S. Bradley, Mark A. Baker, Thomas A. Bell, Halina Rubinsztein-Dunlop, Matthew J. Davis, and Tyler W. Neely, “Giant vortex clusters in a two-dimensional quantum fluid,” Science 364, 1264–1267 (2019).
  * Johnstone  _et al._ (2019) Shaun P. Johnstone, Andrew J. Groszek, Philip T. Starkey, Christopher J. Billington, Tapio P. Simula, and Kristian Helmerson, “Evolution of large-scale flow from turbulence in a two-dimensional superfluid,” Science 364, 1267–1271 (2019).
  * Reeves  _et al._ (2022) Matthew T. Reeves, Kwan Goddard-Lee, Guillaume Gauthier, Oliver R. Stockdale, Hayder Salman, Timothy Edmonds, Xiaoquan Yu, Ashton S. Bradley, Mark Baker, Halina Rubinsztein-Dunlop, Matthew J. Davis, and Tyler W. Neely, “Turbulent relaxation to equilibrium in a two-dimensional quantum vortex gas,” Phys. Rev. X 12, 011031 (2022).
  * Panico  _et al._ (2023) R Panico, P Comaron, M Matuszewski, AS Lanotte, D Trypogeorgos, G Gigli, M De Giorgi, V Ardizzone, D Sanvitto, and D Ballarini, “Onset of vortex clustering and inverse energy cascade in dissipative quantum fluids,” Nature Photonics , 1–6 (2023).
  * Montgomery and Joyce (1974) David Montgomery and Glenn Joyce, “Statistical mechanics of “negative temperature” states,” Physics of Fluids 17 (1974).
  * Edwards and Taylor (1974) Samuel Frederick Edwards and John Bryan Taylor, “Negative temperature states of two-dimensional plasmas and vortex fluids,” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 336, 257–271 (1974).
  * Williamson (1977) J. H. Williamson, “Statistical mechanics of a guiding centre plasma,” Journal of Plasma Physics 17, 85–92 (1977).
  * Pointin and Lundgren (1976) Y. B. Pointin and T. S. Lundgren, “Statistical mechanics of two‐dimensional vortices in a bounded container,” The Physics of Fluids 19, 1459–1470 (1976).
  * Campbell and O’Neil (1991) LJ Campbell and Kevin O’Neil, “Statistics of two-dimensional point vortices and high-energy vortex states,” Journal of Statistical physics 65, 495–529 (1991).
  * Smith (1989) Ralph A. Smith, “Phase-transition behavior in a negative-temperature guiding-center plasma,” Phys. Rev. Lett. 63, 1479–1482 (1989).
  * Smith and O’Neil (1990) Ralph A. Smith and Thomas M. O’Neil, “Nonaxisymmetric thermal equilibria of a cylindrically bounded guiding‐center plasma or discrete vortex system,” Physics of Fluids B: Plasma Physics 2, 2961–2975 (1990).
  * Yatsuyanagi  _et al._ (2005) Yuichi Yatsuyanagi, Yasuhito Kiwamoto, Hiroyuki Tomita, Mitsusada M. Sano, Takeshi Yoshida, and Toshikazu Ebisuzaki, “Dynamics of two-sign point vortices in positive and negative temperature states,” Phys. Rev. Lett. 94, 054502 (2005).
  * Esler and Ashbee (2015) J. G. Esler and T. L. Ashbee, “Universal statistics of point vortex turbulence,” Journal of Fluid Mechanics 779, 275–308 (2015).
  * Billam  _et al._ (2014) T. P. Billam, M. T. Reeves, B. P. Anderson, and A. S. Bradley, “Onsager-kraichnan condensation in decaying two-dimensional quantum turbulence,” Phys. Rev. Lett. 112, 145301 (2014).
  * Siggia and Aref (1981) Eric D Siggia and Hassan Aref, “Point-vortex simulation of the inverse energy cascade in two-dimensional turbulence,” Physics of Fluids (1958-1988) 24, 171–173 (1981).
  * Reeves  _et al._ (2013) Matthew T. Reeves, Thomas P. Billam, Brian P. Anderson, and Ashton S. Bradley, “Inverse energy cascade in forced two-dimensional quantum turbulence,” Phys. Rev. Lett. 110, 104501 (2013).
  * Reeves  _et al._ (2012) M. T. Reeves, B. P. Anderson, and A. S. Bradley, “Classical and quantum regimes of two-dimensional turbulence in trapped bose-einstein condensates,” Phys. Rev. A 86, 053621 (2012).
  * Billam  _et al._ (2015) T. P. Billam, M. T. Reeves, and A. S. Bradley, “Spectral energy transport in two-dimensional quantum vortex dynamics,” Phys. Rev. A 91, 023615 (2015).
  * Simula  _et al._ (2014) Tapio Simula, Matthew J. Davis, and Kristian Helmerson, “Emergence of order from turbulence in an isolated planar superfluid,” Phys. Rev. Lett. 113, 165302 (2014).
  * Gurarie and Chow (2004) D. Gurarie and K. W. Chow, “Vortex arrays for sinh-Poisson equation of two-dimensional fluids: Equilibria and stability,” Physics of Fluids 16, 3296–3305 (2004).
  * Yu  _et al._ (2016) Xiaoquan Yu, Thomas P. Billam, Jun Nian, Matthew T. Reeves, and Ashton S. Bradley, “Theory of the vortex-clustering transition in a confined two-dimensional quantum fluid,” Phys. Rev. A 94, 023602 (2016).
  * Salman and Maestrini (2016) Hayder Salman and Davide Maestrini, “Long-range ordering of topological excitations in a two-dimensional superfluid far from equilibrium,” Phys. Rev. A 94, 043642 (2016).
  * Valani  _et al._ (2018) Rahil N Valani, Andrew J Groszek, and Tapio P Simula, “Einstein–bose condensation of onsager vortices,” New Journal of Physics 20, 053038 (2018).
  * Han and Tsubota (2018) Junsik Han and Makoto Tsubota, “Onsager vortex formation in two-component bose–einstein condensates,” Journal of the Physical Society of Japan 87, 063601 (2018).
  * Han and Tsubota (2019) Junsik Han and Makoto Tsubota, “Phase separation of quantized vortices in two-component miscible bose-einstein condensates in a two-dimensional box potential,” Phys. Rev. A 99, 033607 (2019).
  * Skaugen and Angheluta (2016) Audun Skaugen and Luiza Angheluta, “Vortex clustering and universal scaling laws in two-dimensional quantum turbulence,” Phys. Rev. E 93, 032106 (2016).
  * Skaugen and Angheluta (2017) Audun Skaugen and Luiza Angheluta, “Origin of the inverse energy cascade in two-dimensional quantum turbulence,” Phys. Rev. E 95, 052144 (2017).
  * Kanai and Guo (2021) Toshiaki Kanai and Wei Guo, “True mechanism of spontaneous order from turbulence in two-dimensional superfluid manifolds,” Phys. Rev. Lett. 127, 095301 (2021).
  * Patil and Dunkel (2021) Vishal P. Patil and Jörn Dunkel, “Chiral edge modes in helmholtz-onsager vortex systems,” Phys. Rev. Fluids 6, 064702 (2021).
  * Pandit  _et al._ (2017) Rahul Pandit, Debarghya Banerjee, Akshay Bhatnagar, Marc Brachet, Anupam Gupta, Dhrubaditya Mitra, Nairita Pal, Prasad Perlekar, Samriddhi Sankar Ray, Vishwanath Shukla, and Dario Vincenzi, “An overview of the statistical properties of two-dimensional turbulence in fluids with particles, conducting fluids, fluids with polymer additives, binary-fluid mixtures, and superfluids,” Physics of Fluids 29 (2017), 10.1063/1.4986802.
  * Joyce and Montgomery (1973) Glenn Joyce and David Montgomery, “Negative temperature states for the two-dimensional guiding-centre plasma,” Journal of Plasma Physics 10, 107–121 (1973).
  * Ambegaokar  _et al._ (1980) Vinay Ambegaokar, B. I. Halperin, David R. Nelson, and Eric D. Siggia, “Dynamics of superfluid films,” Phys. Rev. B 21, 1806–1826 (1980).
  * Aref  _et al._ (1999) H Aref, PL Boyland, MA Stremler, and DL Vainchtein, “Turbulent statistical dynamics of a system of point vortices,” in  _Fundamental Problematic Issues in Turbulence_ (Springer, 1999) pp. 151–161.
  * Newton (2001) Paul K Newton,  _The N-vortex problem: analytical techniques_ , Vol. 145 (Springer Science & Business Media, 2001).
  * Eyink and Spohn (1993) GL Eyink and H Spohn, “Negative-temperature states and large-scale, long-lived vortices in two-dimensional turbulence,” Journal of statistical physics 70, 833–886 (1993).
  * Caglioti  _et al._ (1995) Emanuele Caglioti, Pierre-Louis Lions, Carlo Marchioro, and M Pulvirenti, “A special class of stationary flows for two-dimensional euler equations: a statistical mechanics description. part ii,” Communications in Mathematical Physics 174, 229–260 (1995).
  * Lin (1941) C. C. Lin, “On the motion of vortices in two dimensions,” Proceedings of the National Academy of Sciences 27, 570–575 (1941).
  * Kraichnan (1975) Robert H. Kraichnan, “Statistical dynamics of two-dimensional flow,” Journal of Fluid Mechanics 67, 155–175 (1975).
  * Menon (2008) G. Menon, “Statistical theories of turbulence,”  (2008), http://www.dam.brown.edu/people/menon/ publications/notes/turb.pdf.


# MultiSiam: A Multiple Input Siamese Network For Social Media Text
Classification And Duplicate Text Detection

Sudhanshu Bhoi1, Swapnil Markhedkar2, Shruti Phadke3, and Prashant Agrawal4
Department of Computer Engineering
Pune Institute of Computer Technology,
Pune, India  Email: { 1sudh.bhoi.30, 2swapnilmarkhedkar, 3shrutidphadke,
4prashantsa255<EMAIL_ADDRESS>
###### Abstract

Social media accounts post increasingly similar content, creating a chaotic
experience across platforms, which makes accessing desired information
difficult. These posts can be organized by categorizing and grouping
duplicates across social handles and accounts. There can be more than one
duplicate of a post, however, a conventional Siamese neural network only
considers a pair of inputs for duplicate text detection. In this paper, we
first propose a multiple-input Siamese network, MultiSiam. This condensed
network is then used to propose another model, SMCD (Social Media
Classification and Duplication Model) to perform both duplicate text grouping
and categorization. The MultiSiam network, just like the Siamese, can be used
in multiple applications by changing the sub-network appropriately.

##  I Introduction

A significant contributor to the millions of bytes of data produced daily are
the myriad of social media networks in use today [1]. A sizable proportion of
these social media posts are comparable, or duplicates. This, combined with
the increasingly popular infinite feed format makes accessing desirable
information difficult. We aim to optimize the user’s feed in two ways: first,
by dividing posts into different categories based on their content, and
second, by aggregating duplicate posts on and across different social media
into one post.

The two natural language processing problems, categorization, and duplicate
detection involve some common steps. Duplicate posts will belong to the same
category, which means word embeddings calculated for categorization can be
used for duplicate detection. Developing a common model which can handle both
the tasks and understand the relation between them, will avoid repetitive
calculations and improve the application pipeline.

For a user, more than 2 duplicate posts can be found across different handles
or various social media accounts. This requires more than a pair of posts to
be compared. Conventional Siamese models only consider a pair of inputs, which
means, to compare groups of inputs, pairs need to be prepared through careful
batching to ensure no two pairs from the same batch belong to the same group.
Moreover, making all possible pairs from a group makes learning redundant.

In this paper, we propose a modified Siamese neural network, MultiSiam, with a
condensed architecture for duplicate detection, which can handle multiple
inputs at a time. We further use this to create a model, SMCD (Social Media
Classification and Duplication Model), for the proposed application of social
media text classification and duplicate text detection

##  II Background and Related Work

The Siamese network, first introduced in [2], has been successfully
implemented across computer vision [3] and natural language processing [4]
applications. Various distance measures, like Euclidean distance [5], Jaccard
coefficient [6], Cosine similarity [7] have been used to find similarity
between 2 outputs of the sub-networks. This paper [8] uses Cosine similarity
and the widely used Triplet loss. Other losses include Contrastive loss [9],
one of the oldest, which considers the distance between 2 inputs, while
Quadruplet loss [10] considers 4 inputs. In 2020, [11] introduced Circle loss
which aims to maximize within-class similarity and minimize between-class
similarity.

Text categorization is a classical problem in NLP. [12] reviews more than 150
deep learning models including those based on RNN, CNN, and the attention
mechanism. For duplicate question-pair detection, the Siamese network has out-
performed multi-layer perceptron architecture and the traditional shingling
algorithm [13]. Siamese can also be used for text categorization [14]. It is
especially helpful when the number of categories is large, not known during
training and examples per category are small [15].

Social media related tasks can be improved by considering different
modalities, text and image [16], hyperlinks [17], and hashtags [18] from a
post.

##  III Network Architecture

###  III-A Multiple Input Siamese Network (MultiSiam)

A Siamese network has 2 input fields with 2 separate sub-networks to act on
them [2]. Many applications like signature verification [2] or duplicate
question detection [19] require the sub-networks to be identical, having the
same parameters and weights. This paper proposes a network that can accept 2
or more input fields with a sub-network to act on them.

Figure 1: The MultiSiam network architecture with example sub-networks.

The MultiSiam network, shown in Figure 1, expects a batch of groups of input
fields. In case of duplicate detection, a group contains all texts that are
duplicates of each other. Each batch should have groups of equal length. Each
group should have texts of equal length. Text can be padded with zeros or a
special token from vocabulary, while a group can be padded by sampling a
required number of texts from the same group. Thus the input has dimension of
$batch\\_size$ x $group\\_size$ x $text\\_size$.

This input is then reshaped to $(batch\\_size*group\\_size)$ x $text\\_size$
for training. It is then passed through a sub-network. In case of duplicate
detection, the sub-network consists of a Word Embedding layer and an LSTM
layer [19]. Note that this can be replaced with any subnetwork based on the
application as shown in Figure 1.

The output has dimension of $(batch\\_size*group\\_size)$ x
$embedding\\_size$. This is finally reshaped into $batch\\_size$ x
$group\\_size$ x $embedding\\_size$. Thus the network produces embeddings for
each of the input fields. These embeddings are then used for calculating the
triplet loss [8].

####  III-A1 Triplet Loss

The triplet loss aims to minimize the distance between the anchor and positive
and maximizes the distance between the anchor and negative [8]. From the
network output, if one text is considered as an anchor then all the other
texts in the same group are positives while those belonging to other groups
are negatives. This can be represented as a matrix as explained in [20] for
$group\\_size$ of 2. Here we generalize it to dimension $batch\\_size_{1}$ x
$batch\\_size_{2}$ x $...$ x $batch\\_size_{group\\_size}$.

Each cell in the matrix represents the cosine similarity of the corresponding
embeddings. The diagonal of such a matrix represents the cosine similarity
with positives, while off-diagonal represents that with negatives. Considering
the first in each group to be the anchor, the mean of all the negatives
corresponding to it and the highest off-diagonal value is calculated,
resulting in 2 $batch\\_size$ vectors $mean\\_neg$ and $closest\\_neg$
respectively [20].

Triplet Loss is then calculated as follows:

$dist_{i}^{ap}=cos(f(x_{i}^{a}),f(x_{i}^{p}))$ (1)
$dist_{ij}^{an}=cos(f(x_{i}^{a}),f(x_{ij}^{n}))_{j\neq i}$ (2)
$mean\\_neg_{i}=(\sum_{j}^{N-1}dist_{ij}^{an})_{j\neq i}/(N-1)$ (3)
$closest\\_neg_{i}=max_{j}(dist_{ij}^{an})_{j\neq i}$ (4)
$cost_{i}^{1}=max(-dist_{i}^{ap}+mean\\_neg_{i}+\alpha,0)$ (5)
$cost_{i}^{2}=max(-dist_{i}^{ap}+closest\\_neg_{i}+\alpha,0)$ (6)
$cost_{i}=cost_{i}^{1}+cost_{i}^{2}$ (7) $TripletLoss=\sum_{i}^{N}cost_{i}$
(8)

Where $x^{a},x^{p},x^{n}$ are anchor, positive, and negative inputs
respectively. $f(x)$ is the embedding from the network. $\alpha$ is the margin
between positive and negative pairs. $N$ is the $batch\\_size$.

####  III-A2 Inference

Inference can be of two types. One, groups need to be formed from given
inputs, shown in Figure 2. In this, the Reshape layer is disabled so the
MultiSiam produces embeddings for the inputs. This is similar to using only
one of the sub-networks in Siamese during evaluation. Distance measures like
cosine similarity and a threshold value can be used between the embeddings to
group inputs together. Two, already formed groups need to be verified. In
this, MultiSiam outputs groups of embeddings as in Figure 1. A distance
measure between embeddings of a group and a threshold value then can determine
whether the group constitutes duplicates.

Figure 2: MultiSiam, during inference, with distance measure and example
groups of given inputs based on the embeddings produced.

###  III-B Social Media Classification and Duplication Model (SMCD)

In the proposed application, duplicate texts are likely to be found in a
single category. So a single embedding can be used for categorization as well
as further processed for duplicate detection. This model leverages the above
fact to combine the two tasks using the concise MultiSiam network.

First, the input is reshaped according to the MultiSiam network. Then it is
passed to the Word Embeddings and LSTM layers to produce embeddings E. These
embeddings, E, are then passed through an LSTM layer, a Dense layer, and a
Softmax function to get a category. The same embeddings, E, are also passed to
another LSTM layer [21] and then reshaped (as in the MultiSiam) to get the
embeddings for duplicate detection. Refer Figure 3.

Cross-entropy loss is calculated on predicted and actual categories of
dimension $(batch\\_size*group\\_size)$ x $num\\_categories$. Triplet loss is
calculated on the final embeddings as mentioned in the MultiSiam above. Thus,
for a given batch of texts, this proposed model finds their duplicate
embeddings as well as their categories. A distance measure, like cosine
similarity, can be used between these embeddings to get groups of duplicate
texts as shown in Figure 2.

Figure 3: SMCD model for text categorization and duplicate text detection with
respective losses.

##  IV Experiments and Results

The experiments were performed using TensorFlow v2.5 on a 2.3 GHz Quad-Core
Intel Core i5 CPU. These experiments aimed to show the similarity between
proposed MultiSiam and existing Siamese for a pair of inputs under the same
conditions. They also demonstrate the model’s usefulness with more than 2
inputs in a group and integration with the categorization task in the SMCD
model.

###  IV-A MultiSiam

We trained a MultiSiam and a Siamese model on the Quora duplicate questions
dataset [22]. The dataset consists of 404348 pairs of questions with a
corresponding label indicating whether they are duplicates. Out of these
149306 are duplicate pairs and were considered for training with a 90:10
split. Further 14931 non-duplicates pairs were added to the test set. Each
model was then trained for 10 epochs with batches of 256 examples. Triplet
Loss had a margin of 0.25 and a pair was considered duplicate if its cosine
similarity was above 0.7. Although MultiSiam is a little slower due to
dimensional operations, it performs very similarly to the Siamese model in
terms of accuracy. Refer Table I.

TABLE I: Model Comparison Model | Accuracy | CPU Time
---|---|---
Siamese | 0.74 | 1 s/step
MultiSiam | 0.72 | 3 s/step

For more than 2 inputs, experiments were performed on a custom dataset in the
SMCD model.

###  IV-B SMCD

We selected 13 different categories to classify posts into, namely: Business,
Politics, Celebrity, Entertainment, Finance, Health & Wellness, Motivation,
Sports, News, Promotions, Gaming, Travel, and Technology.

TABLE II: Sample from the custom dataset prepared to train the SMCD model before any processing Text | Category | Group ID
---|---|---
Person A officially takes over as Company A CEO from Person B | Business | 1
Today is Day 1 for Person A as former Company B CEO moves to Company A corner office URL by MENTION | Business | 1
Person A takes over the role of Company A CEO from Person B as the company faces an antitrust push. MENTION reports. | Business | 1
Company B Software A beta is here and ready for you to use | Tech | 2
Software A public beta hands-on: These 5 features make all the difference | Tech | 2
Here’s what’s new in Software A beta 2 | Tech | 2
The best new features and changes in Software A beta 2 | Tech | 2

To create the dataset, we sourced posts from 10-12 social media handles per
category (from Reddit, Facebook, and Twitter) which were determined to
exclusively post content belonging to that particular category, and labeled
them accordingly. The labeling errors, if any, were corrected manually.
Duplicate posts were sourced by querying multiple handles for the same content
topic, both within and across the three social media networks.

The dataset thus collected had around 15 groups per category. Each group
consisted of, on average, 2 to 4 texts from the posts. Table II shows a sample
of the dataset (anonymized for this paper). As mentioned in section III-A,
these groups were padded or truncated to have a group size of 4. Thus the
dataset had 201 groups in total and 804 texts after processing. Grouped texts
were used to train the duplicate-detection-half of the model while
corresponding grouped category labels were used to train the categorization-
half. The model was trained for 20 epochs with a batch size of 16.

During the evaluation, we provided a list of texts for which we needed
duplicate embeddings and the category the inputs belonged to. These embeddings
were then used to group the inputs. Table III demonstrates an output from the
model.

TABLE III: Duplicate Grouped and Categorized texts from the SMCD model output obtained by providing list of texts as input Group ID |  Texts | Categories
---|---|---
1 |  Former President announces law suit against Company A | News
Lawsuit announced against Company A by President Person A | News
Former President plans to sue Company A | News
2 |  TV Show A predicted Pandemic A outbreak in 1993? | Entertainment
Fans of ’TV Show A’ convinced show predicted the pandemic | Entertainment

###  IV-C Application

The application [23], for every user, aggregates their feed (across different
social media networks) into one optimized feed, with posts separated into
different categories, and duplicate posts identified and clubbed into one.
This is done using the SMCD model proposed above.

##  V Conclusion and Future Scope

The accuracy in section IV can be improved by various techniques like transfer
learning, text pre-processing, hyper-parameter tuning, etc. Further, a
standard multi-input duplicate dataset can be created for various tasks to
evaluate the models. Different loss functions, sub-networks can be also
experimented with in future work.

Currently, with the help of the proposed MultiSiam network, we were able to
create a condensed model, SMCD, that could categorize texts and group
duplicate texts together. In the future, with NLP techniques like text
summarization, social media feeds can be made more organized. We believe such
a system will help improve users’ mental health as well as increase their
productivity.

The MultiSiam network can also be used for computer vision, natural language
processing, and various other tasks by using an appropriate sub-network. It
brings all the goodness of Siamese networks to more than 2 inputs at a time.

## Acknowledgment

The authors would like to thank Prof. R. V. Bidwe and Vedang Mandhana for
their continued guidance and support throughout the research and development
process.

## References

  * [1] Domo.com (2021) Becoming A Data-Driven CEO, https://www.domo.com/solution/data-never-sleeps-8
  * [2] Bromley, J., Guyon, I., LeCun, Y., Säckinger, E., & Shah, R. (1993). Signature verification using a ”siamese” time delay neural network. Advances in neural information processing systems, 6, 737-744.
  * [3] Pang, H., Xuan, Q., Xie, M., Liu, C., & Li, Z. (2021). Research on Target Tracking Algorithm Based on Siamese Neural Network. Mobile Information Systems, 2021.
  * [4] Wang, Z., Zhang, X., & Tan, Y. (2021). Chinese Sentences Similarity via Cross-Attention Based Siamese Network. arXiv preprint arXiv:2104.08787.
  * [5] Vinayakumar, R., & Soman, K. P. (2020). Siamese neural network architecture for homoglyph attacks detection. ICT Express, 6(1), 16-19.
  * [6] Wang, Z., Hamza, W., & Florian, R. (2017). Bilateral multi-perspective matching for natural language sentences. arXiv preprint arXiv:1702.03814.
  * [7] Xia, P., Zhang, L., & Li, F. (2015). Learning similarity with cosine similarity ensemble. Information Sciences, 307, 39-52.
  * [8] Schroff, F., Kalenichenko, D., & Philbin, J. (2015). Facenet: A unified embedding for face recognition and clustering. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 815-823).
  * [9] Chen, T., Kornblith, S., Norouzi, M., & Hinton, G. (2020, November). A simple framework for contrastive learning of visual representations. In International conference on machine learning (pp. 1597-1607). PMLR.
  * [10] Chen, W., Chen, X., Zhang, J., & Huang, K. (2017). Beyond triplet loss: a deep quadruplet network for person re-identification. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 403-412).
  * [11] Sun, Y., Cheng, C., Zhang, Y., Zhang, C., Zheng, L., Wang, Z., & Wei, Y. (2020). Circle loss: A unified perspective of pair similarity optimization. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 6398-6407).
  * [12] Minaee, S., Kalchbrenner, N., Cambria, E., Nikzad, N., Chenaghlu, M., & Gao, J. (2021). Deep Learning–based Text Classification: A Comprehensive Review. ACM Computing Surveys (CSUR), 54(3), 1-40.
  * [13] Addair, T. (2017). Duplicate question pair detection with deep learning. Stanf. Univ. J.
  * [14] C. Shih, B. Yan, S. Liu and B. Chen, ”Investigating Siamese LSTM networks for text categorization,” 2017 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), 2017, pp. 641-646, doi: 10.1109/APSIPA.2017.8282104.
  * [15] Chopra, S., Hadsell, R., & LeCun, Y. (2005, June). Learning a similarity metric discriminatively, with application to face verification. In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05) (Vol. 1, pp. 539-546). IEEE.
  * [16] Duong, C. T., Lebret, R., & Aberer, K. (2017). Multimodal classification for analysing social media. arXiv preprint arXiv:1708.02099.
  * [17] Kinsella, S., Wang, M., Breslin, J. G., & Hayes, C. (2011, May). Improving categorisation in social media using hyperlinks to structured data sources. In Extended Semantic Web Conference (pp. 390-404). Springer, Berlin, Heidelberg.
  * [18] Stavrianou, A., Brun, C., Silander, T., & Roux, C. (2014). NLP-based feature extraction for automated tweet classification. Interactions between Data Mining and Natural Language Processing, 145.
  * [19] Shrestha, M., Dhungana, M., Panta, S., Pudasaini, U., & Giri, S. Duplicate Question Detection (DQD) using Siamese LSTM. NCE JOURNAL OF SCIENCE AND ENGINEERING (NJSE), 29.
  * [20] Mourri, Kaiser, Shyu. “Siamese Networks”. NLP with Sequence Models. deeplearning.ai. Coursera. https://www.coursera.org/learn/sequence-models-in-nlp
  * [21] Hochreiter, S., & Schmidhuber, J. (1997). Long short-term memory. Neural computation, 9(8), 1735-1780.
  * [22] Shankar, I., Nikhil, D., & Kornél, C. (2017) First Quora Dataset Release: Question Pairs. Quora. https://quoradata.quora.com/First-Quora-Dataset-Release-Question-Pairs
  * [23] Prashant, A., Shruti, P., Swapnil, M., Sudhanshu, B. (2021) Social-Media-Aggregator [Source Code]. https://github.com/prashant255/Social-Media-Aggregator