# X-Ray Scattering from Light-Driven Spin Fluctuations in a Doped Mott
Insulator

Yao Wang  Department of Physics and Astronomy, Clemson University, Clemson,
South Carolina 29631, USA  Yuan Chen  Department of Applied Physics, Stanford
University, Stanford, California 94305, USA  Stanford Institute for Materials
and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill
Road, Menlo Park, California 94025, USA.  Thomas P. Devereaux  Stanford
Institute for Materials and Energy Sciences, SLAC National Accelerator
Laboratory, 2575 Sand Hill Road, Menlo Park, California 94025, USA.
Department of Materials Science and Engineering, Stanford University,
Stanford, California 94305, USA.  Geballe Laboratory for Advanced Materials,
Stanford University, Stanford, California 94305, USA.  Brian Moritz  Stanford
Institute for Materials and Energy Sciences, SLAC National Accelerator
Laboratory, 2575 Sand Hill Road, Menlo Park, California 94025, USA.  Matteo
Mitrano  Department of Physics, Harvard University, Cambridge, Massachusetts
02138, USA

###### Abstract

ABSTRACT
Manipulating spin fluctuations with ultrafast laser pulses is a promising
route to dynamically control collective phenomena in strongly correlated
materials. However, understanding how photoexcited spin degrees of freedom
evolve at a microscopic level requires a momentum- and energy-resolved
characterization of their nonequilibrium dynamics. Here, we study the
photoinduced dynamics of finite-momentum spin excitations in two-dimensional
Mott insulators on a square lattice. By calculating the time-resolved resonant
inelastic x-ray scattering cross-section, we show that an ultrafast pump above
the Mott gap induces a prompt softening of the spin excitation energy,
compatible with a transient renormalization of the exchange interaction. While
spin fluctuations in a hole-doped system (paramagnons) are well described by
Floquet theory, magnons at half filling are found to deviate from this
picture. Furthermore, we show that the paramagnon softening is accompanied by
an ultrafast suppression of $d$-wave pairing correlations, indicating a link
between the transient spin excitation dynamics and superconducting pairing far
from equilibrium.

##  I Introduction

Following the early demonstration of ultrafast demagnetization in ferromagnets
Beaurepaire et al. (1996), the optical manipulation of spin degrees of freedom
emerged as an effective strategy for steering electronic properties in quantum
materials Kirilyuk et al. (2010). Over the years, light excitation protocols
evolved from purely thermal effects to more sophisticated optical excitations
involving Zeeman interaction Kampfrath et al. (2011), nonlinear phonon
excitation Nova et al. (2016); Disa et al. (2020), or the renormalization of
local interactions Mikhaylovskiy et al. (2015); Ron et al. (2020). While most
of these experimental efforts aim to manipulate magnetic orders, driving
fluctuating spins with light could represent an effective dynamical control
strategy for superconductivity in strongly-correlated electron systems.

In the copper oxides, experimental Rossat-Mignod et al. (1991); Mook et al.
(1993); Fong et al. (1999) and theoretical works Scalapino et al. (1986); Gros
et al. (1987); Kotliar and Liu (1988); Schrieffer et al. (1989); Scalapino
(1995); Tsuei and Kirtley (2000); Maier et al. (2006) suggest that spin
fluctuations around the antiferromagnetic wavevector
$\textbf{Q}_{\mathrm{AF}}=(\pi,\pi)$ could contribute to the pairing
interaction and determine the superconducting critical temperature $T_{c}$.
Whether optically driving these spin fluctuations may explain the recent
observation of light-induced superconductivity in certain copper oxides Fausti
et al. (2011); Hu et al. (2014); Kaiser et al. (2014); Nicoletti et al. (2014)
or be conducive to further optimized light-induced coherence is an open
question.

Figure 1:  Probing light-driven spin fluctuations with resonant inelastic
x-ray scattering.  a Sketch of a pump-probe resonant inelastic x-ray
scattering (RIXS) experiment with $\pi$-$\sigma$ polarization. b Magnon
dispersion in the 2D Hubbard model in (blue solid line) and out of equilibrium
(red dashed line, pump field amplitude $A_{0}=0.6$) as determined through
linear spin wave and Floquet linear spin wave (FLSW) theory [see Eqs. (7)].
The gray shaded region indicates momenta which are kinematically inaccessible
to Cu $L$-edge RIXS.

Directly probing the effect of ultrafast light pulses on spin fluctuations at
finite momentum $\mathbf{q}$ requires a tool beyond optical probes, which are
only sensitive to long-wavelength ($\mathbf{q}\sim 0$) dynamics Zhao et al.
(2011); Dal Conte et al. (2012); Batignani et al. (2015); Yang et al. (2020).
Resonant inelastic x-ray scattering (RIXS) at transition-metal $L$-edge RIXS
with cross-polarized incident and scattered x-rays (see Fig. 1a) is sensitive
to magnetic excitations, and is able to map their dispersion thanks to the
high photon momentum Ament et al. (2009); Haverkort (2010); Ament et al.
(2011). Pioneering time-resolved resonant inelastic x-ray scattering (trRIXS)
experiments in the hard x-ray regime provided the first snapshots of
pseudospin dynamics near $\textbf{Q}_{\mathrm{AF}}$ in layered iridates Dean
et al. (2016); Mazzone et al. (2021) and future soft X-ray experiments at the
Cu and Ni $L$-edges will further expand the reach of this new ultrafast
technique. Based on these developments, it becomes critical to theoretically,
and possibly experimentally, investigate how dispersive spin fluctuations and
superconducting pairing evolve in a prototypical correlated system driven far
from equilibrium.

In this work, we present a trRIXS calculation of finite-momentum spin
fluctuations in the paradigmatic single-band Hubbard model. We calculate the
full trRIXS cross section for a direct $L_{3}$-edge absorption for
$\pi$-$\sigma$ polarizations and for two different hole doping levels. We show
that an ultrafast resonant pump induces a prompt softening of the short-ranged
spin excitation spectrum, compatible with a light-induced renormalization of
the exchange interaction [sketched in Fig. 1b], and a simultaneous suppression
of $d$-wave pairing correlations. Our findings in a canonical Hubbard model
demonstrate that light pulses can be used to tailor finite-momentum spin
fluctuations and pave the way to novel strategies for the light control of
unconventional superconductors.

##  II Results and Discussion

###  II.1 Microscopic model

With a particular focus on the high-$T_{c}$ cuprates, we consider valence and
conduction electrons described by the 2D single-band Hubbard model and for two
different hole dopings ($x=0$ and $x=0.167$). Since the x-ray scattering
explicitly involves an intermediate core level, the full Hamiltonian reads as

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle-\sum_{\textbf{i},\textbf{j},\sigma}\\!t_{h}^{\textbf{i}\textbf{j}}(d_{\textbf{j}\sigma}^{\dagger}d_{\textbf{i}\sigma}\\!+\\!h.c.)\\!+U\sum_{i}n^{(d)}_{\textbf{i}\uparrow}n^{(d)}_{\textbf{i}\downarrow}$
(1) $\displaystyle+\sum_{\textbf{i}\alpha\sigma}E_{\rm
edge}(1-n^{(p)}_{\textbf{i}\alpha\sigma})-U_{c}\sum_{\textbf{i}\alpha\sigma\sigma^{\prime}}n^{(d)}_{\textbf{i}\sigma}(1-n^{(p)}_{\textbf{i}\alpha\sigma^{\prime}})$
$\displaystyle+\lambda\sum_{i\alpha\alpha^{\prime}\atop\sigma\sigma^{\prime}}p_{\textbf{i}\alpha\sigma}^{\dagger}\chi_{\alpha\alpha^{\prime}}^{\sigma\sigma^{\prime}}p_{\textbf{i}\alpha^{\prime}\sigma^{\prime}}.$

Here, $d_{\textbf{i}\sigma}$ ($d_{\textbf{i}\sigma}^{\dagger}$) annihilates
(creates) a $3d_{x^{2}-y^{2}}$ electron at site i with spin $\sigma$, while
$p_{\textbf{i}\alpha\sigma}$ ($p_{\textbf{i}\alpha\sigma}^{\dagger}$)
annihilates (creates) a $2p_{\alpha}$ electron ($\alpha\\!=\\!x,y,z$).
$t_{h}^{\textbf{i}\textbf{j}}$ is the hopping integral between site $i$ and
$j$,
$n^{(d)}_{\textbf{i}}\\!=\\!\sum_{\sigma}d_{\textbf{i}\sigma}^{\dagger}d_{\textbf{i}\sigma}$
and
$n^{(p)}_{\textbf{i}\alpha\sigma}\\!=\\!p_{\textbf{i}\alpha\sigma}^{\dagger}p_{\textbf{i}\alpha\sigma}$
are the valence and core-level electron density operators, respectively. The
nearest-neighbor hopping amplitude is set to $t_{h}=300$ meV, the next-nearest
neighbor hopping to $t^{\prime}_{h}\\!=\\!-0.3\,t_{h}$, and the on-site
repulsion to $U\\!=\\!8t_{h}$ Jia et al. (2014, 2016). The core-hole potential
$U_{c}$ is instead fixed at 4$t_{h}$ and regarded identical for all 2$p$
orbitals Tsutsui et al. (2000); Jia et al. (2014, 2016). $E_{\rm edge}=938$ eV
represents instead the Cu $L$-edge absorption energy, i.e., the energy
difference between the $3d$ and $2p$ orbitals without spin-orbit coupling.
Finally, the spin-orbit coupling with the x-ray-induced core hole in the
degenerate $2p$-orbitals is accounted by the last term in Eq. (1) with
$\lambda=13$ eV Tsutsui et al. (2000); Kourtis et al. (2012).

Figure 2:  Snapshots of the time-resolved resonant inelastic x-ray scattering
signal. Cross-polarized time-resolved resonant inelastic x-ray scattering
(trRIXS) spectra for selected pump-probe time delays as function of the
incident photon energy ${\omega_{\rm in}}$ for a $x=0$ and b $x=0.167$ hole
doping. $\omega={\omega_{\rm in}}-{\omega_{\rm s}}$ indicates the energy loss,
while the momentum transfer is kept fixed at $\mathbf{q}=(\pi/2,\pi/2)$. The
peak at $\omega\sim 0.3$ eV corresponds to the collective spin fluctuations of
the system. The red arrow marks the center of mass of the resonance along the
${\omega_{\rm in}}$ axis. The color bar indicates the spectral intensity in
arbitrary units (a.u.).

The trRIXS cross-section at a direct absorption edge can be written as Chen et
al. (2019); Wang et al. (2020) (see Supplementary note 1 for further details)

$\displaystyle\mathcal{I}(\mathbf{q},{\omega_{\rm s}},{\omega_{\rm in}},t)\\!$
$\displaystyle\mkern-6.0mu=$ $\displaystyle\\!\mkern-12.0mu\iiiint
e^{i{\omega_{\rm in}}(t_{2}-t_{1})-i{\omega_{\rm
s}}(t_{2}^{\prime}-t_{1}^{\prime})}g(t_{1};t)g(t_{2};t)$ (2)
$\displaystyle\mkern-42.0mu\times\big{\langle}\hat{\mathcal{D}}_{{\mathbf{q}_{\rm
i}}{{\bm{\varepsilon}}_{\rm
i}}}^{\dagger}(t_{2})\hat{\mathcal{D}}_{{\mathbf{q}_{\rm
s}}{{\bm{\varepsilon}}_{\rm
s}}}(t_{2}^{\prime})\hat{\mathcal{D}}_{{\mathbf{q}_{\rm
s}}{{\bm{\varepsilon}}_{\rm
s}}}^{\dagger}(t_{1}^{\prime})\hat{\mathcal{D}}_{{\mathbf{q}_{\rm
i}}{{\bm{\varepsilon}}_{\rm i}}}(t_{1})\big{\rangle}$
$\displaystyle\mkern-42.0mu\times
l(t_{1}^{\prime}-t_{1})l(t_{2}^{\prime}-t_{2})dt_{1}dt_{2}dt_{2}^{\prime}dt_{1}^{\prime}\,,$

where ${\bf q}={\mathbf{q}_{\rm i}}-{\mathbf{q}_{\rm s}}$ is the momentum
transfer between incident and scattered photons with energies ${\omega_{\rm
in}}$ and ${\omega_{\rm s}}$ respectively, and
$l(\tau)\\!=\\!e^{-\tau/\tau_{\rm core}}\theta(\tau)$ the core-hole lifetime
decay, with $\tau_{\rm core}\\!=\\!0.5\,t_{h}^{-1}$ as its characteristic
timescale. For a direct transition (e.g. Cu $L$-edge), the dipole operator is
explicitly given by
$\mathcal{D}_{\mathbf{q}{\bm{\varepsilon}}}=\sum_{\textbf{i}\alpha\sigma}e^{-i\mathbf{q}\cdot\textbf{r}_{\textbf{i}}}(A_{\alpha{\bm{\varepsilon}}}d_{\textbf{i}\sigma}^{\dagger}p_{\textbf{i}\alpha\sigma}+h.c.)$
where $A_{\alpha{\bm{\varepsilon}}}$ is the matrix element of the dipole
transition associated with a ${\bm{\varepsilon}}$-polarized photon. As shown
in Fig. 1a, we employed the $\pi-\sigma$ polarization configuration
(${{\bm{\varepsilon}}_{\rm i}}$ parallel and ${{\bm{\varepsilon}}_{\rm s}}$
perpendicular to the scattering plane) and fixed $\alpha_{\rm in}+\alpha_{\rm
out}=50^{\circ}$. This probe condition maximizes the spin-flip response in
presence of the spin-orbit coupling term Ament et al. (2009); Jia et al.
(2014). Calculations for the $\pi$-$\pi$ polarization configuration are
reported in Supplementary note 2 for completeness. Due to the $O(N_{t}^{4})$
complexity of the trRIXS calculation, we adopt the 12D Betts cluster as a
compromise between complexity and finite-size.

Figure 3: Time-dependent renormalization of the spin fluctuations. a,b Time
evolution of the cross-polarized time-resolved resonant inelastic x-ray
scattering (trRIXS) spectra for pump energy $\Omega=10t_{h}$ and amplitude
$A_{0}=0.6$ at a $x=0$ and b $x=0.167$ doping. The momentum transfer is fixed
at $\mathbf{q}=(\pi/2,\pi/2)$. The thin dashed line tracks the instantaneous
peak positions of the spin fluctuation peak. The color bar indicates the
spectral intensity in arbitrary units (a.u.). trRIXS spectra at $t=0$ (center
of the pump pulse) for variable pump amplitudes at c $x=0$ and d $x=0.167$
doping. Spin excitation energy $\omega_{{\bf q}}$ (the peak of the trRIXS
intensity) as function of the pump amplitude at e $x=0$ and f $x=0.167$
doping. The energy width of these spectra includes a broadening due to the
finite x-ray probe duration. The dashed line denotes pure photodoping effects
($\gamma\omega^{\mathrm{eq}}_{\bf q}$, where $\omega^{\mathrm{eq}}_{\bf q}$ is
the numerically evaluated magnon energy at equilibrium and $\gamma$ is the
instantaneous local moment defined in the text), while the solid line
represents a Floquet linear spin wave prediction of the magnon energies
(including photodoping effects) ($\gamma\omega^{\mathrm{FLSW}}_{\bf q}$).

In our calculations, the vector potential of the optical pump field is

$\mathbf{A}^{\rm(pump)}(t)=A_{0}\,e^{-t^{2}/2\sigma_{t}^{2}}\cos(\Omega
t)\,\hat{\mathbf{e}}_{\rm pol},$ (3)

with frequency $\Omega=10\,t_{h}\sim 3.0$ eV (above the Mott gap
$\sim\\!4-5\,t_{h}$) and duration $6\sigma_{t}=18\,t_{h}^{-1}=39.49$ fs, while
the probe has a Gaussian envelope $g(\tau;t)=e^{-(\tau-t)^{2}/2\sigma_{\rm
pr}^{2}}/{\sqrt{2\pi}\sigma_{\rm pr}}$ with $\sigma_{\rm
pr}\\!=\\!1.5\,t_{h}^{-1}$. $A_{0}$ is measured in natural units
$ea_{0}E_{0}/\Omega$, where $A_{0}=0.1$ corresponds to a $E_{0}=79$
mV/$\mathrm{\AA}$ peak electric field. The pump polarization
$\hat{\mathbf{e}}_{\rm pol}$ is linear along the $x$ direction.

Here, we numerically calculate the trRIXS cross section at each time step $t$
by scanning the full four-dimensional hyperplane defined by the time-
integration variables in eq. (2) at zero temperature. We first determine the
equilibrium ground-state wavefunction with the parallel Arnoldi method Lehoucq
et al. (1998); Jia et al. (2018) and, then, calculate the time evolution of
the system through the Krylov subspace technique Manmana et al. (2007); Balzer
et al. (2011).

###  II.2 trRIXS and transient spin fluctuations

In Figure 2, we show selected snapshots of the cross-polarized trRIXS spectrum
at $\mathbf{q}=(\pi/2,\pi/2)$, where the magnon dispersion is maximized [see
Supplementary note 3 for a comparison with $\mathbf{q}=(0,2\pi/3)$). At
resonance (${\omega_{\rm in}}\sim 931$ eV for Cu $L_{3}$-edge in Fig. 2], the
loss spectra below the Mott gap are dominated by a single intense peak at
$\omega\sim 2J=0.3$ eV, where $J=4t_{h}^{2}/U$ is the exchange interaction.
The energy and shape of this peak reflect the distribution of spin
excitations: for an undoped Hubbard model, these manifest as a coherent magnon
due to the strong antiferromagnetic (AFM) order; in contrast, in the $x=0.167$
case, the spin excitations damp into a short-ranged paramagnon Le Tacon et al.
(2011); Dean et al. (2013a, b); Lee et al. (2014); Ishii et al. (2014). At
negative time delays, when the system is close to equilibrium, our
calculations reveal a splitting of the resonance along the ${\omega_{\rm in}}$
axis at $x=0.167$ doping. This reflects the presence of two different chemical
environments, as the valence states can be either half-occupied or empty.

At the center of the pump pulse, the resonance undergoes (i) an overall
suppression, and (ii) a redshift along the energy loss axis for both
compositions. Since the center of the absorption peak does not exhibit a pump-
induced shift along ${\omega_{\rm in}}$, one can fix ${\omega_{\rm in}}=931.5$
eV and examine the time evolution of the magnon/paramagnon peak at ${\bf
q}=(\pi/2,\pi/2)$ [see Figs. 3a and b]. The intensity drop and the peak shift
closely follow the pump envelope, and then persist long after the pump arrival
due to the lack of dissipation. The photoexcited RIXS spectrum is visibly
broadened at $x=0$ and to a lesser extent at $x=0.167$. Next, we analyze the
dependence of the spin fluctuation peak as function of the pump field strength
$A_{0}$ up to a maximum value of 1. This physical regime is relevant to recent
experiments on Sr2IrO4 Dean et al. (2016) and Sr3Ir2O7 Mazzone et al. (2021),
which employed fields strengths $A_{0}$ up to 0.35 and 0.52, respectively. As
$A_{0}$ increases, the peak position for both $x=0$ and $x=0.167$ shifts
monotonically to lower energies [see Figs. 3c and d] by $~{}80\%$ and
$~{}50\%$, respectively. At the same time, the magnon/paramagnon peak broadens
due to the generation of high-order fluctuations beyond simple spin-wave
excitations and of photo-doped carriers which perturb the spin background.
These observations imply that optical laser pulses with negligible momentum
have tangible effects on the spin fluctuations at large momenta.

###  II.3 Magnon/paramagnon energy renormalization

While we report here an exact calculation of the full trRIXS cross-section, it
is useful to discuss the microscopic physics at play in more intuitive terms.
The magnon/paramagnon energy softening in our single-band Hubbard model can be
attributed to photodoping, a dynamical Floquet-type renormalization of the
spin exchange interaction $J$ Mentink and Eckstein (2014); Mentink et al.
(2015), or a combination of the two. Other experimentally relevant mechanisms,
such as magnetophononics Fechner et al. (2018); Nova et al. (2016); Afanasiev
et al. (2021), would require extending the Hubbard Hamiltonian through the
inclusion of additional degrees of freedom.

We first consider the photodoping contribution. Due to the light-induced
holon-doublon excitations, the local spin moment $\langle m_{z}^{2}\rangle$ is
diluted. In the simplest linear spin wave picture, this will induce a
homogeneous softening of the spin structure factors across the entire
Brillouin zone. To quantify this effect, we calculate the instantaneous local
moment at $t=0$, defined as
$\gamma=\sqrt{\frac{\langle\mathbf{S}^{2}(t=0)\rangle}{\langle\mathbf{S}^{2}(t=-\infty)\rangle}}\,$,
and plot the renormalized magnon/paramagnon peak energy
$\gamma\omega_{\mathbf{q}}^{\rm eq}$ as dashed lines in Fig. 3e-f. Here,
$\omega_{\mathbf{q}}^{\rm eq}$ is the numerically-evaluated magnon energy at
equilibrium. As visible from the comparison, the local photodoping effect,
i.e. the dilution of the local spin moment, is too small to account for the
observed energy softening.

The other possible mechanism to explain the peak shift hinges on a light-
induced renormalization of the effective Hamiltonian. It has been demonstrated
that, for an ideal system driven by an infinitely long, nonresonant pump, the
effective spin exchange energy $J$ is modified by a Floquet dressing of the
intermediate doubly-occupied states Mentink et al. (2015). In order to
quantitatively understand our trRIXS spectrum, we apply Floquet linear spin
wave (FLSW) theory at different pump field amplitudes. In this framework, we
assume the system at the center of the pump pulse to be in a steady-state,
thus leading to a closed form for the spin excitations dispersion
$\omega_{\mathbf{q}}^{\rm FLSW}$ (see Methods section.)

By including both the Floquet renormalization and photodoping effects, we
calculate the transient magnon/paramagnon energy
$\gamma\omega_{\mathbf{q}=(\pi/2,\pi/2)}^{\rm FLSW}$ and multiply it by a
constant factor to match the equilibrium peak positions (see Fig. 3e-f). The
scale factors are 1.15 for $x=0$ and 1.02 for $x=0.167$. At $x=0$, the FLSW
prediction tracks the magnon softening for weak pump fields but grossly
deviates at higher pump strength. This is likely due to mobile photo-carriers
which frustrate the AFM ground state Tsutsui et al. (2021) and lower the
energy cost of spin flip events. At $x=0.167$ hole-doping, however, FLSW
theory closely aligns with the trRIXS calculation for all pump strengths. We
tentatively attribute the difference between these two trends to intrinsic
three-site correlated hopping. In the low-energy projection of the Hubbard
Hamiltonian, holes can hop to next- and next-next-nearest neighbors (and thus
lower the total energy) provided that the middle site has an opposite spin
configuration Spałek (1988); Stephan and Horsch (1992); Bal et al. (1995),
which is otherwise inhibited. At equilibrium, this term partially cancels the
doping-induced energy softening Jia et al. (2014); Pärschke et al. (2019) and
is essential to reproduce the spin excitation spectrum of copper oxides.
Active at $x=0.167$, these correlations are silent in the AFM ground state and
build up too slowly to compensate the spin frustration. While both Floquet and
linear spin wave theories over-simplify the nonequilibrium many-body physics
at play, their agreement with the paramagnon energy shift provides evidence of
an ultrafast modification of effective model parameters beyond pure
photodoping.

This striking light-induced softening represents a genuine emergent
phenomenon, which is quantitatively distinct from the behavior of
magnons/paramagnons at equilibrium. Short-range spin excitations in La2CuO4
show negligible dispersion changes for temperature changes well above their
Néel temperature Matsuura et al. (2017), while hole-doped cuprates exhibit at
most a redshift of order $20-30\%$ for a doping change $\Delta x=0.4$ (undoped
to overdoped regime) Le Tacon et al. (2011); Dean et al. (2013a, b); Lee et
al. (2014); Ishii et al. (2014). For a direct comparison, we compute the
transient hole concentrations $\langle n_{h}(t=0)\rangle$ [see Figs. 4a and
Supplementary Note 4], defined as
$n_{h}=\sum_{\textbf{i}}(1-n_{\textbf{i}\uparrow})(1-n_{\textbf{i}\downarrow})/N$.
Due to the presence of quantum fluctuations, the $\langle n_{h}\rangle$ at
equilibrium is 0.06 higher than the nominal hole concentration (0.167). Out of
equilibrium at $t=0$, the strongest pump gives only an extra 0.07 holes per
site. In both the doped and undoped systems, the number of photoinduced holes
is less than the full doping excursion between underdoped and overdoped
regimes, and yet produces a larger softening than what is observed in cuprates
at equilibrium.

Figure 4:  Pump-induced modification of spin excitation energy and pairing
correlations. a Equilibrium (gray) and maximal transient (colored) change of
hole concentration and spin excitation energy (extracted from time-resolved
resonant inelastic x-ray scattering spectral peaks) for the pump conditions in
panel Fig. 3 at half-filling and 16.7% hole doping, respectively. b
Corresponding dynamics of the $d$-wave pairing correlations.

###  II.4 $d$-wave pairing correlations

Since spin fluctuations have been discussed as a possible source of pairing in
the cuprates Scalapino et al. (1986); Gros et al. (1987); Kotliar and Liu
(1988); Schrieffer et al. (1989); Scalapino (1995); Tsuei and Kirtley (2000);
Maier et al. (2006); Scalapino (2012), we examine whether the paramagnon
energy renormalization affects the superconducting pairing. We calculate the
$d$-wave pairing correlation function
$P_{d}=\langle\Delta_{d}^{\dagger}\Delta_{d}\rangle$ at $x=0.167$ for various
pump strengths [see Fig. 4b], with the factor $\Delta_{d}$ defined as
$\Delta_{d}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}d_{\mathbf{k}\uparrow}d_{-\mathbf{k}\downarrow}\left[\cos
k_{x}-\cos k_{y}\right]$. At the pump arrival, $P_{d}$ undergoes a rapid
suppression, thus indicating a decrease of pairing in the $d$-wave channel.
The drop in $P_{d}$ is found to be monotonic with the increase of the pump
strength $A_{0}$. Intuitively, the value of $P_{d}$ is determined by both the
quasiparticle density and the pairing interaction strength. In our
calculations, the pump simultaneously increases quasiparticles via photodoping
and reduces the spin excitation energy. While the former effect is favorable
for superconductivity, it is the latter that dominates and reduces the
nonequilibrium pairing correlations.

##  III Conclusions

We studied the light-induced dynamics of finite-momentum spin excitations in
the prototypical 2D Hubbard model through a full calculation of the
$\pi$-$\sigma$ polarized trRIXS spectrum. By exciting a Mott insulator with an
ultrashort pump pulse, we find that magnons and paramagnons undergo a dramatic
light-induced softening, which cannot be exclusively explained by photodoping.
Paramagnon excitations are quantitatively described by a Floquet
renormalization of the effective exchange interactions, while magnons at high
field strengths deviate from this picture. At optimal doping, we also observe
a sizeable reduction of $d$-wave pairing correlations. While here we explore a
paradigmatic optical excitation at the Mott gap, alternative mechanisms such
as phonon Floquet Hübener et al. (2018), ligand manipulation Ron et al.
(2020), multi-band dynamical effects Tancogne-Dejean et al. (2018); Golež et
al. (2019). and coupling to cavities Gao et al. (2020); Sentef et al. (2020)
might lead to harder spin fluctuations and increased pairing Mentink et al.
(2015); Coulthard et al. (2017). Future polarization-dependent trRIXS
experiments and theoretical calculations will be key to investigate these
photoinduced spin dynamics in a wide variety of quantum materials Cao et al.
(2019); Mitrano and Wang (2020).

##  IV Methods

###  IV.1 Floquet Linear Spin Wave Theory

To describe the evolution of magnetic excitations at finite momentum, we
generalize the local Floquet renormalization of the spin-exchange interaction
Mentink and Eckstein (2014); Mentink et al. (2015) through Floquet linear spin
wave theory. This approach is based on the assumption that the driven system
can be instantaneously approximated by an effective Hamiltonian
$\mathcal{H}_{\rm F}$ and a steady state wavefunction $|\psi_{F}\rangle$. In
the Floquet framework, the hopping $t_{ij}$ in the effective Hamiltonian
$\mathcal{H}_{\rm F}$ is renormalized by a factor
$\mathcal{J}_{m}(\mathbf{A}\cdot\mathbf{r}_{ij})$, where $\mathcal{J}_{m}(x)$
is the Bessel function of the first kind and $\mathbf{A}$ is the vector
potential of an infinitely long driving field. Following second-order
perturbation theory, the effective spin-exchange interaction between two
coordinates is renormalized into Mentink et al. (2015)

$\displaystyle J_{ij}^{\rm
eff}=\frac{4t_{ij}^{2}}{U}\sum_{m=-\infty}^{\infty}\frac{\mathcal{J}_{m}(\mathbf{A}\cdot\mathbf{r}_{ij})^{2}}{1+m\Omega/U}\,.$
(4)

Then, through the standard Holstein–Primakoff transformation and large-$S$
expansion, one can obtain the linearized bosonic Hamiltonian

$\displaystyle\mathcal{H}_{FLSW}=\sum_{\mathbf{q}}\sqrt{A_{\mathbf{q}}^{2}-B_{\mathbf{q}}^{2}}\,a_{\mathbf{q}}^{\dagger}a_{\mathbf{q}}=\sum_{\mathbf{q}}\omega_{\mathbf{q}}^{\rm
FLSW}a_{\mathbf{q}}^{\dagger}a_{\mathbf{q}}\,,$ (5)

where $a_{\mathbf{q}}$ is the magnon annihilation operator and the
coefficients $A_{\mathbf{q}}$ and $B_{\mathbf{q}}$ in the two-dimensional
plane are Delannoy et al. (2009)

$\displaystyle A_{\mathbf{q}}$ $\displaystyle=$ $\displaystyle J_{x}^{\rm
eff}+J_{y}^{\rm eff}-J_{xy}^{\rm eff}\left[1-\cos(q_{x}+q_{y})\right]$ (6)
$\displaystyle-J_{yx}^{\rm eff}\left[1-\cos(q_{x}-q_{y})\right]\,,$
$\displaystyle B_{\mathbf{q}}$ $\displaystyle=$ $\displaystyle J_{x}^{\rm
eff}\cos q_{x}+J_{y}^{\rm eff}\cos q_{y}\,.$ (7)

Here, unlike the equilibrium counterpart, the coefficients are anisotropic due
to the electric field polarization:

$\displaystyle J_{x/y}^{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{4t_{h}^{2}}{U}\sum_{m=-\infty}^{\infty}\frac{\mathcal{J}_{m}(\mathbf{A}\cdot\mathbf{e}_{x/y})^{2}}{1+m\Omega/U}\,,$
(8) $\displaystyle J_{xy/yx}^{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{4t_{h}^{\prime
2}}{U}\sum_{m=-\infty}^{\infty}\frac{\mathcal{J}_{m}(\mathbf{A}\cdot(\mathbf{e}_{x}\pm\mathbf{e}_{y}))^{2}}{1+m\Omega/U}\,.$
(9)

This derivation can be extended to higher order following the strategy of
equilibrium LSW theory Delannoy et al. (2009).

## Data Availability

The numerical data that support the findings of this study are available from
the corresponding authors upon reasonable request.

## Code Availability

The relevant scripts of this study are available from the corresponding
authors upon reasonable request.

## Acknowledgements

We thank D.R. Baykusheva, M. Buzzi, C.-C. Chen, M. P. M. Dean, A. A. Husain,
C. J. Jia, D. Nicoletti, and M. Sentef for insightful discussions. Y.C.,
T.P.D., and B.M. were supported by the US Department of Energy, Office of
Basic Energy Sciences, Division of Materials Sciences and Engineering, under
Contract No. DE-AC02-76SF00515. M. M. was supported by the William F. Milton
Fund at Harvard University. This research used resources of the National
Energy Research Scientific Computing Center (NERSC), a U.S. Department of
Energy Office of Science User Facility operated under Contract No. DE-
AC02-05CH11231.

## Author contributions

Y.W. and M.M. conceived the project. Y.W. and Y.C. performed the calculations.
Y.W. and M.M. performed the data analysis. Y.W. and M.M. wrote the paper with
the help from Y.C., B.M., and T.P.D.

## Competing interests

The authors declare no competing interests.

## References

## References

  * Beaurepaire et al. (1996) E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Phys. Rev. Lett. 76, 4250 (1996).
  * Kirilyuk et al. (2010) A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * Kampfrath et al. (2011) T. Kampfrath, A. Sell, G. Klatt, A. Pashkin, S. Mährlein, T. Dekorsy, M. Wolf, M. Fiebig, A. Leitenstorfer, and R. Huber, Nat. Photonics 5, 31 (2011).
  * Nova et al. (2016) T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, Nat. Phys. 13, 132 (2016).
  * Disa et al. (2020) A. S. Disa, M. Fechner, T. F. Nova, B. Liu, M. Först, D. Prabhakaran, P. G. Radaelli, and A. Cavalleri, Nat. Phys. 16, 937 (2020).
  * Mikhaylovskiy et al. (2015) R. V. Mikhaylovskiy, E. Hendry, A. Secchi, J. H. Mentink, M. Eckstein, A. Wu, R. V. Pisarev, V. V. Kruglyak, M. I. Katsnelson, T. Rasing, et al., Nat. Commun. 6, 8190 (2015).
  * Ron et al. (2020) A. Ron, S. Chaudhary, G. Zhang, H. Ning, E. Zoghlin, S. D. Wilson, R. D. Averitt, G. Refael, and D. Hsieh, Phys. Rev. Lett. 125, 197203 (2020).
  * Rossat-Mignod et al. (1991) J. Rossat-Mignod, L. P. Regnault, C. Vettier, P. Bourges, P. Burlet, J. Bossy, J. Henry, and G. Lapertot, Physica C: Superconductivity 185, 86 (1991).
  * Mook et al. (1993) H. A. Mook, M. Yethiraj, G. Aeppli, T. E. Mason, and T. Armstrong, Phys. Rev. Lett. 70, 3490 (1993).
  * Fong et al. (1999) H. Fong, P. Bourges, Y. Sidis, L. P. Regnault, A. Ivanov, G. D. Gu, N. Koshizuka, and B. Keimer, Nature 398, 588 (1999).
  * Scalapino et al. (1986) D. Scalapino, E. Loh Jr, and J. Hirsch, Phys. Rev. B 34, 8190 (1986).
  * Gros et al. (1987) C. Gros, R. Joynt, and T. Rice, Zeitschrift für Physik B Condensed Matter 68, 425 (1987).
  * Kotliar and Liu (1988) G. Kotliar and J. Liu, Phys. Rev. B 38, 5142 (1988).
  * Schrieffer et al. (1989) J. Schrieffer, X. Wen, and S. Zhang, Phys. Rev. B 39, 11663 (1989).
  * Scalapino (1995) D. J. Scalapino, Phys. Rep. 250, 329 (1995).
  * Tsuei and Kirtley (2000) C. Tsuei and J. Kirtley, Rev. Mod. Phys. 72, 969 (2000).
  * Maier et al. (2006) T. A. Maier, M. S. Jarrell, and D. J. Scalapino, Phys. Rev. Lett. 96, 047005 (2006).
  * Fausti et al. (2011) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * Hu et al. (2014) W. Hu, S. Kaiser, D. Nicoletti, C. R. Hunt, I. Gierz, M. C. Hoffmann, M. Le Tacon, T. Loew, B. Keimer, and A. Cavalleri, Nat. Mater. 13, 705–711 (2014).
  * Kaiser et al. (2014) S. Kaiser, C. R. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Y. Liu, M. Le Tacon, T. Loew, D. Haug, and B. Keimer, Phys. Rev. B 89, 184516 (2014).
  * Nicoletti et al. (2014) D. Nicoletti, E. Casandruc, Y. Laplace, V. Khanna, C. R. Hunt, S. Kaiser, S. S. Dhesi, G. D. Gu, J. P. Hill, and A. Cavalleri, Phys. Rev. B 90, 100503 (2014).
  * Zhao et al. (2011) H. Zhao, D. Talbayev, X. Ma, Y. Ren, A. Venimadhav, Q. Li, and G. Lüpke, Phys. Rev. Lett. 107, 207205 (2011).
  * Dal Conte et al. (2012) S. Dal Conte, C. Giannetti, G. Coslovich, F. Cilento, D. Bossini, T. Abebaw, F. Banfi, G. Ferrini, H. Eisaki, M. Greven, et al., Science 335, 1600 (2012).
  * Batignani et al. (2015) G. Batignani, D. Bossini, N. Di Palo, C. Ferrante, E. Pontecorvo, G. Cerullo, A. Kimel, and T. Scopigno, Nat. Photonics 9, 506 (2015).
  * Yang et al. (2020) J. A. Yang, N. Pellatz, T. Wolf, R. Nandkishore, and D. Reznik, Nat Commun. 11, 2548 (2020).
  * Ament et al. (2009) L. J. Ament, G. Ghiringhelli, M. M. Sala, L. Braicovich, and J. van den Brink, Phys. Rev. Lett. 103, 117003 (2009).
  * Haverkort (2010) M. W. Haverkort, Phys. Rev. Lett. 105, 167404 (2010).
  * Ament et al. (2011) L. J. P. Ament, M. van Veenendaal, T. P. Devereaux, J. P. Hill, and J. van den Brink, Rev. Mod. Phys. 83, 705 (2011).
  * Dean et al. (2016) M. Dean, Y. Cao, X. Liu, S. Wall, D. Zhu, R. Mankowsky, V. Thampy, X. Chen, J. Vale, D. Casa, et al., Nat. Mater. 15, 601 (2016).
  * Mazzone et al. (2021) D. G. Mazzone, D. Meyers, Y. Cao, J. G. Vale, C. D. Dashwood, Y. Shi, A. J. A. James, N. J. Robinson, J. Lin, V. Thampy, et al., Proceedings of the National Academy of Sciences 118 (2021).
  * Jia et al. (2014) C. J. Jia, E. A. Nowadnick, K. Wohlfeld, Y. F. Kung, C. C. Chen, S. Johnston, T. Tohyama, B. Moritz, and T. P. Devereaux, Nat. Commun. 5, 3314 (2014).
  * Jia et al. (2016) C. Jia, K. Wohlfeld, Y. Wang, B. Moritz, and T. P. Devereaux, Physical Review X 6, 021020 (2016).
  * Tsutsui et al. (2000) K. Tsutsui, H. Kondo, T. Tohyama, and S. Maekawa, Physica B: Condensed Matter 284, 457 (2000).
  * Kourtis et al. (2012) S. Kourtis, J. van den Brink, and M. Daghofer, Phys. Rev. B 85, 064423 (2012).
  * Chen et al. (2019) Y. Chen, Y. Wang, C. Jia, B. Moritz, A. M. Shvaika, J. K. Freericks, and T. P. Devereaux, Phys. Rev. B 99, 104306 (2019).
  * Wang et al. (2020) Y. Wang, Y. Chen, C. Jia, B. Moritz, and T. P. Devereaux, Phys. Rev. B 101, 165126 (2020).
  * Lehoucq et al. (1998) R. B. Lehoucq, D. C. Sorensen, and C. Yang, _ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods_ (Siam, 1998).
  * Jia et al. (2018) C. Jia, Y. Wang, C. Mendl, B. Moritz, and T. Devereaux, Comput. Phys. Commun. 224, 81 (2018).
  * Manmana et al. (2007) S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, Phys. Rev. Lett. 98, 210405 (2007).
  * Balzer et al. (2011) M. Balzer, N. Gdaniec, and M. Potthoff, J. Phys. Condens. Matter 24, 035603 (2011).
  * Le Tacon et al. (2011) M. Le Tacon, G. Ghiringhelli, J. Chaloupka, M. M. Sala, V. Hinkov, M. Haverkort, M. Minola, M. Bakr, K. Zhou, S. Blanco-Canosa, et al., Nat. Phys. 7, 725 (2011).
  * Dean et al. (2013a) M. Dean, G. Dellea, R. Springell, F. Yakhou-Harris, K. Kummer, N. Brookes, X. Liu, Y. Sun, J. Strle, T. Schmitt, et al., Nat. Mater. 12, 1019 (2013a).
  * Dean et al. (2013b) M. Dean, A. James, R. Springell, X. Liu, C. Monney, K. Zhou, R. Konik, J. Wen, Z. Xu, G. Gu, et al., Phys. Rev. Lett. 110, 147001 (2013b).
  * Lee et al. (2014) W. Lee, J. Lee, E. Nowadnick, S. Gerber, W. Tabis, S. Huang, V. Strocov, E. Motoyama, G. Yu, B. Moritz, et al., Nat. Phys. 10, 883 (2014).
  * Ishii et al. (2014) K. Ishii, M. Fujita, T. Sasaki, M. Minola, G. Dellea, C. Mazzoli, K. Kummer, G. Ghiringhelli, L. Braicovich, T. Tohyama, et al., Nat. Commun. 5, 3714 (2014).
  * Mentink and Eckstein (2014) J. H. Mentink and M. Eckstein, Phys. Rev. Lett. 113, 057201 (2014).
  * Mentink et al. (2015) J. H. Mentink, K. Balzer, and M. Eckstein, Nat. Commun. 6, 6708 (2015).
  * Fechner et al. (2018) M. Fechner, A. Sukhov, L. Chotorlishvili, C. Kenel, J. Berakdar, and N. A. Spaldin, Phys. Rev. Materials 2, 064401 (2018).
  * Afanasiev et al. (2021) D. Afanasiev, J. Hortensius, B. Ivanov, A. Sasani, E. Bousquet, Y. M. Blanter, R. V. Mikhaylovskiy, A. V. Kimel, and A. D. Caviglia, Nat. Mater. 20, 607 (2021).
  * Tsutsui et al. (2021) K. Tsutsui, K. Shinjo, and T. Tohyama, Phys. Rev. Lett. 126, 127404 (2021).
  * Spałek (1988) J. Spałek, Phys. Rev. B 37, 533 (1988).
  * Stephan and Horsch (1992) W. Stephan and P. Horsch, Int. J Mod. Phys. B 6, 589 (1992).
  * Bal et al. (1995) J. Bal, A. M. Oleś, and J. Zaanen, Phys. Rev. B 52, 4597 (1995).
  * Pärschke et al. (2019) E. M. Pärschke, Y. Wang, B. Moritz, T. P. Devereaux, C.-C. Chen, and K. Wohlfeld, Phys. Rev. B 99, 205102 (2019).
  * Matsuura et al. (2017) M. Matsuura, S. Kawamura, M. Fujita, R. Kajimoto, and K. Yamada, Phys. Rev. B 95, 024504 (2017).
  * Scalapino (2012) D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Hübener et al. (2018) H. Hübener, U. De Giovannini, and A. Rubio, Nano Lett. 18, 1535 (2018).
  * Tancogne-Dejean et al. (2018) N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, Phys. Rev. Lett. 121, 097402 (2018).
  * Golež et al. (2019) D. Golež, M. Eckstein, and P. Werner, Phys. Rev. B 100, 235117 (2019).
  * Gao et al. (2020) H. Gao, F. Schlawin, M. Buzzi, A. Cavalleri, and D. Jaksch, Phys. Rev. Lett. 125, 053602 (2020).
  * Sentef et al. (2020) M. A. Sentef, J. Li, F. Künzel, and M. Eckstein, Phys. Rev. Research 2, 033033 (2020).
  * Coulthard et al. (2017) J. Coulthard, S. R. Clark, S. Al-Assam, A. Cavalleri, and D. Jaksch, Phys. Rev. B 96, 085104 (2017).
  * Cao et al. (2019) Y. Cao, D. Mazzone, D. Meyers, J. Hill, X. Liu, S. Wall, and M. Dean, Philosophical Transactions of the Royal Society A 377, 20170480 (2019).
  * Mitrano and Wang (2020) M. Mitrano and Y. Wang, Commun. Phys. 3, 184 (2020).
  * Delannoy et al. (2009) J.-Y. Delannoy, M. Gingras, P. Holdsworth, and A.-M. Tremblay, Phys. Rev. B 79, 235130 (2009).


2021

[1]Jagdish Chand Bansal

1] South Asian University, India

2]Liverpool Hope University, United Kingdom

# Drone Flocking Optimization using NSGA-II and Prinicipal Component Analysis

<EMAIL_ADDRESS>Nikhil Sethi<EMAIL_ADDRESS>Ogbonnaya Anicho
<EMAIL_ADDRESS>Atulya Nagar<EMAIL_ADDRESS>[  [

###### Abstract

Individual agents in natural systems like flocks of birds or schools of fish
display a remarkable ability to coordinate and communicate in local groups and
execute a variety of tasks efficiently. Emulating such natural systems into
drone swarms to solve problems in Defence, agriculture, industry automation
and humanitarian relief is an emerging technology. However, flocking of aerial
robots while maintaining multiple objectives, like collision avoidance, high
speed etc. is still a challenge. In this paper, optimized flocking of drones
in a confined environment with multiple conflicting objectives is proposed.
The considered objectives are collision avoidance (with each other and the
wall), speed, correlation, and communication (connected and disconnected
agents). Principal Component Analysis (PCA) is applied for dimensionality
reduction, and understanding the collective dynamics of the swarm. The control
model is characterised by 12 parameters which are then optimized using a
multi-objective solver (NSGA-II). The obtained results are reported and
compared with that of the CMA-ES algorithm. The study is particularly useful
as the proposed optimizer outputs a Pareto Front representing different types
of swarms which can applied to different scenarios in the real world.

###### keywords:

Drone swarm, Multi-Objective Optimization, PCA, NSGA-II, Drone swarm
simulator, Collective dynamics

##  1 Introduction

Collective behaviour is pervasive in nature and is frequently observed in
diverse organisms ranging from microscopic bacteria allison_bacterial_1991  to
large scale flocking of birds and insectsnagy_hierarchical_2010 annurev-
ento-010814-020627 . While researchers still hypothesize the underlying
mechanisms behind such behaviour, moving in groups can offer several
advantages like avoiding predators or carrying collective cargoron2018bi .
Emulating these natural systems has gained popularity in the past few years
and the development of a robust, fault-tolerant and generalised swarm of
robots is now a widely regarded problem among researchers saffre2021design
coppola2020survey . Aerial swarms, owing to their high maneuverability and
speed, find a number of applications in various industries. They can be
deployed as counterdrone measures brust_defending_2018  or basic surveillance
operations in a defence scenario. In abraham2019swarm , the authors utiltise
the sensing capability of multiple robots to yield topographical and
population density maps of a disaster afflicted area. In
tosato_autonomous_2019 , a centralised swarm architecture was proposed for
measuring air pollution. A system like this could potentially reduce
measurement error due to the bigger sample size and distributed data points
over the coverage volume. Mixed aerial and ground swarms have also been used
for automation in construction taskskrizmancic_cooperative_2020 . In
ju2018multiple , multiple UAVs have been shown to outperform a single UAV for
tasks like agricultural sensing and monitoring by measuring multiple metrics
like energy consumption, flight time, and area coverage. In general, drone
swarms can be classified into three categories in order of increasing
complexity:

  * •

Coordinated: This refers to the collective movement with basic environmental
awareness and collision avoidance.

  * •

Cooperative: Here, the robots start to work together to achieve a certain goal
using lesser resources than that of a single drone.

  * •

Collaborative: This refers to multiple drones working together irrespective of
their nature, i.e. heterogeneous collaboration.

In this paper, we propose a methodology to solve the drone swarm coordination
problem with multiple conflicting objectives. In this document, the terms
drone, UAV and aerial robot are used interchangeably.

Developing a robust velocity controller that allows multiple drones to self-
organise comes with it’s own set of challenges. According to the taxonomy
defined in coppola2020survey  the control of velocities comes under Swarming
Behaviour i.e. deciding a high level control policy with shared information
across each agent’s neighbors. The challenge is to take this shared
information (the agent’s own state as well states of neighbors) and come up
with functions (policies) that output an instantaneous velocity vector for
each agent. Over time, each agent’s velocity gives rise to various patterns
and mutual interactions that can potentially emerge into self-organising
behaviour. Conventionally, a first of it’s kind algorithm by Reynolds was
based on simple rules for each agent and has been successfully applied in many
fields hauert2011reynolds  dewi_simulating_2011  moere_time-varying_2004 . In
vasarhelyi2018optimized , the authors address this problem by defining a
single fitness function and optimizing it through the CMA-ES algorithm.
However, the study doesn’t take into account multiple conflicting objectives,
the priority of which can vary depending on the scenario. In marquez2021multi
, a multi-objective solution for quad-rotors is proposed but the swarm size is
limited and the full range of solutions that take into account the relations
among the fitnesses is not explored. We explore these relationships using
unsupervised learning and extend our findings to highlight the use of
obtaining a non-dominating solution set for drone flocking.

Simulation is a powerful tool while studying these systems as it allows risk
free experimentation and many programming platforms have been leveraged to
come up with such solutions shah2018airsim soria2020swarmlab 6920950 . In
shah2018airsim , the authors create a 3D simulator
(AirSim111https://github.com/microsoft/AirSim) for autonomous vehicles written
in C++ with growing support for multi-drone scenarios. It is easy to setup a
simulation through Python scripting. The goal of AirSim however, is graphical
realism and it is cumbersome to extend the same for the domain of swarm
intelligence. RobotSim 222https://github.com/csviragh/robotsim, a fast
simulator targeting flocking scenarios, is written in C but lacks sufficient
abstraction and extension to multi-objective optimization which is required
for the work carried out in this paper. In soria2020swarmlab , the authors
develop SwarmLAB333https://github.com/lis-epfl/swarmlab, a dedicated drone
swarm simulator capable of handling non-linear quadcopter simulation. While
the simulator itself could meet requirements of speed and accuracy, it was
insufficient for extending future work in the direction of hardware testing.
We use Python to develop such a simulator as it can leverage existing
libraries and algorithms. The inter-operability with libraries like Pymavlink,
and the ease of testing protocols like the MultiWii Serial Protocol (MSP) on
Raspberry pi on board computers was found to be particularly easy with Python
based on past experience with similar projects.

Modelling of natural processes through simulation often needs to be
complemented by an in depth qualitative understanding of the performance
measures. Unsupervised learning can help in understanding and clustering data
especially in high dimensional spaces which can’t be visualised. It is widely
used in experiments where abundant data is available such as mapping
vulnerability indices abson2012using , understanding relationships between
economical and environmental objectives in a chemical supply chain pozo2012use
, understanding global motions of atoms in proteins loeffler2009collective ,
and most commonly for dimensionality reduction in evolutionary algorithms
deb2006searching . Similar to many real systems, the solution of an
optimization problem depends on various factors. Often, these factors or
objectives are conflicting in nature and they cannot be solved simultaneously
without compromising on the overall fitness. In case of flocking we consider
six objectives:

  * •

Collision avoidance with the wall.

  * •

Collision avoidance with each other.

  * •

Average speed of the swarm.

  * •

Average velocity alignment or correlation.

  * •

Total number of connected agents.

  * •

Total number of disconnected agents.

We use PCA for both understanding collective dynamics of multi-agent systems
and therefore reducing the objective space for the multi-objective optimizer.
To the best of our knowledge, this work is the first attempt that involves
using PCA to reduce the objective functions for a drone flocking optimization
problem.

These objectives are then optimized via a well-established multi-objective
optimizer (NSGA-II) to yield a Pareto front which can be used to guide
decision making and trade-offs under various situations. We report the results
and show that the results at extremities of the Pareto front perform better
than that of the CMA-ES algorithm. We conclude by giving some practical
examples of the use of such abstract mathematical formalism for real time
decision making with a flock of UAVs.

In short, in this research, we create a drone swarm simulator integrated with
a multi-objective solver, use PCA to understand the collective dynamics of
swarms, and give a Pareto front that represents different swarms that can be
used in real-world scenarios. The rest of the paper is organised as follows:
Section 2 presents the background of Principal Component Analysis (PCA) and
multi-objective optimizer (NSGA-II). A Drone flocking optimization problem is
formulated in Section 3. In Section 4 PCA is used to reduce the number of the
objective functions and a discussion on the correlations is followed. Section
5 presents the experimental setup while the numerical results and discussions
are given in Section 6. The research is concluded by giving some potential
use-cases and possible future work.

##  2 Background

###  2.1 Principal Component Analysis

A high dimensional objective space is known to suffer from problems like poor
selection pressure and convergencedeb2006searching . It is also difficult to
visualize the space and gain intuition which is often required for appropriate
decision making. Principal component analysis, a technique under the domain of
unsupervised learning may be useful to understand the underlying structure of
the data without explicit labels. The idea is to search for the eigenvectors
of an $m$-dimensional covariance matrix ($K$) which is then used to decide the
redundant objectives. Here, $m$ is the number of objectives. This covariance
matrix (often called correlation matrix when scaled) is symmetric and it’s
elements give the relations between the design variables on which the analysis
has been run. Such an analysis on the objectives of a optimization problem can
give insights about their correlations and can help in understanding their
qualitative aspects.

Let $X$ be an $n$ x $m$ design matrix with $n$ rows as the samples and $m$
columns as objectives. A pre-processing step often carried out is the
normalisation of design matrix to bring the mean of samples for each objective
to 0 and the variance to 1.0 (Eq. (1)). The covariance matrix is then
calculated by taking the mean of all samples of the pairwise products for each
objective (Eq. (2)). In a vectorized format, this is equivalent to taking the
matrix product of the design matrix $X$ with it’s transpose (Eq. (3)).

$\displaystyle X_{ij}^{norm}$
$\displaystyle=\frac{X_{ij}-\mu_{j}}{\sigma_{j}}$ (1) $\displaystyle K_{ij}$
$\displaystyle=\frac{1}{n}\sum^{n}_{k=1}X_{ki}X_{kj}$ (2) $\displaystyle K$
$\displaystyle=\frac{1}{n}(X^{norm})^{T}X^{norm}$ (3)

where,

$X_{ij}\ \equiv$ Element of $X$ at $i^{th}$ row and $j^{th}$ column

$X_{ij}^{norm}\ \equiv X_{ij}$ normalised to 0.0 mean and 1.00 standard
deviation

$\mu_{j}\ \equiv$ Mean of all $n$ samples of $j^{th}$ objective

$\sigma_{j}\ \equiv$ Standard deviation of all $n$ samples of $j^{th}$
objective

$n\ \equiv$ Number of samples

$m\ \equiv$ Number of objectives

$K\equiv$ Covariance matrix

$K_{ij}\ \equiv$ Element of $K$ at $i^{th}$ row and $j^{th}$ column

###  2.2 Non-Dominating Sorting Genetic Algorithm-II

NSGA-II is a multi-objective optimization algorithm that is based on ranking
each solution in the population according to their fitness and progressively
producing better solutions using genetic operators like reproduction and
mutation. The entire algorithm is explained in detail in deb_fast_2002 .
However, a brief explanation covering the salient features of NSGA-II is
explained here.

Let $P_{o}$ be a N sized initial random population. This population is sorted
based on non-domination according to the following rules: An individual
$X_{1}$ in the population is said to be dominated by individual $X_{2}$ if
satisfies both of the following conditions:

  * •

All fitnesses of $X_{1}$ must be less than or equal to that of $X_{2}$
particle.

  * •

At least one fitness of $X_{1}$ must be strictly less than that of $X_{2}$.

Mathematically, individual $X_{1}$ dominates $X_{2}$ if $d=1$, and the
individuals are non-dominated if $d=0$.

Where, $d=\\{\forall m\ F(X_{1})^{m}\leq F(X_{2})^{m}\\}\cap\\{\exists m\
F(X_{1})^{m}<F(X_{2})^{m}\\}$

This method divides the population into dominating and non-dominating
solutions which is a heuristic used to guide the population towards better
solutions through the generations. Each solution in this population is also
ranked based on the number of other members it is dominated by and accordingly
it is assigned a front rank. Next, an offspring population Q is created from
the sorted population by applying tournament selection, recombination and
mutation operators. To ensure that the best solution across generations
remains (elitism), a new 2N sized population is made using $P\cup Q$ which is
again sorted and ranked. To make the next population $P_{t+1}$ from this
combined set, solutions are taken in order of their front ranking. In case the
number of solutions belonging to a front exceed the amount that can be
accommodated into the new N sized population, the remaining solutions in that
front are ranked based on a crowding operator as follows:

Let $\mathcal{F}^{k}$ be the set of solutions on the $k^{th}$ ranked pareto
front. The crowding distance ($c^{m}_{i}$) for the $m^{th}$ objective for
$i^{th}$ solution on this front is defined as the normalised distance between
the two nearest solutions i.e. $(i+1)^{th}$ and $(i-1)^{th}$ (Eq. (4)). The
overall crowding distance ($c_{i}$) is the sum taken for each objective (Eq.
(5)).

$\displaystyle\forall X_{i}\in\mathcal{F}^{k}:c^{m}_{i}$
$\displaystyle=\frac{F^{m}(X_{i+1})-F^{m}(X_{i-1})}{F^{m}_{max}-F^{m}_{min}}$
(4) $\displaystyle c_{i}$ $\displaystyle=\sum^{M}_{m=1}c^{m}_{i}$ (5)

This crowding operator ensures that the Pareto Front is uniformly distributed
and the range of each each objective value is maximised as the search
progresses. The remaining solutions are ranked according to $c_{i}$ and the
new population $P_{t+1}$ moves forward to the next generation. NSGA-II is
faster than NSGA-I and has a worst case complexity of $O(MN^{2})$.

##  3 Drone flocking optimization problem

A completely decentralized flocking swarm is based on certain simple rules
like Separation, Alignment, Cohesion. These rules when defined using a
velocity control algorithm has certain parameters which can be tuned to flock
optimally. In this section these parameters are introduced and a simulation
framework capable of handling artificial sensor noise is created. The
algorithm used for flocking is based on the work in vasarhelyi2018optimized
and Reynold’s Flocking model 10.1145/37402.37406 . Some subtle modifications
have been incorporated in order to handle a multi-objective optimization
framework. We use vectorized versions of the equations to leverage fast
computation with matrix computation libraries.

To simulate a multi-agent system, there must be a mechanism to share
information across the agents. In case of a decentralised system this
information is shared in each agent’s neighborhood $\mathcal{N}_{o}$.
Moreover, real-systems are characterised by a stochastic uncertainty and noise
which are incorporated into the position(r) and velocity(v) vectors of the
drones. The model for simulating the noise is taken from viragh2014flocking .
The relative position ($\textbf{r}_{ji}$) and velocity ($\textbf{v}_{ji}$) at
time $t$ is then found using the following equations:

$\displaystyle\textbf{r}_{ji}(t)$
$\displaystyle=(\textbf{r}_{j}(t-t_{del})+\textbf{r}_{j}^{gps})-\textbf{r}_{i}(t)-\textbf{r}_{i}^{gps}$
(6) $\displaystyle\textbf{v}_{ji}(t)$
$\displaystyle=(\textbf{v}_{j}(t-t_{del})+\textbf{v}_{j}^{gps})-\textbf{v}_{i}(t)-\textbf{v}_{i}^{gps}$
(7) $\displaystyle R^{rel}_{j}$ $\displaystyle=\textbf{r}_{ji}(t)$ (8)
$\displaystyle V^{rel}_{j}$ $\displaystyle=\textbf{v}_{ji}(t)$ (9)

where,

$\textbf{r}_{ji}\ \equiv$ Relative position vector of $j^{th}$ agent with
respect to $i^{th}$ agent at time $t$

$\textbf{v}_{ji}\ \equiv$ Relative velocity vector of $j^{th}$ agent with
respect to $i^{th}$ agent at time $t$

$R^{rel}_{j}\ \equiv$ $j^{th}$ row of the Relative position matrix for agent
$i\ \forall\ j={1,2,...,\mathcal{N}_{o}}$

$V^{rel}_{j}\ \equiv$ $j^{th}$ row of the Relative velocity matrix for agent
$i\ \forall\ j={1,2,...,\mathcal{N}_{o}}$

$t_{del}\ \equiv$ Simulated communication delay

$\textbf{r}^{gps}\ \equiv$ Simulated GPS noise for position

$\textbf{v}^{gps}\ \equiv$ Simulated GPS noise for velocity

###  3.1 Decision variables

Based on the above modification for the relative position and velocities the
flocking rules are explained the following sections. These rules give rise to
certain parameters which are used as decision variables for the drone flocking
optimization problem. Note that all the flocking operations are carried out
for all $N$ agents .

####  3.1.1 Separation

To flock effectively without collisions, the agents must have a mechanism for
repulsion. A spring-like mechanism is used which is activated at short ranges
of inter-agent distance in the flock. The following two equations (11) and
(12) are used to find a repulsion vector for agent $i$ after scaling it
according to the relative distances in $\textbf{r}_{i}^{mag}$ and a gain
$p^{rep}$.

$\displaystyle\textbf{r}^{mag}$ $\displaystyle=\|R^{rel}\|^{\bot r_{0}^{rep}}$
(10) $\displaystyle V^{rep}$
$\displaystyle=p^{rep}.(\textbf{r}^{mag}-r_{0}^{rep}).\frac{R^{rel}}{\textbf{r}^{mag}}$
(11) $\displaystyle\textbf{v}^{rep}_{i}$
$\displaystyle=\sum_{j=1}^{\mathcal{N}_{o}}V^{rep}_{j}$ (12)

where,

$\ a^{\bot c}_{\top b}=min(max(a,b),c))$

$\textbf{r}^{mag}\equiv\mathcal{N}_{o}\ \text{x}\ 1$ sized vector containing
inter-agent distances

$r_{0}^{rep}\ \equiv$ Repulsion cutoff distance (user dependent parameter)

$p^{rep}\ \equiv$ Repulsion gain (user dependent parameter)

$V^{rep}\ \equiv\mathcal{N}_{o}\ \text{x}\ 2$ sized matrix of scaled repulsion
velocities

$V^{rep}_{j}\ \equiv$ Repulsion velocity of $j^{th}$ neighbor

$\textbf{v}^{rep}_{i}\ \equiv$ Desired collective repulsion vector

Note that the upper bound of $\textbf{r}^{mag}$ is the parameter $r_{0}^{rep}$
to enable short range effects. The matrix norm in Eq. (10) is only taken along
the row axis, i.e. for each neighbor. $V^{rep}$ contains all the the
corresponding scaled repulsion velocities and the division and multiplication
in Eq. (11) is done element wise.

####  3.1.2 Alignment

Vásárhelyi et al. vasarhelyi2018optimized  realised that effective control of
_both_ the magnitude and direction of velocities as a function of inter-agent
distances can yield the best alignment with scalable velocities. The equations
for alignment are similar to repulsion with one major difference: the upper
bound for the velocity magnitude ($\textbf{v}^{frictmax}_{i}$) is now
calculated dynamically with decay function $D$ in Eq. (13) which is dependent
on the inter-agent distance vasarhelyi2018optimized . Eqs. (14) - (16)
describe the process of finding out the combined alignment vector for the
agent.

$\displaystyle\textbf{v}^{frictmax}$
$\displaystyle=D(\textbf{r}_{i}^{mag}-r_{0}^{frict}-r_{0}^{rep},a^{frict},p^{frict})_{\top{v}^{frict}}$
(13) $\displaystyle\textbf{v}^{mag}$
$\displaystyle=\|V^{rel}\|_{\top\textbf{v}^{frictmax}}$ (14) $\displaystyle
V^{frict}$
$\displaystyle=c^{frict}.(\textbf{v}^{mag}-\textbf{v}^{frictmax}).\frac{V^{rel}}{\textbf{v}_{mag}}$
(15) $\displaystyle\textbf{v}^{frict}_{i}$
$\displaystyle=\sum_{j=1}^{\mathcal{N}_{o}}V^{frict}_{j}$ (16)

where,

$D$ is a vectorized version of the velocity decay function taken from
vasarhelyi2018optimized

$p^{frict}\equiv$ Slope for the linear part of the decay curve (user dependent
parameter)

$a^{frict}\equiv$ Acceleration for the non-linear part of the decay curve
(user dependent parameter)

$c^{frict}\equiv$ Overall Gain for alignment (user dependent parameter)

$v^{frict}\equiv$ Velocity slack for alignment (user dependent parameter)

$r_{0}^{frict}\ \equiv$ Alignment cutoff distance for maximum alignment (user
dependent parameter)

$\textbf{v}^{frict}_{i}\ \equiv$ Desired collective alignment vector

Eq. (13) gives a vector composed of the maximum allowable velocity difference
for each neighbor. The maximum is proportional to the inter-agent distance.
This ensures that the alignment for two agents that are in close proximity is
larger and vice-versa. Also, the maximum allowable difference is lower bound
by an optimization parameter $v^{frict}$ so that the agents do not strive for
perfect alignment and there is some slack. Eq. (15) compensates the velocity
difference for each neighbor and Eq. (16) sums the alignment velocities for
each neighbor.

####  3.1.3 Wall collisions

To account for collisions at walls, the authors in vasarhelyi2018optimized
have proposed virtual "shill" agents at the walls which the actual agents can
try to align their velocities with. These shill agents have no gain and
therefore repulsion at walls takes place to the maximum extent
($c^{shill}=1$). This makes sense while flocking in confined environments,
because one of the primary goals is to avoid the wall at any cost. In our
research however, while seeking for non dominated set of solutions(ref.
section 6), we can characterise the _elasticity_ of the virtual geo-fence
using a shill gain ($c^{shill}$) parameter. The following equations are used
to find a shill velocity vector from each wall so as to align with it.
$\textbf{r}_{ci}$ is the relative position vector from the agent to the
arena’s center $\textbf{r}_{c}$. This vector is used to find the distances to
the walls in Eq. (18). Eq. (20) gives a $m\ \text{x}\ m$ sized matrix which
has the rows as the shill vector from each wall.

$\displaystyle\textbf{r}_{ci}$ $\displaystyle=\textbf{r}_{c}-\textbf{r}_{i}$
(17) $\displaystyle\textbf{r}_{s}^{mag}$
$\displaystyle=L_{c}/2-|\textbf{r}_{ci}|$ (18)
$\displaystyle\textbf{v}^{shillmax}_{i}$
$\displaystyle=D(\textbf{r}^{mag}_{s}-r_{0}^{shill},a^{shill},p^{shill})$ (19)
$\displaystyle V_{s}$
$\displaystyle=(v^{shill}.\frac{\textbf{r}_{ci}}{|\textbf{r}_{ci}|})\odot I$
(20) $\displaystyle\textbf{v}_{s}^{mag}$
$\displaystyle=\|V_{s}-\textbf{v}_{i}\|\ _{\top\textbf{v}^{shillmax}_{i}}$
(21) $\displaystyle V^{shill}$
$\displaystyle=c^{shill}.(\textbf{v}_{s}^{mag}-\textbf{v}^{shillmax}_{i}).\frac{V_{s}}{\textbf{v}_{s}^{mag}}$
(22) $\displaystyle\textbf{v}^{shill}_{i}$
$\displaystyle=\sum_{k=1}^{m}V^{shill}_{k}$ (23)

where,

$\textbf{r}_{c}$ Absolute position of the center of the arena

$L_{c}$ Side length of the arena

$\textbf{r}_{ci}$ Relative position of the center with respect to the agent

$p^{shill}\equiv$ Slope for the linear part of the decay curve (user dependent
parameter)

$a^{shill}\equiv$ Acceleration for the non-linear part of the decay curve
(user dependent parameter)

$c^{shill}\equiv$ Overall Gain for shilling alignment (user dependent
parameter)

$v^{shill}\equiv$ Speed of shilling agents (user dependent parameter)

$r_{0}^{shill}\ \equiv$ Alignment cutoff distance for maximum alignment (user
dependent parameter)

$\textbf{v}^{frict}_{i}\ \equiv$ Desired collective alignment vector

We assume a square geo-fence in our research but trivial modifications to Eq.
(17) and (18) can generalise it other shapes as well. Here, $I$ is the
identity matrix and $\odot$ is the Hadamard product. Eqs. (21) - (23) have the
same velocity alignment procedure done in section 3.1.2 but here it’s done for
each wall’s shill velocity instead of each agent.

The above three velocities (3.1.1 \- 3.1.3) along with the normalised flocking
velocity are summed up and normalised again to give the desired velocity for
the respective agent.

$\displaystyle\textbf{v}_{i}^{desired}$
$\displaystyle=\frac{\textbf{v}_{i}}{\|\textbf{v}_{i}\|}v^{flock}+\textbf{v}_{i}^{rep}+\textbf{v}_{i}^{frict}+\textbf{v}_{i}^{shill}$
(24) $\displaystyle\textbf{v}_{i}^{desired}$ $\displaystyle\longleftarrow
min\\{v^{max},\|\textbf{v}_{i}^{desired}\|\\}\
\frac{\textbf{v}_{i}^{desired}}{\|\textbf{v}_{i}^{desired}\|}$ (25)

Finally, the set of resulting 12 parameters to optimize is:

$x=\\{r_{0}^{sep},p^{rep},r_{0}^{frict},a^{frict},p^{frict},v^{frict},c^{frict},r_{0}^{shill},v^{shill},a^{shill},p^{shill},c^{shill}\\}$

###  3.2 Fitness functions

To measure the performance of one simulation run, order parameters are defined
and passed through transfer functions to get the
fitnessesvasarhelyi2018optimized .

$\displaystyle\begin{split}F^{speed}&=\textbf{F}_{1}(\phi^{vel},v^{flock},v^{tol})\\\
F^{coll}&=\textbf{F}_{3}(\phi^{coll},a^{tol})\\\
F^{wall}&=\textbf{F}_{2}(\phi^{wall},r^{tol})\\\
F^{corr}&=\Theta(\phi^{corr})\phi^{corr}\\\
F^{disc}&=\textbf{F}_{3}(\phi^{disc},N/5)\\\
F^{cluster}&=\textbf{F}_{3}(\phi^{cluster},N/5)\end{split}$ (26)

Here, the order parameters $\phi^{vel},\phi^{coll},\phi^{corr},\phi^{wall}$
and transfer functions $\textbf{F}_{1}$, $\textbf{F}_{2}$, $\textbf{F}_{3}$
are taken from vasarhelyi2018optimized  and $\Theta$ is a the heave-side step
function. Parameters $r^{tol},a^{tol},\text{and }v^{tol}$ are explained in
section 5. Order parameters for disconnected agents ($\phi^{disc}$) and the
minimum connected agents ($\phi^{cluster}$) are explained below. These
parameters are calculated locally in $r^{cluster}$ sized clusters:

$\displaystyle
r^{cluster}=r^{rep}+r^{frict}+\tilde{D}(v^{flock},a^{frict},p^{frict})$ (27)

$\tilde{D}$ is the braking distance r for which D(r, a, p) = v for any agent.

#### Disconnected agents

This parameter measures the average number of completely disconnected agents
throughout the simulation. Eq. (30) gives the number of agents within
$r^{cluster}$ distance of each agent at any given moment. Eq. (31) is then
used to determine the number of agents throughout the simulation with zero
connected agents, i.e. disconnected.

$\displaystyle\Theta(x)$ $\displaystyle=\begin{cases}1&\text{if $x\geq 0$}\\\
0&\text{if $x<0$}\end{cases}$ (28) $\displaystyle n_{i}^{cluster}(t)$
$\displaystyle=\sum_{j\neq i}^{N-1}\Theta(r^{cluster}-r_{ij}(t))$ (30)
$\displaystyle\phi^{disc}$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\sum_{i=1}^{N}\Theta(n_{i}^{cluster}(t)-1)$
(31)

#### Minimum connected agents

This parameter measures the minimum number of connected agents averaged
throughout the simulation and is therefore dependent on time. Since the drones
start at random positions, it was observed that keeping this parameter time
dependent instead of steady state (the global minimum throughout the
simulation) gave a better idea of the robustness of the communication graph
throughout the simulation.

$\displaystyle\phi^{cluster}(t)$
$\displaystyle=\frac{1}{T}\int_{0}^{T}min\\{n_{1}^{cluster},n_{2}^{cluster}...\
\ n_{i}^{cluster}\\}(t)$ (32)

$\forall\ i:1,2,...,N$

Finally, we need to optimize these fitness functions given in Eq. (26)
simultaneously To optimize the above fitness functions simultaneously, the six
objectives must be analyzed for correlations among them so that the system can
be represented with fewer objectives, preferably two. In the next section,
Principal Component Analysis (PCA) is used for the dimensionality reduction so
that the multi-objective optimizer NSGA-II can be used, effectively.

##  4 Dimensionality reduction using PCA

To reduce the number of objective functions, a data set of the six objectives
discussed in section 3 is collected. This data is just the result of 500
random simulations without any heuristic so as to cover the entire search
space. Note that the data used for PCA is for the fitness values after being
passed through the transfer functions. This can also be done directly on the
order parameters as well. Both processes would give different correlations
depending on the nature of the transfer function. We prefer the former method
as it gives a more accurate representation of the matrix components and the
exact fitnesses functions used for optimization. This data is used to create
the covariance matrix and principal components shown in section 2 which is
followed by a qualitative discussion on the correlations.

M

M

$\begin{array}[]{c}F^{wall}\\\ F^{speed}\\\ F^{corr}\\\ F^{coll}\\\
F^{disc}\\\ F^{cluster}\end{array}$

$\ \ 1^{st}$ component

$(w)$

$\left[\begin{array}[]{c}0.191\\\ 0.329\\\ -0.495\\\ -0.285\\\ -0.509\\\
-0.518\end{array}\right]$

Covariance matrix

$(K)$

$\left[\begin{array}[]{cccccc}1.002&0.2821&-0.172&-0.115&-0.094&-0.131\\\
0.282&1.002&-0.285&-0.283&-0.301&-0.259\\\
-0.172&-0.285&1.002&0.253&0.5815&0.682\\\
-0.115&-0.283&0.253&1.002&0.278&0.204\\\
-0.094&-0.301&0.581&0.278&1.002&0.748\\\
-0.131&-0.259&0.682&0.204&0.748&1.002\end{array}\right]$

Figure 1: Matrices obtained from Principal component analysis

In Fig. 1, the matrices obtained from the application of PCA on the objective
space is given. The matrices show some interesting results. Some insights are
discussed as follows:

$K_{23}$ is negative implying that a higher velocity doesn’t necessarily imply
higher correlation. This might be false in situations where the UAVs have very
high velocity magnitudes while travelling long distances or have a large turn
radius (as in the case of fixed wing drones). But in a confined environment,
to maintain correlation at the edges (where the flock gets broken up most),
the speeds must be reduced. This is also a consequence of a limited
acceleration which aligns with the actual physical systems.

The above statement regarding confined environments is also confirmed by
$K_{13}$. To maintain correlation at walls, the UAVs can either slow down or
skip the wall altogether. A combination of slowing down and breaching the geo-
fence makes the above movement the most efficient. Note that intuition would
suggest that as speed increases, it would be easier to decrease the wall
fitness as there is indeed a limited acceleration/deceleration available. Upon
running numerous simulations and making correlation matrices, it was found
that this is because $F_{wall}$ itself is time dependent. This means that the
fitness is inversely proportional to the amount of time frames that the drones
spend outside the wall. Since the goal of $F_{2}$ is to maintain correlation
and connectivity wherever possible at the expense of $\phi^{vel}$ and
$F^{wall}$, whenever the drones slow down they naturally spend more time
frames outside and turn slowly irrespective of the acceleration, This makes
$F^{wall}$ and $F^{speed}$ directly correlated with each other on average. The
elements $w_{11}$ and $w_{12}$ and $K_{12}$ are representative of this very
fact.

Drones naturally collide less with each other when their velocities are
aligned and they are well connected. This is because the time it takes for
velocity changes to travel throughout the communication network is much
lesser. Although, when this network is strongly connected, the agent has to
sum up through many velocity differences in it’s neighborhood. While this is
advantageous when the neighbors are moving in similar directions, it can be
detrimental when there is a lot of noise and the inter-agent velocity
differences point in different directions. As a result the summed up alignment
velocity for the concerned agent gets dampened by cancelling out. This results
in inter-agent friction and makes the entire flock sluggish (slow to react).
This is very clearly shown by the elements $w_{2}$ and $w_{3}$ which are
strongly uncorrelated. It is also worth pointing out that elements $w_{3}$
through $w_{6}$ are strongly correlated which confirms the association of
correlation and collisions with the communication network. The cluster
parameters for disconnection and minimum number of connected UAVs are strongly
correlated ($K_{65}$) as expected as they are both direct functions of the
communication network.

Using a single objective can result in loss of important information as the
final fitness is just the collective product or weighted sum. Particularly, in
noisy dynamical systems such as multi-agent robotics, efforts need to be made
to retain as much information as possible and use it intelligently to guide
the decision making process. We propose a multi-objective methodology for
optimization of the swarm’s fitness to tackle this problem.

The principal component ($w$) for the maximum variance captures all the above
relations and shows them how they relate with each other on average. The sign
of the elements indicate correlation which gives rise to the following
features/objectives:

$\displaystyle F_{1}$ $\displaystyle=F^{wall}\ .\ F^{speed}$ (33)
$\displaystyle F_{2}$ $\displaystyle=F^{corr}\ .\ F^{coll}\ .\ F^{disc}\ .\
F^{cluster}$ (34)

Unlike traditional PCA, we don’t use just the non-redundant objectives. Each
objective captures tangible physical information about the simulation and
therefore we multiply the two sets individually to retain that information and
also make it easier to draw a comparison with the single objective CMA-ES
optimizer as given in section 6.

It is worth mentioning here that the above correlation matrix is dependent on
the number of agents and the size of the confined arena. While some parameters
like the cluster connectivity and correlation still remain the same because
they are independent of the above parameters, a different non redundant set of
objectives was obtained upon changing the size of the geo-fence. The
simulations dictated that the same number of agents in a larger space took
more time to align with the shill agents due to the stronger inter-agent
alignment over long distances. The correlation matrix for the same is not
shown here for the sake of brevity. In vasarhelyi2018optimized , there is a
certain ambiguity in the size of the geo-fence. While the authors mentioned
that they used a side length of 250m for the square arena, the averaged
results on their open-source simulator were closer to the claimed ones when a
radius of 250m (or side length of 500m) was used for the arena. To make
comparisons easier we continue with the latter definition for our study as
well.

The reduced objectives are passed to the multi-objective solver NSGA-II
deb_fast_2002  and the results are summarized below.

##  5 Numerical Experiments

Figure 2: Class diagram

To test the propsed algorithm and for future work as well, a custom simulator
MOflock was created in the Python programming language. The ease of use in
setting up multiple processes, leveraging optimization and machine learning
libraries was a major influence in choosing Python. The simulator is highly
object oriented and modular. It has the drones abstracted at various levels
and allows experimenting with both single (Bot) and multiple collaborative
agents (CoBot). The class diagram for the same is given in Fig. 2. It was kept
in mind that error between RobotSimvasarhelyi2018optimized  and the current
work should remain under a threshold of 5-10%. The link of the repository for
the code is given in supplementary material (S1) and a screenshot of the
simulation is shown in Fig. 3. All the experiments are carried out with a
flocking velocity ($v^{flock}$) and maximum velocity ($v^{max}$) of 6 m/s but
no changes were made in the algorithm so as to disrupt the scalability in
velocity. Artificial GPS noise is added using the Brownian noise model used in
viragh2014flocking . Communication delays are integral to the result of
optimization as they simulate a kind of inertia at the walls and with
neighbours as well. Without these delays and noises it is observed that the
drones favour high gain and short range repulsion as opposed to the model
optima.

Figure 3: MOflock simulation screenshot

After analysing the covariance matrix for correlations (Section 4), the
objectives are combined accordingly and passed to the multi-objective
optimizer. A good multi-objective optimization algorithm should contain the
following characteristics:

  1. 1.

Guide the solutions to an optimal Pareto front

  2. 2.

Maintain solution diversity

NSGA-II is proven to be one of the best performing algorithms in this regard.
The ‘pymoo’ blank2020pymoo  python library is used for the same with default
parameters.

It is imperative to setup the optimization problem in such a way that there is
enough diversity in the search space so as to find "good enough" solutions
through heuristic methods. For the sake of exploration a test run is conducted
using the CMA-ES algorithm without any bounds on the parameters. The solution
for this setup revealed that the flock only moves in circles around the centre
and does not interact with the walls at all. While such a solution is
mathematically the most optimal, it does not encapsulate the physical
limitations and logical constraints on the variable bounds. This happens
because the correlation and wall fitnesses become abnormally high.

To avoid such false positives in the simulation, either explicitly known
bounds can be set on the variables which are realistic and relevant to the
physics of a UAV or another objective which maximises the search area covered
in minimum time can also be incorporated in the optimization process. For this
study, the former approach is used without any loss of generality. The bounds
used for the variables are shown in Table 2 and some miscellaneous simulation
parameters including certain tolerance parameters $r^{tol}$, $a^{tol}$, and
$v^{tol}$ for the transfer functions in section 3.2 are given in Table 2.
Appropriate values for these tolerance parameters promote better search of
solutions and gradient directions.

All the experiments were performed on a machine with the AMD Ryzen 7 4800H 16
core CPU and 16 GB of RAM. The results are reported in Section 6.

Table 1: Optimization bounds  | Lower bound | Upper bound
---|---|---
$r_{0}^{rep}$ | 30.8 | 51
$p^{rep}$ | 0.02 | 0.10
$r_{0}^{frict}$ | 58.5 | 100
$a^{frict}$ | 5.04 | 10.0
$p^{frict}$ | 0.38 | 9.67
$v^{frict}$ | 0.3 | 2.7
$c^{frict}$ | 0.03 | 0.22
$r_{0}^{shill}$ | -10 | 0
$v^{shill}$ | 10.0 | 15.0
$a^{shill}$ | 1.54 | 6.55
$p^{shill}$ | 0.48 | 9.96
$c^{shill}$ | 0.3 | 1
|  |

Table 2: Simulation parameters  Parameter | Value
---|---
$v^{flock}$ | 6 m/s
$v^{max}$ | 6 m/s
$t_{del}$ | 0.2 s
$N$ | 30
$L^{arena}$ | 500 m
$\sigma_{inner}$ | 0.005 $m^{2}/s^{2}$
$t_{del}$ | 1 s
$r^{coll}$ | 3 m
$v^{tol}$ | 3.75 m/s
$a^{tol}$ | 0.0003
$r^{tol}$ | 5 m
|

##  6 Results and Discussions

The results of the optimization procedure are analyzed and discussed in this
section.

Figure 4: Comparison of different configurations with respect to the optimal
Pareto front.

Fig. 4 shows statistical evaluations (mean $\pm$ standard deviation) of 100
simulations for different points. The blue and red curves are comparisons of
our simulator and RobotSim at the model optima for $v^{flock}=6m/s$.
$F(X_{a})|\emph{RobotSim}$ is the multi-objective fitness for the model optima
$X_{a}$ taken from vasarhelyi2018optimized  and evaluated on RobotSim itself.
$F(X_{a})|\emph{MOflock}$ is the fitness for $X_{a}$ evaluated on MOflock. As
targeted, the error on mean fitnesses between both simulator at the model
optima in vasarhelyi2018optimized  is 4.28%. $F^{CMA-
ES}(X_{opt})|\emph{MOflock}$ is the optimized fitness result on out simulator
using the CMA-ES algorithm. Note that this fitness wasn’t evaluated using a
multi-objective algorithm but was separated into $F_{1}$ and $F_{2}$ according
to Section 4. This is done so that comparisons can be drawn easily between the
single objective and multi-objective results. The Pareto front for the last
generation using the NSGA-II algorithm is also shown.
$F_{1}^{NSGA2},F_{2}^{NSGA2}|\emph{MOflock}$ are mean values for the extreme
points on this Pareto Front.

(a) Last 10 even Pareto fronts

(b) Point A (blue) and Point B (red)

Figure 5: Optimization Results

Since the single objective fitness is just a scalar product of all individual
fitnesses, it follows that neither of the six fitnesses can be close to zero
or even guaranteed to be maximum in case there exists a negative correlation
between some of them. As a result, when optimizing a single objective function
a ’best of both’ situation is sought after. In the case of multiple
conflicting objectives, however, this can be forgiven for better performance
on the separated fitness functions. This also explains why the CMA-ES point
lies around the knee of the Pareto fronts. It should be noted however, that
the CMA-ES optima on our simulator outperforms the Pareto front at it’s knee.
This is owed to the high degree of automation and robustness of the CMA-ES
algorithm.

While the user can now choose amongst any of the points depending on the
scenario and relative importance, there are two interesting points on the
optimal front corresponding to the extreme situations when either one of the
two solutions is compromised for the other. They are given by point A and B in
Fig. 4. The values of the variables and fitnesses at the above points are
summarised in Table 3. Fig. 5(a) shows the last 10 even Pareto fronts ranked
in descending order. The even ones were only chosen to display the spread and
convergence in a neat manner. A snapshot of the relevant simulations for both
points is also shown in Fig. 5(b) along with the graphs for their order
parameters in Fig. 6. They can be qualitatively understood as follows:

Point A: A weaker cluster dependent fitness shows that multiple clusters can
coexist in the same environment when correlation and speed is sacrificed.

Point B: Similarly, the other point clearly skips the geo-fence and/or slows
down to maintain a good correlation and make up for the damping caused by
inter-agent friction and pressure at the walls.

The generation of the above two points is a direct consequence of the physical
and environmental restrictions imposed on the swarm. The limited acceleration
does not allow the entire swarm to turn sharply without slowing down. The
confined walls don’t allow agents to flock together when moving at high speeds
without losing on some amount of correlation. These statements are a testament
to the complex dynamics that multi-agent systems exhibit. A video showing the
above interactions is given in supplementary material (S2). Better
mathematical formalism and high fidelity simulations can be developed to
realise such intertwined relationships.

(a) Point A

(b) Point B

Figure 6: Cumulative order parameters for points A and b.

The trend in the order parameters in Fig. 6 also confirms the elements of the
covariance matrix in Section 4. Note that the graph is scaled to the (0,1)
interval with the relevant maximum feasible values for each parameter and
cumulative values are shown for the curves.

Table 3: Optimization results | Point A | Point B
---|---|---
$F_{1}$($\mu\pm\sigma$) | $0.890\pm 0.039$ | $0.065\pm 0.039$
$F_{2}$($\mu\pm\sigma$) | $0.112\pm 0.09$ | $0.896\pm 0.179$
$r_{0}^{rep}$ | 33.69 | 33.45
$p^{rep}$ | 0.023 | 0.028
$r_{0}^{frict}$ | 59.26 | 58.95
$a^{frict}$ | 5.38 | 8.223
$p^{frict}$ | 4.62 | 2.67
$v^{frict}$ | 1.73 | 3.00
$c^{frict}$ | 0.035 | 1.84
$r_{0}^{shill}$ | -2.45 | -0.21
$v^{shill}$ | 12.93 | 12.93
$a^{shill}$ | 4.84 | 2.57
$p^{shill}$ | 4.83 | 1.30
$c^{shill}$ | 0.55 | 0.43
|  |

Further statistical analysis on the data from the optimization shows that
there is a lot of redundancy in the decision variables. The following
observations indicate this finding:

  * •

Even though points A and B are far apart on the Pareto front, their respective
parameters for repulsion are very similar.

  * •

It was observed that the right combination of $r_{0}^{shill}$ and $a^{shill}$
gives similar fitnesses and order parameters even with a constant shilling
velocity.

  * •

The introduced shill gain ($c^{shill}$) doesn’t take it’s maximum possible
value (1.0) even when seeking the best $F_{1}$ which is highly dependent on
this parameter.

Note that a full PCA correlation analysis on the decision variables can be
performed to confirm the above observation, and reduce the dimension of the
input space as well.

The above results are more consequential than just a Pareto front. Real life
missions and the inherent stochastic nature of the environment demands a range
of potential solutions from which a human in the loop can choose in an ad-hoc
manner. A typical mission profile consists of cruise, loiter, surveillance,
and occasionally a payload drop. A brief description of the use of the
practical applications of the Pareto optimal points are shown below.

  * •

Target search and loitering is a common phase in surveillance missions. A
snapshot of an extreme case where the target is located at a corner of the
geo-fence is shown in Fig. 7(a). The flock breaks at corners and walls to
loiter around the target. To make this observation mathematically sound,
another order parameter called $\phi^{target}$ is created.

$\displaystyle\textbf{x}^{COM}(t)$
$\displaystyle=\frac{\sum_{i=1}^{N}\textbf{r}_{i}(t)}{N}$ (35)
$\displaystyle\bar{d}^{target}(t)$
$\displaystyle=\|\textbf{x}^{target}-\textbf{x}^{COM}(t)\|$ (36)

where,

$\textbf{x}^{COM}(t)\equiv$ Center of mass of the swarm at time t

$\bar{d}^{target}(t)\equiv$ Mean distance to target (over all $N$ agents) at
time t

This parameter includes two performance measures- the closeness of the entire
flock to the target on average (Eq. 36) and the ‘Loiter Frequency ($\omega$)’.
This frequency measures how fast the flock can loiter around the target and
turn around as a whole. As opposed to the other parameters, the steady state
version of this parameter is measured. Since the motion is circular and
periodic, the time series is fit to a sinusoidal wave similar to an audio
signal.

$\displaystyle\phi^{target}(t)$
$\displaystyle=a.\sin(\omega.\bar{d}^{target}(t)+\psi)+c$ (38) $\displaystyle
F$ $\displaystyle=FFT(\phi^{target})$ (39)
$\displaystyle\overline{\bar{d}^{target}}$
$\displaystyle=\frac{1}{T}\sum_{t=0}^{T-1}\bar{d}^{target}(t)$ (40)
$\displaystyle a_{o}$
$\displaystyle=\sqrt{\frac{2}{T}\sum_{t=0}^{T-1}(\bar{d}^{target}(t)-\overline{\bar{d}^{target}})^{2}}$
(41) $\displaystyle f_{o}$ $\displaystyle=|f^{s}_{argmax(|A_{k}|)}|$ (42)
$\displaystyle\psi_{o}$ $\displaystyle=0$ (43) $\displaystyle c_{o}$
$\displaystyle=\overline{\bar{d}^{target}}$ (44) $\displaystyle
a,\omega,\psi,c$
$\displaystyle=LSF(\phi^{target},\bar{d}^{target},a_{o},f_{o},\psi_{o},c_{o})$
(45)

where,

$F\ \equiv$ Fourier transform output

$\overline{\bar{d}^{target}}\equiv$ Mean of the mean distance throughout the
simulation

$f^{s}\ \equiv$ Sample frequencies for the time series data

$\bar{d}^{target}\ \equiv$ Average distance throughout simulation

$f_{o}\ \equiv$ Initial guess of frequency for $\phi^{target}(t)$
corresponding to the maximum $F$

$\psi_{o}\ \equiv$ Initial guess of phase for $\phi^{target}(t)$

$a_{o}\ \equiv$ Initial guess of amplitude for $\phi^{target}(t)$

$c_{o}\ \equiv$ Initial guess of offset for $\phi^{target}(t)$

This is done by first getting an estimate of the initial coefficients, namely
amplitude ($a_{o}$), phase ($\psi_{o}$), offset ($c_{o}$), and frequency
($f_{o}=\omega_{o}/2\pi)$) via a Fast Fourier Transform ($FFT$) on the data
(Eq. (38)-(39)) and then passing this estimate for Least Squares curve Fit
represented by $LSF$ (Eq. (45)). The final order parameter is just the angular
frequency divided by the amplitude.

$\displaystyle F^{target}$ $\displaystyle=\omega/a$ (46)

The analysis shows that point A on the Pareto front has a lower loiter
frequency because of the extra inter-agent friction created to maintain the
flock correlation. Point B on the other hand has almost half the amplitude and
double the frequency because of the higher velocity, loosely correlated flock
with more collisions. These curves and an accompanying simulation screenshot
are shown in Fig. 7. A full video showing the target tracking and fitness
analysis is given in supplementary material S3.

(a) Target following simulation

(b) $\bar{d}^{target}(t)$ for Point A and B

Figure 7: Target order parameter comparison

  * •

There have been recent studies in which collisions are handled explicitly
through physical boundaries and mechanisms rather than by an explicit
algorithm mulgaonkar2017robust . The idea is to allow for some amount of
collisions as long as agility is maintained and the drones reach their target.
Point A on the front is akin to such a situation. The flock doesn’t give much
attention to inter-agent separation or correlation in local clusters. Rather,
the speed is given a higher priority. This is especially useful when tiny
drones need to overcome narrow passages and crevices without acting as a fully
connected flock but get through the region as fast as possible with each drone
acting for themselves.

  * •

Point B naturally resembles a good flock where connectivity and correlation is
concerned. The decentralised neighbor architecture makes the flock very
desirable where robustness and swarm health is an absolute requirement and the
entire swarm needs to travel long distances as a fully connected cluster.

While developing the methodology for this work, there were a number of nuanced
characteristics of collective behaviour noticed in the multi-agent
simulations. For instance, the parameters which characterise the swarm changed
drastically based on factors like communication delay and the arena size.
These two variables affect the swarm as a whole because any control action for
an agent close to the wall is propagated throughout the swarm with the
appropriate communication delay. Naturally, every PCA analysis with different
simulation parameters yielded unique objectives and therefore a different
Pareto front. The advantage of separating the objective function into multiple
grouped objectives is that global swarm behaviour can be controlled by
choosing a point on the Pareto front instead of tuning parameters manually or
running an offline optimization for each possible situation that the swarm
would encounter. It follows therefore, that the behaviour of the swarm can be
controlled by a supervisor with access to the appropriate Pareto front. This
problem of generalising a semi-autonomous swarm based on various scenarios is
often difficult to handle with just online learning algorithms. A compromise
between both, wherein, we can control large scale behaviour through multiple
objectives, and individual decision making through reinforcement learning can
be sought after to solve the generalisation problem.

##  7 Conclusion

In this paper, we proposed a methodology to address the problem of drone
flocking. First, a simulator with an integrated optimizer was designed to test
the algorithm. The decision variables which characterise the flocking
operators, and fitness functions which indicate the performance of the swarm
are defined. Then, to use the multi-objective optimizer effectively, the six
dimensional objective space is reduced to two dimensions using Principal
Component Analysis. The correlation analysis showed that fitness functions for
both speed and wall avoidance can be treated separately from the cohesive
movement of the entire flock. This process also gave insight into the various
complex relationships that multi-agent systems can exhibit. Further, the so
formed two objective optimization problem is optimized using NSGA-II and the
results are compared with the single objective CMA-ES optimization algorithm.
It is found that while CMA-ES performs better with respect to the knee of the
Pareto front, NSGA-II outperforms CMA-ES on the extreme points and offers an
entire range of solutions to choose from. The study also discussed the use
cases of such a Pareto front to guide the decision-making process in real-
world scenarios. Incorporating algorithms like Reinforcement Learning with the
proposed methodology can be future research agenda.

## Supplementary Material

  * •

S1: Repository: Multi-objective flocking simulator (https://github.com/nikhil-
sethi/MOflock)

  * •

S2: Comparison of the swarms from the two extreme points on the Pareto front
(https://youtu.be/MIgc80M1dHg

  * •

S3: Target following for the two extreme points on the Pareto front
(https://youtu.be/Fl71fg-DU_c)

## Data Availability

The datasets generated during and/or analysed during the current study are
available from the corresponding author on reasonable request.

## Declarations

This work was supported by the funding received in the international joint
research project titled “Design and Applications of Swarm Intelligence based
algorithms for drone swarm and COVID19 Spread Prediction” funded by Liverpool
Hope University UK. The authors declare that there is no conflict of interest.

## References

  *   * (1) Allison, C. & Hughes, C.  Bacterial swarming: an example of prokaryotic differentiation and multicellular behaviour.  _Science Progress_ 75 (298 Pt 3-4), 403–422 (1991) .
  * (2) Nagy, M., Akos, Z., Biro, D. & Vicsek, T.  Hierarchical group dynamics in pigeon flocks.  _arXiv:1010.5394 [physics, q-bio]_ (2010).  URL http://arxiv.org/abs/1010.5394.  10.1038/nature08891, arXiv: 1010.5394 .
  * (3) Czaczkes, T. J., Grüter, C. & Ratnieks, F. L.  Trail pheromones: An integrative view of their role in social insect colony organization.  _Annual Review of Entomology_ 60 (1), 581–599 (2015).  10.1146/annurev-ento-010814-020627, pMID: 25386724 .
  * (4) Ron, J. E., Pinkoviezky, I., Fonio, E., Feinerman, O. & Gov, N. S.  Bi-stability in cooperative transport by ants in the presence of obstacles.  _PLoS computational biology_ 14 (5), e1006068 (2018) .
  * (5) Saffre, F., Hildmann, H. & Karvonen, H.  _The design challenges of drone swarm control_ , 408–426 (Springer, 2021).
  * (6) Coppola, M., McGuire, K. N., De Wagter, C. & de Croon, G. C.  A survey on swarming with micro air vehicles: Fundamental challenges and constraints.  _Frontiers in Robotics and AI_ 7, 18 (2020) .
  * (7) Brust, M. R. _et al._ Defending against Intrusion of Malicious UAVs with Networked UAV Defense Swarms.  _arXiv:1808.06900 [cs]_ (2018).  URL http://arxiv.org/abs/1808.06900.  10.1109/LCN.Workshops.2017.71, arXiv: 1808.06900 .
  * (8) Abraham, L. _et al._ _Swarm robotics in disaster management_ , 1–5 (IEEE, 2019).
  * (9) Tosato, P. _et al._ _An Autonomous Swarm of Drones for Industrial Gas Sensing Applications_ , 1–6 (2019).
  * (10) Krizmancic, M., Arbanas, B., Petrovic, T., Petric, F. & Bogdan, S.  Cooperative Aerial-Ground Multi-Robot System for Automated Construction Tasks.  _IEEE Robotics and Automation Letters_ 5 (2), 798–805 (2020).  10.1109/LRA.2020.2965855, conference Name: IEEE Robotics and Automation Letters .
  * (11) Ju, C. & Son, H. I.  Multiple uav systems for agricultural applications: control, implementation, and evaluation.  _Electronics_ 7 (9), 162 (2018) .
  * (12) Hauert, S. _et al._ _Reynolds flocking in reality with fixed-wing robots: communication range vs. maximum turning rate_ , 5015–5020 (IEEE, 2011).
  * (13) Dewi, M., Hariadi, M. & Purnomo, M. H.  _Simulating the movement of the crowd in an environment using flocking_ , 186–191 (IEEE, Bandung, West Java, Indonesia, 2011).  URL http://ieeexplore.ieee.org/document/6108638/.
  * (14) Moere, A.  _Time-Varying Data Visualization Using Information Flocking Boids_ , 97–104 (2004).  ISSN: 1522-404X.
  * (15) Vásárhelyi, G. _et al._ Optimized flocking of autonomous drones in confined environments.  _Science Robotics_ 3 (20) (2018) .
  * (16) Márquez-Vega, L. A., Aguilera-Ruiz, M. & Torres-Treviño, L. M.  Multi-objective optimization of a quadrotor flock performing target zone search.  _Swarm and Evolutionary Computation_ 60, 100733 (2021) .
  * (17) Shah, S., Dey, D., Lovett, C. & Kapoor, A.  _Airsim: High-fidelity visual and physical simulation for autonomous vehicles_ , 621–635 (Springer, 2018).
  * (18) Soria, E., Schiano, F. & Floreano, D.  Swarmlab: A matlab drone swarm simulator.  _arXiv preprint arXiv:2005.02769_ (2020) .
  * (19) Kate, B., Waterman, J., Dantu, K. & Welsh, M.  _Simbeeotic: A simulator and testbed for micro-aerial vehicle swarm experiments_ , 49–60 (2012).
  * (20) Abson, D. J., Dougill, A. J. & Stringer, L. C.  Using principal component analysis for information-rich socio-ecological vulnerability mapping in southern africa.  _Applied Geography_ 35 (1-2), 515–524 (2012) .
  * (21) Pozo, C., Ruiz-Femenia, R., Caballero, J., Guillén-Gosálbez, G. & Jiménez, L.  On the use of principal component analysis for reducing the number of environmental objectives in multi-objective optimization: Application to the design of chemical supply chains.  _Chemical Engineering Science_ 69 (1), 146–158 (2012) .
  * (22) Loeffler, H. H. & Kitao, A.  Collective dynamics of periplasmic glutamine binding protein upon domain closure.  _Biophysical journal_ 97 (9), 2541–2549 (2009) .
  * (23) Deb, K. & Saxena, D.  _Searching for pareto-optimal solutions through dimensionality reduction for certain large-dimensional multi-objective optimization problems_ , 3352–3360 (2006).
  * (24) Deb, K., Pratap, A., Agarwal, S. & Meyarivan, T.  A fast and elitist multiobjective genetic algorithm: NSGA-II.  _IEEE Transactions on Evolutionary Computation_ 6 (2), 182–197 (2002).  URL http://ieeexplore.ieee.org/document/996017/.  10.1109/4235.996017 .
  * (25) Reynolds, C. W.  Flocks, herds and schools: A distributed behavioral model.  _SIGGRAPH Comput. Graph._ 21 (4), 25–34 (1987).  URL https://doi.org/10.1145/37402.37406.  10.1145/37402.37406 .
  * (26) Virágh, C. _et al._ Flocking algorithm for autonomous flying robots.  _Bioinspiration & biomimetics_ 9 (2), 025012 (2014) .
  * (27) Blank, J. & Deb, K.  pymoo: Multi-objective optimization in python.  _IEEE Access_ 8, 89497–89509 (2020) .
  * (28) Mulgaonkar, Y., Makineni, A., Guerrero-Bonilla, L. & Kumar, V.  Robust aerial robot swarms without collision avoidance.  _IEEE Robotics and Automation Letters_ 3 (1), 596–603 (2017) .