# A spectral-neighbour representation for vector fields: machine-learning
potentials including spin

M. Domina  M. Cobelli  S. Sanvito  School of Physics and CRANN Institute,
Trinity College Dublin, Ireland

###### Abstract

We introduce a translational and rotational invariant local representation for
vector fields, which can be employed in the construction of machine-learning
energy models of solids and molecules. This allows us to describe, on the same
footing, the energy fluctuations due to the atomic motion, the longitudinal
and transverse excitations of the vector field, and their mutual interplay.
The formalism can then be applied to physical systems where the total energy
is determined by a vector density, as in the case of magnetism. Our
representation is constructed over the power spectrum of the combined angular
momentum describing the local atomic positions and the vector field, and can
be used in conjunction with different machine-learning schemes and data taken
from accurate ab initio electronic structure theories. We demonstrate the
descriptive power of our representation for a range of classical spin
Hamiltonian and machine-learning algorithms. In particular, we construct
energy models based on both linear Ridge regression, as in conventional
spectral neighbour analysis potentials, and gaussian approximation. These are
both built to represent a Heisenberg-type Hamiltonian including a longitudinal
energy term and spin-lattice coupling.

Machine learning, Descriptors, Magnetism

##  I Introduction

The modelling of the structural, electronic and magnetic properties of
materials at finite temperature requires the exploration of complex energy
surfaces, a task that is usually performed in the configuration space or
through extensive time-dependent simulations. The gold standard is set by ab
initio methods, where one solves directly an electronic problem, for instance
through density functional theory. The accuracy of the method is then given by
the accuracy of the underline electronic structure theory, which often can be
determined by the level of approximation taken. Notably, also the
computational overheads are set by the electronic structure theory, so that
there is always a trade off between the accuracy of the prediction, the
duration of the simulation and the maximum size of the system to simulate.
Real-time time-dependent simulations of purely electronic quantities, such as
the dynamics of spins [1, 2, 3], are limited to a few atoms and a few hundreds
of femtoseconds, while ab initio molecular dynamics (MD) simulations can reach
well within the picosecond range and may involve several hundreds of atoms.

The general strategy for extending the range of dynamics simulations across
both the time and length scale is to abandon completely the ab initio
description and replace the solution of the electronic problem with some
parametric functions, constructed to reproduce the ab-initio potential energy
surface, namely classical force fields [4]. In their most canonical form one
describes the interaction among the ions by introducing energy contributions
that account for the various physical forces at play (covalent bond,
dispersive forces, etc.). Once the force field is defined, only the atomic
positions determine the total energy. A similar approach has been recently
introduced for spin-dynamics. In this case one associates to each atom a
classical spin vector, $\mathbf{S}_{i}$, so that the total energy is defined
over a continuous vector field with values at the atomic positions [5, 6]. The
total energy then takes the form of a classical Heisenberg model and may
include both anisotropy and friction terms. Furthermore, the formalism can be
extended to include spin-lattice coupling [7], longitudinal spin fluctuations
[8], and possibly the effects of a spin current [9].

In general, force fields constructed in this way have two main drawbacks. On
the one hand, their accuracy is significantly inferior to that of an ab initio
electronic structure theory, although this varies depending on the class of
compounds one wants to study. On the other hand, they tend to be specific to
the particular type of bond they describe. In addition, spin-type force fields
may not be able to describe entire excitations types. For instance, magnetic
Stoner excitations are not part of the spectrum of a classical Heisenberg
model.

Recently, a new class of force fields, named machine-learning force fields
(MLFFs), have been shown to solve both the accuracy and specificity issues.
The general philosophy of MLFFs is quite different from that of their
classical counterparts, since one does not pretend to construct an energy
function of the atomic coordinates with terms baring a physical
interpretation, but instead tries to reproduce extremely accurately the ab
initio potential energy surface. MLFFs comprise of two parts, an abstract
representation of the atom density distribution [10] and a machine-learning
model that correlates such representation to the system total energy. For the
two parts to work together the representation should be translational, atom-
permutational and rotational invariant. The first two conditions are usually
met by local representations (also known as atomic-neighbour descriptions),
where the energy is expressed in terms of atomic contributions, while the last
condition is constructed in. Symmetry functions [11] and bispectra [12] are
two examples of such rotational invariant local representations.

Importantly, the spin degrees of freedom are not explicitly included in the
representation, which typically describes an atomic density only. For this
reason MLFFs are currently unable to describe the energy difference between
inequivalent magnetic phases, for instance between a ferromagnetic and an
antiferromagnetic ground state, unless the different phases are also
associated to different structures. This limitation can be addressed by
combining the MLFF with a classic spin Hamiltonian in order to create a model
able to predict energies dependent on both the atomic positions and the spin
order. This route was recently explored by Nikolov and co-workers [13], who
equipped a spectral neighbour analysis potential (SNAP) with a classical
Heisenberg Hamiltonian to compute thermodynamical properties (e.g. the Curie
temperature) of $\alpha$-Fe.

Another possible strategy is to make the MLFF aware of the spin configuration
by including such information in the input features of the model. This is a
non trivial task, since such features need to retain the aforementioned
symmetries in order for the model to perform well. A recent attempt along this
direction consists in the introduction of a novel definition of symmetry
functions carrying spin information [14]. Such reformulation was constructed
for the spin-collinear case and for fixed spin magnitude (no longitudinal spin
information is available). Similarly, the atomic-cluster expansion method was
recently extended to vector fields [15], in a way that enables the description
of non-collinear spin configurations. Notably, by design both methods require
a large number of features for an accurate description of the magnetic
environment, making them particularly data hungry. Thus, a compact
representation describing a vector field, hence able to compute atomic and
magnetic excitations on the same footing, still remains at large. Our paper
aims to fill the gap.

Here, we propose a new local representation for vector fields, which can be
used with either linear and non-linear machine-learning models. This is based
on the power spectrum of the combined angular momentum describing the local
atomic positions and the vector field. The representation is rotationally
invariant and can be further generalised to tensorial densities. In order to
test its descriptive power, such representation is combined with either linear
Ridge regression and gaussian approximation to construct MLFFs describing the
potential energy surface of a Heisenberg model with longitudinal fluctuations
and spin-lattice coupling. Our results show that extremely accurate energy
predictions can be obtained with a rather moderate number of training data.

##  II Methods

###  II.1 Density for a vector field

The starting point of any atomic-neighbour description of solids and molecules
consists in defining the local particle density associated to the $i$-th atom,

$\rho_{i}(\mathbf{r})=\sum_{r_{ai}<r_{\text{cut}}}w_{a}h_{a}(\bm{r}-\bm{r}_{ai})\>.$
(1)

Here, $\bm{r}_{ai}=\bm{r}_{a}-\bm{r}_{i}$ is the distance between the atoms at
the position $\bm{r}_{a}$ and $\bm{r}_{i}$, $r_{ai}=|{\bm{r}_{ai}}|$, so that
the coordinates of the $i$-th atom define the origin of the local reference
frame. The sum in Eq. (1) runs over all the atoms inside a sphere of radius
$r_{\text{cut}}$ with center at $\bm{r}_{i}$, while $w_{a}$ are weights
usually associated to the atomic species of the $a$-th atom. In this expansion
$h_{a}$ is a localisation function, such a Gaussian or a Dirac-delta, centred
at $\bm{r}_{a}$, whose specific shape, in general, can depend on the $a$-th
atom type. Atomic-neighbour descriptions are then constructed by defining
rotationally invariant combinations of the coefficients of expansion of
$h_{a}$ over an appropriate local basis [10].

In the same spirit, we can now define a local _vector_ density,
$\bm{\rho}(\mathbf{r})$, through Eq. (1) by associating a vector $\bm{v}_{a}$
to each position $\bm{r}_{a}$, namely

$\bm{\rho}(\mathbf{r})=\sum_{a}w_{a}h_{a}(\bm{r}-\bm{r}_{a})\bm{v}_{a}\>,$ (2)

where for simplicity we have dropped the index $i$. In this formulation
$\bm{v}_{a}$ may, for instance, represent the local moment of the ions in a
magnetic compound, so that $\bm{\rho}(\mathbf{r})$ describes the local
magnetisation field. In Eq (2) the vector $\bm{v}_{a}$ is defined through its
cartesian components

$\bm{v}_{a}=v_{a,x}\hat{\bm{e}}_{x}+v_{a,y}\hat{\bm{e}}_{y}+v_{a,z}\hat{\bm{e}}_{z}=\sum_{i=x,y,z}v_{a,i}\hat{\bm{e}}_{i}\>,$
(3)

with $\hat{\bm{e}}_{i}$ being the unit vector along $i=x,y,z$. However, it is
convenient to replace the decomposition of Eq. (3) with one using the
spherical versors [16] ,

$\hat{\bm{e}}_{\pm 1}=\mp\dfrac{1}{\sqrt{2}}(\hat{\bm{e}}_{x}\pm
i\hat{\bm{e}}_{y})\quad\text{and}\quad\hat{\bm{e}}_{0}=\hat{\bm{e}}_{z}\>,$
(4)

so $\bm{v}_{a}$ becomes,

$\bm{v}_{a}=\sum_{q=0,\pm
1}v_{a,q}\hat{\bm{e}}_{q}\quad\text{with}\quad\begin{dcases}v_{a,\pm
1}=\mp\dfrac{1}{\sqrt{2}}(v_{a,x}\mp iv_{a,y})\>,\\\
v_{a,0}=v_{a,z}\>.\end{dcases}$ (5)

This decomposition is a particular case of the more general one for a tensor
of order $\lambda$ in its irreducible spherical components. Therefore, the
spherical components $v_{a,q}$ transforms under rotations as the spherical
harmonic $Y^{q}_{1}$ [16, 17].

In order to construct covariant descriptors for the local vector density of
Eq. (2), one first needs to expand the spatial part, $h_{a}$, over an
orthonormal radial basis. Here we use the product between a radial basis,
$R_{nl}(r)$, and the three-dimensional spherical harmonics,
$Y_{l}^{m}(\hat{\bm{r}})$, where as usual $n$, $l$ and $m$ are, respectively,
the principal, the angular momentum and the third component of the angular
momentum quantum numbers. The local vector density then becomes,

$\bm{\rho}(\bm{r})\simeq\sum_{n=0}^{n_{\text{max}}}\sum_{l=0}^{n}\sum_{m=-l}^{l}\sum_{q=0,\pm
1}c_{nlmq}R_{nl}(r)Y_{l}^{m}(\hat{\bm{r}})\hat{\bm{e}}_{q}\>,$ (6)

where the equality holds only for a complete basis but not for the one
truncated at $n_{\text{max}}$, and where the coefficients of expansion are
calculated as

$c_{nlmq}=\sum_{a}w_{a}v_{a,q}\int\differential
r\differential\hat{\bm{r}}\>r^{2}R_{nl}(r){Y_{l}^{m}}^{*}(\hat{\bm{r}})h_{a}(\bm{r}-\bm{r}_{a})\>.$
(7)

In this work we choose the radial-basis set introduced for the Spherical-
Bessel descriptors [18], namely $R_{nl}(r)=g_{n-l,l}(r)$. These are
orthonormal on the sphere and smoothly vanish at the cut-off radius. Note that
the choice of the radial basis set is not unique or crucial and alternative
basis can be selected. The only practical criterion is that they should
approximate completeness with a relative small number of basis functions,
namely the convergence must be rapid. Note also that, in what follows, we will
always assume orthonormality within the radial basis set, although the
expressions derived could be easily generalised to the non-orthogonal case.

The use of the spherical components of a tensor has been already exploited in
the construction of covariant kernels for vectorial and a tensorial properties
related to an atomic environment [19, 20, 21]. Here, we follow the same basic
idea and formalism, which stems from the decomposition of a tensorial object
into spherical components. The main difference is that we cannot just
associate a tensorial object to our density, since the density itself is the
vectorial field. As such, we will avail of the same concepts and methods that
are typically used to describe the coupling of angular momenta and vector
fields in atomic physics [16]. Another difference with existing literature is
that our primary target is the construction of invariant quantities instead of
covariant ones. A similar goal has been already pursued in Ref. [15], where a
vector field was described using the same strategy employed in dealing with
the atomic positions. More explicitly, in [15] the magnetic vectors are
encoded in Dirac’s delta distributions, which are then expanded on a suitable
basis set and coupled with the analogous expansion arising from the positions
of the atoms. Here, however we will preserve the vectorial nature of the field
at each step in the derivation, a strategy that results in a simpler coupling
scheme, as it will be shown in detail below. With this in mind, we will first
proceed with deriving an invariant power spectrum for the vector density of
Eq. (2), and then we will implement a linear regression to fit such power
spectrum. This second step is similar to what is commonly done with the
formalism of the spectral neighbour analysis potentials (SNAPs) [22].

For the remaining of the paper we will develop the formalism by using a
conventional Dirac notation, which allows one to appreciate better the
structure of the representation. In fact, as shown in reference [23], the
Dirac notation gives us a natural tool for dealing with local atomic
densities. In this way the expansion of the vector density defined in Eq. (6)
can be written in a compact form as

$\ket{\bm{\rho}}=\sum_{nlmq}c_{nlmq}\ket{nlmq}\>,$ (8)

where

$\innerproduct{\bm{r}}{nlm}=g_{n-l,l}(r)Y^{m}_{l}(\hat{\bm{r}}),\quad\text{and}\quad\ket{q}\equiv\hat{\bm{e}}_{q}\>.$
(9)

We can then express all the relevant quantities over the the basis,
$\ket{nl1JM}$, of the combined angular momenta, $\bm{L}+\bm{1}=\bm{J}$, by
using the standard addition scheme [17],

$\ket{nl1JM}=\sum_{m=-l}^{l}\sum_{q=-1}^{1}C^{JM}_{lm1q}\ket{nlmq},$ (10)

where $C^{JM}_{lm1q}$ are the Clebsch-Gordan coefficients. As usual, $J$ and
$M$, are the quantum numbers for the total angular momentum and its
projection, while ‘1’ refers to the angular momentum of $\bm{1}$, describing
the vector nature of the field. Hence, we have

$\absolutevalue{l-1}\leq J\leq{l+1}\quad\text{and}\quad-J\leq M\leq J\>.$ (11)

The states $\ket{nlJM}$, when projected over the position representation, are
the products of _vector_ spherical harmonics and radial functions. By
inverting equation (10),

$\ket{nlmq}=\sum_{JM}C^{JM}_{lm1q}\ket{nlJM}\>,$ (12)

we can write the vector density as

$\ket{\bm{\rho}}=\sum_{nlJM}u_{nlJM}\ket{nlJM},$ (13)

with

$\displaystyle u_{nlJM}$ $\displaystyle=$
$\displaystyle\innerproduct{nlJM}{\bm{\rho}}=\sum_{mq}C^{JM}_{lm1q}\innerproduct{nlmq}{\bm{\rho}}$
(14) $\displaystyle=$ $\displaystyle\sum_{mq}C^{JM}_{lm1q}c_{nlmq}.$

The last equality follows from the fact that the Clebsch-Gordan coefficient
are real and from the orthonormality of the $g_{n-l,l}$ functions. The form of
$\ket{\bm{\rho}}$ given by Eq. (13) contains the expansion of the vector
density over the combined angular-momenta basis. It should be noted that the
presence of the Clebsch-Gordan coefficients imposes that the values of $J$ and
$M$ must satisfy the conditions (11). The Clebsch-Gordan coefficients impose
also that the non-zero terms in the double sum of (14) are such that $M=m+q$.

Figure 1:  (Color online) The invariant power spectrum (right-hand side panel)
for a pyramid-shaped set of atoms, where a vector field is associated to each
atom (left-hand side image). The power spectrum is evaluated with respect to
the shown centre of rotation (red dot) and it is invariant against the
simultaneous rotation of the system and the vectorial field around such
centre. We chose the value $n_{\text{max}}=6$ for this example. The number of
independent power-spectrum components is given by
$(n_{\text{max}}+1)(3n_{\text{max}}+2)/2=70$. The $x$-axis labels the power-
spectrum components, while the $y$-axis shows their actual value.

###  II.2 Invariant power spectrum for a vector field

In this section we are going to introduce an invariant power spectrum for the
density given in Eq. (13). Let us initially restrict our formulation to the
case in which we ignore the atom at the origin of the local reference frame,
namely we assume that the magnitude of the vector field is zero at the origin.
One way to obtain the power spectrum, $p_{nlJ}$, is through the construction
of the following inner product,

$\innerproduct{\bm{\rho}}{\bm{\rho}}=\sum_{nlJ}p_{nlJ}\>,$ (15)

which explicitly reads

$p_{nlJ}=\sum_{M}\absolutevalue{u_{nlJM}}^{2}\>,$ (16)

where we have used the orthogonality of the $\ket{nlJM}$ basis. Since the
vectors $\ket{nlJM}$ correspond to the coupled angular momenta, they transform
under system rotation, $\hat{R}$, as the spherical harmonics $Y^{M}_{J}$,
namely

$\hat{R}Y^{M}_{J}=\sum_{M^{\prime}}D^{J}_{M^{\prime}M}(\hat{R})Y^{M^{\prime}}_{J}\>,$
(17)

where $D^{J}_{MM^{\prime}}(\hat{R})$ is the Wigner $D$-matrix associated to
the rotation $\hat{R}$. It must be noted that, when one considers the original
angular momentum basis, the rotation $\hat{R}$ appears as a simultaneous
rotation of both the positions of the atoms and the vector field. By applying
this rotation to the density in Eq. (13) we obtain

$\displaystyle\hat{R}\ket{\rho}$ $\displaystyle=$
$\displaystyle\sum_{nlJM}u_{nlJM}\hat{R}\ket{nlJM}=$ (18) $\displaystyle=$
$\displaystyle\sum_{nlJM}u_{nlJM}\sum_{M^{\prime}}D^{J}_{M^{\prime}M}(\hat{R})\ket{nlJM^{\prime}}=$
$\displaystyle=$
$\displaystyle\sum_{nlJM^{\prime}}u_{nlJM^{\prime}}^{\prime}\ket{nlJM^{\prime}}\>,$

from which we can infer the transformation rule for the expansion coefficients

$\hat{R}:u_{nlJM}\rightarrow\sum_{M^{\prime}}u_{nlJM^{\prime}}D^{J}_{MM^{\prime}}(\hat{R})\>.$
(19)

Therefore under rotation $\hat{R}$ the power spectrum,
$p_{nlJ}=\sum_{M}u_{nlJM}^{*}u_{nlJM}$, transforms as

$\displaystyle\hat{R}:p_{nlJ}$ $\displaystyle\rightarrow$
$\displaystyle\sum_{M^{\prime}M^{\prime\prime}}u^{*}_{nlJM^{\prime}}u_{nlJM^{\prime\prime}}\underbrace{\sum_{M}\left(D^{J}_{MM^{\prime}}\right)^{*}D^{J}_{MM^{\prime\prime}}}_{=\delta_{M^{\prime}M^{\prime\prime}}}$
(20)
$\displaystyle=\sum_{M^{\prime}}u_{nlJM^{\prime}}^{*}u_{nlJM^{\prime}}=p_{nlJ}\>,$

where we have used the unitarity of the Wigner $D$-matrices [note that we have
shorten the notation $D^{J}_{MM^{\prime}}(\hat{R})$ into
$D^{J}_{MM^{\prime}}$]. This proves that the power spectrum obtained from Eq.
(16) is rotationally invariant for simultaneous rotations of the atomic
positions and the vector field. Such invariance is shown numerically in Fig.
1, where the power spectra computed for different rotations are shown to
perfectly overlap. In the Appendix, we will briefly discuss the generalized
power spectrum connecting different radial channels. Furthermore, we will also
extend our construction to the more general case of a tensorial density.

If we now take a localization function of the form

$h_{a}(\bm{r}-\bm{r}_{i})=\delta(\bm{r}-\bm{r}_{a})\>,$ (21)

namely a Dirac-delta function centered on the $a$-th atom, then the local
vector density reads

$\bm{\rho}(\bm{r})=\sum_{a}w_{a}\delta(\bm{r}-\bm{r}_{a})\bm{v}_{a}\>.$ (22)

In this case the expansion coefficients of Eq. (14) are readily evaluated by
using Eq. (7) as

$u_{nlJM}=\sum_{a}w_{a}g_{n-l,l}(r_{a})\sum_{mq}C^{JM}_{lm1q}Y^{m*}_{l}(\hat{\bm{r}}_{a})v_{a,q}\>.$
(23)

In what follows, we will use this expression to explicitly evaluate the power
spectrum.

If we now consider a vector field having a vector-baring atom at the origin,
it can be proven (see the Appendix) that the power spectrum is not generally
invariant under rotations. We can interpret this rotational-symmetry breaking
by noting that a vector field at the origin introduces an inner preferential
direction for the local reference frame. One pragmatic solution to recover the
invariance is to always rotate the system so that the vector field at the
origin points along the $z$-axis. After such alignment we obtain a power
spectrum, which is invariant under rotations around the $z$-axis. Another
possible solution, is to choose a suitable radial basis set so that all the
non-invariant terms are automatically removed. As shown in the Appendix, we
proved that the Spherical Bessel functions have this property. In the
following, we will always consider power spectra with a central atom. As an
example of its explicit evaluation, it is useful to obtain the complete
expression for the $p_{n0J}$ power spectra. The component $u_{n0JM}$ is
proportional to

$u_{n0JM}\propto\delta_{J1}\sum_{a}w_{a}g_{n0}(r_{a})v_{a,M}\>,$ (24)

where we used the equalities $q=M$ and $J=1$ enforced through the Clebsh-
Gordan coefficients $C^{JM}_{0010}$, and we did not carry over the spherical
harmonics values and the Clebsh-Gordan coefficients which, in this case, are
unessential constants. The power spectrum is then proportional to

$\displaystyle p_{n0J}$ $\displaystyle\propto$
$\displaystyle\delta_{J1}\sum_{ab}w_{a}w_{b}g_{n0}(r_{a})g_{n0}(r_{b})\sum_{M}v_{a,M}v_{b,M}^{*}\>,$
(25) $\displaystyle\propto$
$\displaystyle\delta_{J1}\sum_{ab}w_{a}w_{b}g_{n0}(r_{a})g_{n0}(r_{b})\bm{v}_{a}\cdot\bm{v}_{b}\>.$

It is interesting to note that, if the vector field is made of local spins,
the $p_{n0J}$ power spectrum component will have a structure similar to that
of an Heisenberg model, namely it depends on the inner product between the
spins with distances-only dependent coefficients.

Having derived a set of invariant power spectrum also for the case of an atom
at the origin, we can now introduce the models used to test our formalism and
the machine-learning scheme that implements the power spectrum.

###  II.3 Training a machine-learning model

The machine-learning model used here is a linear regression constructed over
the power spectrum $\\{p_{nlJ}\\}$, following the same philosophy of SNAP
[22]. Thus, let us assume to have a system of $N$ atoms, each one of them
bearing a local spin. Let us also define the power spectrum vector, $\bm{p}$,
as the one dimensional vector, whose entries are the $p_{nlJ}$ components
describing a specific local neighborhood. Specifically, the $i$-th power
spectrum vector $\bm{p}^{(i)}$ is the vector obtained by centring the local
reference frame on the $i$-th atom, and then by evaluating the $i$-th power
spectrum set $\\{p_{nlJ}^{(i)}\\}$ with respect to that frame. Thus, given our
$N$-atoms system, we obtain $N$ power spectrum vectors, $\bm{p}^{(i)}$.

Our main working hypothesis is that the energy, or any other quantity that we
wish to represent, can be written as the sum of short-ranged contributions,
$\varepsilon_{i}(\bm{q}^{i})$, located on each atoms [12], namely

$E=\sum_{i}\varepsilon_{i}(\bm{q}^{i})\>,$ (26)

where $\bm{q}^{i}$ is a vector describing the local environment of the $i$-th
atom. Then, following the same idea behind the SNAP [22], we further assume
that the power spectrum vectors $\bm{p}^{(i)}$ form a suitable set of
descriptors to represent such decomposition, so that the local energies can be
expressed as a linear combination of power spectrum vectors,

$E\simeq\bm{\theta}\cdot\sum_{i}\bm{p}^{(i)}=\sum_{nlJ}\theta_{nlJ}\sum_{i}p_{nlJ}^{(i)}\>.$
(27)

Here $\\{\theta_{nlJ}\\}$ is an appropriate set of weights. The validity of
these assumptions cannot be determined from the outset and must be tested on a
case-by-case base. Within our formalism, the power spectrum vectors can be
seen as the descriptors of a linear regression problem, where the target is
the energy of the system.

In what follows we will first calculate the energies of several atomic and
spin configurations obtained by displacing the position of the atoms and the
direction and magnitude of the magnetic vectors. Then, we will evaluate the
power spectrum vectors for each atom and for each of the configurations
considered. Finally, we will train a Ridge regression and optimize the weights
vector, $\bm{\theta}$, to predict the total energy. This will allow us to
investigate the descriptive power of our vectorial representation and of the
full method proposed. In the next section we describe the different models
investigated.

Figure 2: The physical system investigated in the present work: a tetrahedral
cluster of bcc iron. The system is made of $7$ stacked square-shaped layers of
$5$ or $6$ atoms per side, as shown in the insert. The total number of atoms
is 219.

###  II.4 The physical system investigated

In this work we consider a rectangular cluster of sites arranged over a bcc
lattice, containing 219 atoms of the same species. The cell is rectangular
with a six-atom wide square base and no periodic boundary conditions, as shown
in Fig. 2. Each atom bares a local spin and can be displaced from the ideal
high-symmetry bcc site. The training data, namely the atomic and spin
configurations and their associated total energies, may come from a suitable
total energy theory. This is usual some ab initio method such as spin-
polarised density functional theory or a quantum-chemistry wave-function
scheme. Since the generation of such dataset is rather time consuming, and our
objective here is simply that of introducing our vector field representation,
we use instead a range of analytical energy models.

In particular, we assume that the total energy is determined by the
Hamiltonian,

$H=H_{\mathrm{H}}+H_{\mathrm{L}}\>,$ (28)

where

$\begin{dcases}H_{\mathrm{H}}=-\dfrac{1}{2}\sum_{\langle
i,j\rangle}J_{ij}(r_{ij})\mathbf{S}_{i}\cdot\mathbf{S}_{j}\>,\\\
H_{\mathrm{L}}=\sum_{i}\left(AS_{i}^{2}+BS_{i}^{4}+CS_{i}^{6}\right)\>.\end{dcases}$
(29)

Here, $H_{\mathrm{H}}$ describes an Heisenberg model, where the exchange
parameter between the pair of atoms $\langle i,j\rangle$, baring spin
$\mathbf{S}_{i}$ and $\mathbf{S}_{j}$, depends on the atoms distance $r_{ij}$.
Note that the spin vectors are in units of $\hbar$, so that
$S_{i}=M_{i}/g_{e}\mu_{\mathrm{B}}$, with $M_{i}$ being the $i$-th local
magnetic moment and $\mu_{\mathrm{B}}$ the Bohr magneton. In particular, in
this work we choose the following functional form [24] for $J_{ij}$

$J_{ij}(r_{ij})=J_{n}(1-\Delta r_{ij}/r_{n})^{3},\qquad\text{with}\qquad\Delta
r_{ij}=r_{ij}-r_{n}\>,$ (30)

where the index $n$ indicates that the atoms $i$ and $j$ form a $n$-th
neighbours pair. The distance $r_{n}$ is that between two $n$-th neighbour
atoms in the undistorted bcc lattice (the $n$-th neighbour equilibrium
distance). Similarly, the constants $J_{n}$ are the Heisenberg coupling
elements between two $n$-th neighbours at equilibrium. It should be noted that
$H_{\mathrm{H}}$ describes coupling between the position and the spin degrees
of freedom by mean of the coupling constants, $J_{ij}(r_{ij})$. The
Hamiltonian is then completed by a Landau-like term, $H_{\mathrm{L}}$, which
describes the dependence of the energy on the longitudinal local magnetization
(the magnitude of local spins) [8], where $A,B$ and $C$ are constants to be
determined.

In this work we set the various parameters to describe bcc iron [8]. Thus, the
Heisenberg exchange interaction extends to second nearest neighbours with
$J_{1}=22.52$ meV and $J_{2}=17.99$ meV, while the Landau parameters are
chosen to be $A=-440.987$ meV, $B=150.546$ meV and $C=50.769$ meV. We will now
proceed to show how our descriptors are able to capture the potential energy
surface of the Hamiltonian of Eq. (28).

##  III Numerical Simulations

The Hamiltonian given in Eq. (28) consists in two qualitatively different
terms, $H_{\mathrm{H}}$ and $H_{\mathrm{L}}$. The first describes transverse
energy excitations and spin-lattice coupling, while the second accounts for
longitudinal excitations. In order to investigate the descriptive power of the
power spectrum and of our linear energy model over these two different types
of excitation, we first consider only the Heisenberg term with fixed magnetic
momenta lengths, $M_{i}=2.2$ $\mu_{\textrm{B}}$, and later the full model.

###  III.1 Representing the Heisenberg Model with spin-lattice interaction

The dataset has been built by displacing the atomic positions from the ideal
equilibrium bcc structure and by choosing different orientations of the local
magnetic moments. For the atomic positions we have chosen three sets
presenting a different maximum displacement of 5%, 10% and 20% of the lattice
constant, respectively. The sampling of the displacements is uniform in space.
In contrast, we have used two different strategies to define the spin
structure. In the first one, we align the majority of the spins along the
$z$-axis, while the remaining magnetic moments point in a random direction.
More specifically, out of the available 219 magnetic moments, we randomly
choose always more than 200 spins (the actual number is between 200 and 219,
and it also randomly selected), to be aligned along the $z$-axis, thus forming
an almost-ferromagnetic structure (this training set is called the
‘ferromagnetic’ one). The second strategy, instead, consists in assigning to
all the magnetic moments a random orientation (‘random’ training set).
Considering the three different choices for the maximum atomic displacement,
and the two for the spins alignment, we have thus built a total of $6$
datasets, each made of $100$ configurations.

In order to test the predictions made by our fitted spin potentials, we build
three further test-sets for each of the $6$ dataset explored. The first set
consists of $219$ configurations; the $n$-th configuration having $n$
randomly-chosen spins aligned along the $z$-axis, while the remaining ones
being randomly oriented in space. In this case the lattice is chosen to be
pristine bcc (no atomic displacements) so to test independently the vectorial
character of the potential. In contrast, the second and the third test-sets
are designed to investigate also the atomic displacements. They consist of
$50$ configurations each, and the atomic displacement has the same maximum
magnitude of that of the dataset used to train the model. The spins
configuration of the second test-set has 200 randomly-chosen spins aligned
along the $z$-axis while the remaining ones are randomly oriented in space.
The third test-set, instead, has all the spins randomly oriented in space. As
such, these sets have been designed to test the predictions on a mostly-
ferromagnetic environment and on a paramagnetic one, respectively.

Figure 3:  (Color online) Predicted against actual energies for a Ridge linear
regression trained on the ferromagnetic dataset with $10$% maximum atomic
displacement. The actual energies are from the Heisenberg model with spin-
lattice coupling. The results are tested against three different test-sets.
The red dots represent the configurations with undistorted bcc atomic
positions and progressively $z$-aligned magnetic moments. The other dots
represent different displacements of the atoms for near to ferromagnetic and
paramagnetic configurations, respectively. The figure demonstrates the good
agreement of the potential also for configurations, which are energetically
far from that used in the training (blue region). Zoom-in around different
energy regions are displayed in the inserts.

In Fig. 3 we shown the results for the ferromagnetic-trained potential with a
10% maximum atomic displacement, for which we explicitly report the procedure
and the results. The results on the other training sets are reported in the
Supplementary Informations. The optimal potential parameters are found to be
$n_{\text{max}}=4$, corresponding to 35 features only, and
$r_{\text{cut}}=1.4$ (lattice units), while the regularization constant of the
Ridge-regression is $\alpha=$3.2\text{\times}{10}^{3}$$. In the cross-
validation procedure, we split the dataset in training- and test-sets five
times with a 80-20 ratio with respect to the total dataset. We obtain an
energy mean absolute error (MAE) of $(4.83\pm 0.15)\times 10^{-5}$ eV/atom on
the training set, and of $(6.8\pm 0.7)\times 10^{-5}$ eV/atom on the test one,
which roughly corresponds to an error smaller than 0.1%. When looking more
specifically at the model predictions, the analysis on the first test-set (red
points in Fig. 3) returns us a MAE of $5.6\text{\times}{10}^{-4}$ eV/atom
($\sim 1$%). Notably we find that the MAE of configurations with an energy
above $-$0.01$$ eV/atom, namely those that are further away from the energy
range of the training set, is $8.2\text{\times}{10}^{-4}$ eV/atom. This means
that the prediction of the model is still effective also in the portion of the
configuration space far from that of the training. The MAEs on the second
(green points in Fig 3) and third test-sets (blue points in the Figure) are
respectively $1.9\text{\times}{10}^{-4}$ eV/atom and
$6.0\text{\times}{10}^{-4}$ eV/atom. This suggests that the model can
extrapolate rather well across the configuration space.

The remaining trained models show us MAEs similar to that reported above, when
the training is performed over the ferromagnetic datasets. In contrast, when
training on the three training sets denoted as “random” spin configurations,
we notice a significantly larger MAEs for large atomic displacements.
Explicitly, the MAE is $1.3\text{\times}{10}^{-4}$ eV/atom for the random
training set with 5% maximum atomic displacement, but this already increases
to $2.5\text{\times}{10}^{-3}$ eV/atom for 10% maximum displacement and
reaches up to $0.01$ eV/atom for 20% maximum displacement. The failure for the
largest maximum displacement can be attributed to the fact that the random
dataset explores a much smaller portion of the energy landscape of the model.
This is because the random configurations are all characterized by a small
total magnetization and hence rather similar energies. For all the other
cases, the good agreement obtained between our model and the true potential
energy surface, demonstrates that the Heisenberg model including spin-lattice
coupling is accurately described by our potential, which is able to
extrapolate the entire energy landscape. It is crucial to remark at this point
that we did not introduce any prior knowledge of the functional dependence of
the coupling constant on the pair-wise distance between the atoms, namely the
model and the descriptors are able to autonomously interpolate the spin-phonon
coupling. Moreover, given the modest size of the dataset, we found quite
remarkable that simultaneous using a small datasets and a reduced number of
features are able to reach the accuracies reported here.

Next we will consider the full Hamiltonian, including $H_{\mathrm{L}}$,
describing both transverse and longitudinal spin excitations.

###  III.2 Heisenberg Model with longitudinal excitations

The investigation of the model described by the complete Hamiltonian of Eq.
(28) follows the same approach used for the analysis on the Heisenberg part.
In this case we build a dataset corresponding only to one maximum displacement
of the atomic positions, namely 10% of the lattice constant. The spin
configurations correspond to the “ferromagnetic” case described in the
previous section. However, having to deal with longitudinal excitations as
well, we also vary the magnetic moment’s length. In particular, the magnetic
moments aligned along the $z$-axis are chosen to be 2.25 $\mu_{\textbf{B}}$,
while the randomly oriented ones have a length randomly chosen in the range
1.8-2.3 $\mu_{\textbf{B}}$.

Figure 4:  (Color online) Predicted against actual energies for both a Ridge
linear regression (red circles) and a GAP (blue triangles) constructed to
predict the potential energy surface of the full Hamiltonian of Eq. (28)
(Heisenberg model with spin-lattice coupling and longitudinal spin
fluctuations). Results are presented for the test set. While the accuracy of
the two models is similar for energies in the training set region, this
becomes significantly different for the high-energy data points. In
particular, the linear model accurately predicts energies up to about -0.26
eV/atom, but then significantly deviates from the parity line (dashed black
line). In contrast, the GAP, while it appears to slightly underestimate the
actual energies, outperforms the linear model at extrapolating away from the
training set region.

When testing the predictions, we build an additional test-set, corresponding
to the first one presented in the previous section, namely containing an
increasing number of aligning spins. In this case, the length of the
$z$-aligned magnetic moments are again fixed to $2.25$ $\mu_{\textbf{B}}$,
while the randomly oriented ones have a length in the range $1.9-2.3$
$\mu_{\textbf{B}}$. Also the cross-validation procedure is similar to the one
employed before with a five-time split of the dataset into training and test
sets, with a 80-20 ratio. The parameters chosen are then $n_{\text{max}}=4$,
$r_{\text{cut}}=1.4$ (lattice constant), $\alpha=$2\text{\times}{10}^{5}$$.

The MAEs obtained in this case are $(4.9\pm 0.1)\times 10^{-4}$ eV/atom and
$(6.0\pm 0.5)\times 10^{-4}$ eV/atom, respectively for the training and test
set. These values are about one order of magnitude larger than those obtained
previously for the Heisenberg model with spin-lattice coupling. We can
understand such accuracy loss by noticing that the descriptors are quadratic
in the spin magnitude, as evident from Eqs. (16), (23) and (25). Therefore, a
linear machine-learning model, as employed here, will not be able to capture
the energy contributions to the fourth and sixth power in the magnetization,
which defines the longitudinal part of the Hamiltonian, $H_{\textrm{L}}$. In
fact, it may be surprising that the model still performs accurately even in
this case. This is because we are exploring a region of the potential energy
surface relatively close to the minimum, where the energy contributions in
$S^{4}$ and $S^{6}$ remain modest.

In order to corroborate this hypothesis, we evaluate the model predictions on
a test set containing progressively-aligned spins, for which we obtain a MAE
of $6.0\text{\times}{10}^{-3}$ eV/atom. These results are shown as red dots in
Fig. 4, where it is clear that the Ridge regression performs poorly as we
progressively explore energy regions away from the training range. Such
behaviour must be associated to the limit of the machine-learning linear model
constructed over our description. In more detail, we find that the low-energy
regions are still well described, with a MAE of $7.1\text{\times}{10}^{-4}$
eV/atom for energies smaller than less than -0.32 eV/atom (in the training
range). In contrast, the potential rapidly departs from the parity line at
higher energies, where we compute a MAE of $1.0\text{\times}{10}^{-2}$ eV/atom
for data above -0.26 eV/atom.

We can improve on the error and go beyond the quadratic nature of our
descriptors by combining the power-spectrum representation of the atomic and
vector field with a non-linear machine-learning model. In particular, we
consider here a Gaussian approximation potential (GAP) [25]. GAP expresses the
atomic energy of the $i$-th atom as

$\varepsilon_{i}=\sum_{t=0}^{N_{\mathrm{train}}}\theta_{t}S(\bm{p}^{i},\bm{p}^{t})=\sum_{t}\theta_{t}e^{-{\frac{1}{2\sigma}(\bm{p}^{i}-\bm{p}^{t}})^{2}}$
(31)

where the sum is extended over all the atoms in the training set, and where
$\bm{p}^{i}$ and $\bm{p}^{t}$ are the power spectrum, respectively of the
$i$-th atom and of the training set. The non linearity of the similarity
kernel, $S$, allows us to describe energy contributions going beyond the
quadratic order in the spin magnitude.

The GAP predictions obtained over the test-set are shown in Fig. 4 as blue
triangles. We notice that the MAE associated with the configurations having
energy smaller than -0.32 eV/atom remains very close to that obtained with
ridge regression. However, the total MAE decreases to
$4.5\text{\times}{10}^{-4}$ eV/atom and, most importantly, the MAE for
energies larger than -0.26 eV/atom is now reduced to
$5.8\text{\times}{10}^{-3}$ eV/atom, namely is halved. In fact, the figure
clearly shows that the non-linear GAP improves the ability of the model to
extrapolate away from the training set range. Since, the actual potential
energy surface for a spin system, as the one obtainable from density
functional theory, is expected to include energy contributions going beyond a
quadratic dependence on the magnetization, we conclude that the best use of
our representation will be in conjunction with non-linear machine-learning
models.

##  IV Conclusion

In this work we have introduced a new invariant power spectrum representation
for vectorial fields. After having presented an in-depth analysis of its
rotational invariance and basic properties, we have designed a linear energy
model, closely following the SNAP [22] approach. Such spin SNAP has then be
put to the test against a bcc iron model described by an Hamiltonian
containing spin-lattice coupling and both transverse and longitudinal energy
excitations. Spin-lattice coupling is introduced by mean of a Heisenberg model
with exchange parameters depending on the interatomic distance, while the
longitudinal spin excitations are described by a simple Landau term containing
even powers of the magnetization.

We have then trained a first potential, linear in the power spectrum, for the
situation where the longitudinal spin excitations are neglected. This was
trained over a dataset obtained by displacing the atomic positions and the
orientations of the atomic-magnetic moments, comprising a total of only 100
different configurations. Our results showed that the power spectrum is able
to describe the entire energy surface, by accurately extrapolating far beyond
the energy range covered by the training set. This proves that a linear model
using the power spectrum is sufficient to describe both the Heisenberg model
and the spin-lattice coupling, already from a small dataset. Crucially, no
prior information on the dependence of the exchange constants on the atomic
position were used by the model.

We have then repeated the exercise for the complete Hamiltonian, containing
also the Landau term, by training over a dataset containing spins of different
magnitude. Our results are highly accurate for configurations with energies
within the range explored by our training set, but the model does not perform
well in extrapolating. We have attributed this result to the inability of the
power spectrum, combined with a linear machine-learning model, to describe
energy contributions scaling beyond a quadratic dependence on the spin
magnitude. Such shortcoming can be recovered by employing a non-linear model.
Thus, we have investigated a Gaussian approximation potential and shown that
extrapolation over a much-larger portion of the energy landscape is indeed
possible.

All in all, our analysis has shown that a power spectrum representation of the
magnetization field can be used, together with non-linear machine-learning
models, as an efficient descriptor of spin potential energy surfaces. This can
now be used in conjunction with training sets obtained from accurate
electronic structure theory to predict finite temperature properties of
magnets.

###### Acknowledgements.

This work has been supported by the Irish Research Council Advanced Laureate
Award (IRCLA/2019/127), and by the Irish Research Council postgraduate program
(MC). We acknowledge the DJEI/DES/SFI/HEA Irish Centre for High-End Computing
(ICHEC) and Trinity Centre for High Performance Computing (TCHPC) for the
provision of computational resources.

##  Appendix A Atom-centered power spectrum

Let us prove that, if the vector field has an atom at the origin, the power
spectrum of (16) will not be generally rotationally invariant. For simplicity
we can consider first the trivial case in which our system is made of just a
single atom, e.g. there are no other atoms inside the cut-off radius. The
component $u_{n100}$ of the coefficients of Eq. (23) reads

$\displaystyle u_{n100}$ $\displaystyle=$ $\displaystyle
w_{0}R_{n1}(0)\sum_{mq}Y^{m*}_{1}(\bm{0})C^{00}_{1m1q}v_{q}\>,$ (32)
$\displaystyle\propto$ $\displaystyle R_{n1}(0)\sum_{q}C^{00}_{101q}v_{q}\>,$
$\displaystyle\propto$ $\displaystyle R_{n1}(0)v_{0}\>,$

where in the first step we have used the fact that the spherical harmonics
along the $z$-axis vanish unless $m=0$, and in the second step we have
considered the equality $q=M-m$ implemented through the Clebsch-Gordan
coefficients (for simplicity we do not carry over unessential constant). The
power spectrum component for this case then is simply

$p_{n10}\propto\absolutevalue{R_{n1}(0)v_{0}}^{2}\>.$ (33)

Equation (33) establishes that the $n10$ component of the power spectrum is
proportional to the magnitude of the $z$ component of the vector field at the
origin. Crucially, for a general rotation of the reference frame, the $z$
component becomes mixed with the other components, so that it changes its
value. We then deduce that the power spectrum is not rotational invariant, if
the vector field does not vanish at the origin of the local reference frame.
It is worth stressing that this proof holds, since the origin is a fixed point
for the rotation. If the center or rotation is not the origin, then also the
argument of the spherical harmonics will rotate, making the $m\neq 0$ terms
relevant too. This proof can be generalized also to the case in which there
are other atoms within the cutoff radius. As pointed out in the main text, we
can however recover a cylindrical symmetry by simply rotate the system so that
the vector fields in the origin points in the same direction as the $z$-axis.
Another possible solution stems from the presence of the radial-basis set in
Eq. (33) which, if carefully chosen, could remove the symmetry breaking terms
of the power spectrum. Spherical Bessel functions are suitable for this
purpose, as we will show now. We first notice that for $l\neq 0$ the function
vanishes at the origin, namely $g_{n-l,l}(0)=0$. This is a consequence of the
fact that the basis set is defined in terms of the spherical-Bessel functions
$j_{l}(x)$, which also vanish at the origin for non-zero $l$. Thus, for $l\neq
0$ all the contributions arising from an atom at the origin are removed,
ensuring the rotational invariance. We are then left to prove that the power
spectrum is invariant also for $l=0$. This is easily done by noticing that the
spherical harmonic $Y^{0}_{0}$ is just a constant, with no angular dependence,
completing the demonstration.

##  Appendix B Power spectrum connecting different radial channels

Following the same approach used in reference [12], we can generalize the
expression for the power spectrum so that different radial channels are
coupled. A generalized expression for the power spectrum may read

$p_{nn^{\prime}lJ}=\sum_{M}u_{nlJM}^{*}u_{n^{\prime}lJM}\>.$ (34)

The rotational invariance is still ensured by the fact that the transformation
rules for the expansion coefficients $u_{nlJM}$ involve only the Wigner
D-matrix belonging to the angular momentum $J$ [see Eq. (18)]. However, it is
apparent that the number of components defining the descriptors are increased
with this coupling choice. Note that we must ensure the same $l$ for the two
factors in the sum above, so that the resulting quantity is real: the only
complex part are found in the angular dependent terms being the radial
functions real.

##  Appendix C Generalization to a tensorial density

In this appendix we generalize our formalism from a vectorial field to a
tensorial one. In this case, the local density reads

$\bm{\rho}(\mathbf{r})=\sum_{a}w_{a}h_{a}(\mathbf{r}-\mathbf{r}_{a})\bm{T}^{n_{a}}_{a}\>,$
(35)

where $\bm{T}^{n_{a}}_{a}$ is a $n_{a}$-rank tensor associated to the $a$-th
atom. Note that tensors associated with different atoms can have different
ranks as, for example, when dealing simultaneously with a scalar field and a
vector field. We can then use the spherical decomposition of the tensors and
write $\bm{T}^{n_{a}}_{a}$ in Dirac notation as

$\ket{\bm{T}^{n_{a}}_{a}}=\sum_{\lambda}^{n_{a}}\sum_{\mu=-\lambda}^{\lambda}(T_{a})^{\mu}_{\lambda}\ket{\lambda\mu}\>,$
(36)

where $(T_{a})^{\mu}_{\lambda}$ is the spherical component of the tensor
$\bm{T}^{n_{a}}_{a}$ relative to the spherical basis $\ket{\lambda\mu}$. This
is analogous to $v_{a,q}$ and $\ket{q}$ for the case of a vectorial field. The
components $(T_{a})^{\mu}_{\lambda}$ transform as the the spherical harmonic
$Y^{\mu}_{\lambda}$ under rotation. We can then write the density as

$\ket{\bm{\rho}}=\sum_{n=0}^{n_{\text{max}}}\sum_{l=0}^{n}\sum_{m=-l}^{l}\sum_{\lambda=0}^{\Lambda}\sum_{\mu=-\lambda}^{\lambda}c_{nlm\lambda\mu}\ket{nlm\lambda\mu}\>,$
(37)

where the expansion has been truncated at $n_{\text{max}}$, and with
$\Lambda=\max_{a}(n_{a})$ being the highest rank of the tensors in the
tensorial field additional zero coefficients can be introduced to have an
homogeneous representation in the highest-tensorial rank. The expansion
coefficients are then obtained by projection as

$c_{nlm\lambda\mu}=\sum_{a}w_{a}(T_{a})^{\mu}_{\lambda}\int\differential
r\differential\hat{\bm{r}}\>r^{2}R_{nl}(r)Y^{m*}_{l}(\hat{\bm{r}})h_{a}(\bm{r}-\bm{r}_{a})\>,$
(38)

with $(T_{a})^{\mu}_{\lambda}=0$ if $\lambda>n_{a}$. Following the same
procedure outlined previously for the case of a vectorial field, we can
express the density over the coupled basis as

$\ket{\bm{\rho}}=\sum_{nl\lambda JM}u_{nl\lambda JM}\ket{nl\lambda JM}\>,$
(39)

by mean of the coupling scheme,

$\ket{nl\lambda
JM}=\sum_{m=l}^{l}\sum_{\mu=-\lambda}^{\lambda}C^{JM}_{lm\lambda\mu}\ket{nlm\lambda\mu}\>,$
(40)

with $\absolutevalue{l-\lambda}\leq J\leq l+\lambda$. The coupled-basis
coefficients are given in terms of the uncoupled ones as

$u_{nl\lambda JM}=\sum_{m\mu}C^{JM}_{lm\lambda\mu}c_{nlm\lambda\mu}\>.$ (41)

Finally, the power spectrum is again given by squaring the coefficients

$p_{nl\lambda J}=\sum_{M}\absolutevalue{u_{nl\lambda JM}}^{2},$ (42)

and it is invariant under simultaneous rotations of the frame of reference and
the tensorial field. A further generalization to multi-channel coupling can be
obtained by using the same argument presented in Appendix B.

## References

  * [1] P. Elliott, T. Müller, J.K. Dewhurst S. Sharma and E.K.U. Gross, Ultrafast laser induced local magnetization dynamics in Heusler compounds, Sci. Rep. 6, 38911 (2016).
  * [2] J. Simoni, M. Stamenova and S. Sanvito, Ab initio dynamical exchange interactions in frustrated antiferromagnets, Phys. Rev. B 96, 054411 (2017).
  * [3] J. Simoni, M. Stamenova and S. Sanvito, Ultrafast demagnetizing fields from first principles, Phys. Rev. B 95, 024412 (2017).
  * [4] D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press (2002).
  * [5] R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis, and R. W. Chantrell, Atomistic spin model simulations of magnetic nanomaterials, J. Phys.: Condens. Matter 26, 103202 (2014).
  * [6] O. Eriksson, A. Bergman, L. Bergqvist and J. Hellsvik, Atomistic Spin Dynamics: Foundations and Applications, Oxford University Press, (2017).
  * [7] P.-W. Ma and S.L. Dudarev, Atomistic Spin-Lattice Dynamics, in W. Andreoni and S. Yip S. (eds) Handbook of Materials Modeling, Springer (2018).
  * [8] P.-W. Ma and S.L. Dudarev, Longitudinal magnetic fluctuations in Langevin spin dynamics, Phys. Rev. B 86, 054416 (2012).
  * [9] M.O.A. Ellis, M. Stamenova and S. Sanvito, Multiscale modeling of current-induced switching in magnetic tunnel junctions using ab initio spin-transfer torques, Phys. Rev. B 96, 224410 (2017).
  * [10] F. Musil, A. Grisafi, A.P. Bartók, C. Ortner, G. Csányi and M. Ceriotti, Physics-Inspired Structural Representations for Molecules and Materials, Chem. Rev. 121, 9759 (2021).
  * [11] J. Behler and M. Parrinello, Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces, Phys. Rev. Lett. 98, 146401 (2007).
  * [12] A.P. Bartók, R. Kondor and G. Csányi, On representing chemical environments, Phys. Rev. B 87, 184115 (2013).
  * [13] S. Nikolov, M.A. Wood, A. Cangi, J. Maillet, M. Marinica, A.P. Thompson, M. Desjarlais and J. Tranchida, Data-driven magneto-elastic predictions with scalable classical spin-lattice dynamics, npj Comput. Mater. 7, 153 (2021).
  * [14] M. Eckhoff and J. Behler, High-dimensional neural network potentials for magnetic systems using spin-dependent atom-centered symmetry functions, npj Comput. Mater. 7, 170 (2021).
  * [15] R. Drautz, Atomic cluster expansion of scalar, vectorial, and tensorial properties including magnetism and charge transfer, Phys. Rev. B 102, 024104 (2020)
  * [16] M. Weissbluth, Atoms and Molecules, Academic Press (1978).
  * [17] D.A. Varshalovich, A.N. Moskalev and V.K. Khersonskii, Quantum Theory of Angular Momentum, World Scientific (1988).
  * [18] E. Kocer, J.K. Mason and H. Erturk, Continuous and optimally complete description of chemical environments using Spherical Bessel descriptors, AIP Adv. 10, 015021 (2020).
  * [19] A. Grisafi, D.M. Wilkins, G. Csányi and M. Ceriotti, Symmetry-Adapted Machine Learning for Tensorial Properties of Atomistic Systems, Phys. Rev. Lett. 120, 036002 (2018).
  * [20] A. Glielmo, P. Sollich and A. De Vita, Accurate interatomic force fields via machine learning with covariant kernels, Phys. Rev. B 95, 214302 (2017).
  * [21] V.H.A. Nguyen and A. Lunghi, Predicting tensorial molecular properties with equivariant machine-learning models, Phys. Rev. B 105, 165131 (2022).
  * [22] A.P. Thompson, L.P. Swiler, C.R. Trott, S.M. Foiles and G.J. Tucker, Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials, J. Comp. Phys. 285, 316 (2015).
  * [23] M.J. Willatt, F. Musil and M. Ceriotti, Atom-density representations for machine learning, J. Chem. Phys. 150, 154110 (2019).
  * [24] P.-W. Ma, C.H. Woo and S.L. Dudarev, Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron, Phys. Rev. B 78, 024434 (2008).
  * [25] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons, Phys. Rev. Lett. 104, 136403 (2010).


# Small UAVs-supported Autonomous Generation of Fine-grained 3D Indoor Radio
Environmental Maps

Ken Mendes, Filip Lemic, Jeroen Famaey Email<EMAIL_ADDRESS>Internet and Data Lab (IDLab), Universiteit Antwerpen - imec, Belgium

###### Abstract

Radio Environmental Maps are a powerful tool for enhancing the performance of
various communication and networked agents. However, generating REMs is a
laborious undertaking, especially in complex 3-Dimensional (3D) environments,
such as indoors. To address this issue, we propose a system for autonomous
generation of fine-grained REMs of indoor 3D spaces. In the system, multiple
small indoor Unmanned Aerial Vehicles are sequentially used for 3D sampling of
signal quality indicators. The collected readings are streamlined to a Machine
Learning (ML) system for its training and, once trained, the system is able to
predict the signal quality at unknown 3D locations. The system enables
automated and autonomous REM generation, and can be straightforwardly deployed
in new environments. In addition, the system supports REM sampling without
self-interference and is technology-agnostic, as long as the REM-sampling
receivers features suitable sizes and weights to be carried by the UAVs. In
the demonstration, we instantiate the system design using two UAVs and show
its capability of visiting 72 waypoints and gathering thousands of Wi-Fi data
samples. Our results also include an instantiation of the ML system for
predicting the Received Signal Strength (RSS) of known Wi-Fi Access Points at
locations not visited by the UAVs.

##  I Introduction

A Radio Environmental Map (REM) documents radio signal properties over a given
geographic area. These properties can among others include the frequency,
protocol and technology of the radio signal, as well as an indication of the
quality of the signal (e.g., Received Signal Strength (RSS), Signal-to-Noise
Ratio (SNR), or Channel State Information (CSI)), and are stored together with
the location where they were measured. These REMs and the data they hold can
then be used for a variety of purposes, for example as an aid in cognitive
radio networks [1], for RF-based localization [2], or for optimizing network
discovery and handover procedures [3].

As a trade-off to their utility, the generation of REMs is a burdensome
process, particularly for 3D environments [4]. Apart from the significant
labor needed for carrying out a measurement campaign, an additional complexity
stems from the fact that the measurements have to be correlated with the
physical locations at which they are collected, suggesting the need for
accurate localization of the entity collecting such measurements. The
complexity of REM generation is further exacerbated by the fact that such
generation often has to be performed periodically, as the REMs can become
obsolete due to long-term changes in the signal propagation [5]. Finally, as a
variety of different wireless technologies often co-exist [6], the generation
of the corresponding REMs is ideally to be done with a single tool for the
coexisting technologies.

From the above discussion, there is a need for a system for autonomous
generation of 3D REMs, which has been recognized in the community. One of the
most promising approaches for such generation is to utilize UAVs as carriers
of the REM-sampling devices [7, 8, 9]. In such systems, the samples generated
by the REM-sampling devices are correlated with the UAV-originating location
at which the measurements are taken, providing the primitive for REM
generation. Such systems have mostly been proposed for outdoor scenarios and
Global Positioning System (GPS)-enabled large outdoor UAVs. These systems are
intuitively not suitable for fine-grained generation of REMs indoors because
the UAVs are not practically utilizable in smaller and more complex
environments. In addition, the GPS-originating location information is not
suitable for the generation of fine-grained REMs due to the relatively large
localization errors that such location information unavoidably features [10].

To address this issue, we propose a system for generating fine-grained 3D REMs
in small indoor environments. The system can be viewed as a toolchain
consisting of small-indoor UAVs with fine-grained localization capabilities
supported by an Ultra Wide-Band (UWB)-based positioning system. The UAVs are
envisioned to serve as the carriers of the REM-sampling receivers and
provisioners of location annotation for the sampled measurements. The
location-annotated measurements are envisioned to be streamlined as training
inputs to an ML-based part of the toolchain, which, once trained, is able to
predict the signal quality at locations not visited by the UAVs.

The most similar existing effort to our work is [11], where the authors
propose a UAV-based system for the generation of a training dataset for 3D Wi-
Fi fingerprinting-based indoor localization. In contrast to their work
focusing solely on Wi-Fi-based fingerprinting, we argue that REMs can be
beneficial and utilized more broadly, for example in optimizing the
positioning of UAVs serving as mobile relays [12] or planning the extensions
of any wireless networking infrastructure by adding Access Points (APs) or
base stations to cover “dark” connectivity regions in an environment of
interest [13]. Hence, our system abides to two additional design requirements
compared to [11]. The first one is a modular design of the interface between
the UAV-based system and a given REM-sampling device. This requirement allows
for a simple integration of different REM-sampling device (e.g., Wi-Fi, LoRa,
BLE, mmWave) with the UAV, extending the REM capabilities beyond the
traditional Wi-Fi. The second is to guarantee no self-interference between the
wireless network used to control the UAVs and the REM-sampling device the UAV
is carrying, allowing highly repeatable measurement collections with minimized
external influences. Finally, in our system additional UAVs can be seamlessly
integrated into the toolchain, allowing for sequential data collection and
scalable REM generation.

We instantiate the proposed system design on two CrazyFlie UAVs to
autonomously gather IEEE 802.11b/g/n Wi-Fi beacon frame data in the 2.4 GHz
Industrial, Scientific and Medical (ISM) band. In a real-world indoor
environment (i.e., a living room in a condo apartment), we demonstrate that
the UAVs are each able to visit 36 waypoints and collectively gather several
thousands of Wi-Fi beacon data samples. Finally, we instantiate the remaining
part of the toolchain on several contemporary ML algorithms, train the system
using the data streamlined by the UAVs, and demonstrate a reasonable
prediction accuracy at locations not visited by the UAVs.

##  II System Overview

The REM generation is envisioned to be initiated from the control station by
providing a set of waypoints to be visited by a fleet of UAVs operating
sequentially. The first UAV in the fleet is envisioned to visit a subset of
the provided points, with the main limitation on the number of points that can
be visited in one run stemming from the constrained battery capabilities of
the UAV. At each location, the UAV instructs the REM-generating receiver to
collect the signal quality indicators of interest, and upon each response it
reports the obtained results to the control station. This procedure is
repeated until all points have been visited by the UAVs or all UAVs in the
fleet have depleted their batteries. In order to visit the instructed
locations, each UAV requires a means for localizing itself, which is supported
through its localization system consisting of a client mounted on the UAV and
a set of infrastructural devices (i.e., anchors) for the client’s
localization.

The main design requirements for the envisioned system include i) accurate
location-annotated sampling for streamlined generation of fine-grained 3D
REMs, ii) straightforward deployment of the system in unknown complex indoor
environments, iii) support for technology-agnostic REM-generating receiver,
iv) guaranteed mitigation of self-interference.

We have compared a range of commercial off-the-shelf UAVs based on the design
requirements. This comparison can be found in the accompanying technical
report [14]. We have decided to utilize BitCraze Crazyflie 2.1 (Figure 4),
primarily as it is an open hardware and software platform. Crazyflie 2.1 UAVs
come with a FreeRTOS-based operating system and a radio and Bluetooth LE
transceivers for control. They also feature an accelerometer, gyroscope,
magnetometer, and a high precision pressure sensor through its 10-DOF Inertial
Measurement Unit (IMU). These capabilities can be extended by adding up to two
expansion boards or decks as shown in Figure 4. In this work, both expansion
slots are used: one for the Loco Positioning Deck (LPD), the other for the
integration of an REM-generating receiver. The Crazyflie provides a set of
interfaces that can be addressed over 20 pins for communication with each
expansion deck. Figure 4 documents the pin allocations, where one can choose
between an I2C, SPI, STM32, and two UART interfaces.

###  II-A Interfacing with REM-generating Receivers

Using the FreeRTOS-flavored Crazyflie 2021.06 firmware release as a basis, a
custom driver is responsible for interfacing with an REM-generating receiver.
The driver should support: i) initializing and ii) checking the state of the
REM-generating receiver, iii) instructing the REM-generating receiver to
collect a measurement, and iv) enabling parsing of the output of the previous
instruction. For integration with the UAV, the user is required to provide the
driver for the REM-generating receiver to react to the four specified
instructions. In terms of hardware integration, the user can choose between
UART and I2C interfaces available on the Crazyflie UAV. We argue that such an
integration procedure is straightforward as it is supported by well-known
hardware interfaces and a four instructions-long C-flavored driver, as long as
the REM-generating receiver features suitable size (i.e., USB-dongle
dimensions) and weigh (up to 20 grams) to be carried by the UAV.

Figure 1: Customized Crazyflie 2.1 UAVs

Figure 2: Crazyflie expansion boards

Figure 3: Crazyflie interfaces and pin allocation

Figure 4: The Crazyflie commander framework

###  II-B UAV Localization

We use the Crazyflie’s Loco Positioning System (LPS) to provide accurate
indoor positioning. This system works with a tag, the LPD, and multiple
anchors distributed over the volume used for localizing the tag. The system
can be deployed by simply positioning of the localization anchors, measuring
their coordinates relative to a chosen origin, and initializing their
automated calibration for synchronizing their transmission schedules. Once the
localization anchors are self-calibrated, they can be used for localizing the
UAVs and consequently for the generation of a 3D REM can be initiated. Given
that the procedure for the REM-generating system is relying solely on physical
deployment, localization, and initiation of their self-calibration, we argue
that it can be utilized for rapid deployment of our system in complex new
environments.

The LPS is based on the Decawave DWM1000 module and uses UWB radio technology
for communication and localization. The tag can estimate its own position
based on the UWB signals received from the localization anchors. The
localization is then performed using either the Two-Way Ranging (TWR)
procedure or different flavors of the Time Difference of Arrival (TDoA)
procedure, the latter featuring slightly better accuracy and supporting
simultaneous localization of multiple UAVs. Regardless of the utilized
localization procedure, the LPS is able to localize the tag in a 3D
environment at the range of about 10 m [11].

A minimum of four localization anchors are required to support UAV
localization in 3D environments. A Crazyflie equipped with an LPS deck makes
use of an extended Kalman filter to estimate its state (orientation and
position), its implementation is based on [15]. An increase in the number of
localization anchors increases the robustness and accuracy of this process.
Hence, Bitcraze advises to use at least six such anchors to mitigate potential
negative effects in the deployment environment (e.g., no Line of Sight (LoS)
to an anchor, fluctuating tag orientation). Chekuri and Won have shown in
their tests [11] that localization with 6 anchors can achieve an accuracy of 9
cm when the UAV is hovering, this is important as we will record the UAV’s
position and scan while the UAV is holding its position and orientation. In
summary, the system is able to generate location-annotated measurements for
REM generation with decimeter-level accuracy at 10 m range.

###  II-C UAV Communication and Control

Software control of the UAV is done through a Python application that makes
use of the Crazyflie Python library. This application can communicate with the
UAV to send instructions to move to a waypoint, scan for signal quality
indicators and other relevant information, parse the results, and store them
for later processing. It runs through the following sequence: i)
initialization and instruct the UAV to take-off. For every configured
waypoint, make the UAV: ii) move to the waypoint defined with the $\langle
x,y,z\rangle$ tuple, iii) initiate an on-demand scan, iv) shutdown the
Crazyradio while the scan is running, v) restart the radio connection after
the scan is done, vi) fetch the scan results, parse and store them, land the
UAV and shut it down, and store the results. In terms of technology,
communicating with and controlling the Crazyflie UAvs remotely can be done
using a custom USB dongle called the Crazyradio. This radio uses a nRF24LU1
chip providing 126 channels and works with a custom protocol to communicate
with the Crazyflie: the Crazyradio RealTime Protocol (CRTP). These 126
channels are uniformly distributed over the 2400 MHz - 2525 MHz frequencies.

Figure 5: The number of APs detected per IEEE 802.11 channel with the
Crazyradio set at different frequencies or completely turned off

Figure 6: Number of samples per UAV and scanned location

Figure 7: Histograms showing the number of samples collected per bin of 0.5 m
along the x and y-axis

Figure 8: RMSE of prediction for different models

There are three potential sources of self-interference from the system toward
the REM-generating receiver. The first two are the positioning system and the
propulsion of UAV’s rotational engines, although [11] showed that these types
of interferences have a negligible effect on the REM-generating receiver
operating in 2.4 GHz Industrial, Scientific, and Medical (ISM) band. The third
potential source of interference is the Crazyradio. The interference generated
by the Crazyradio while the System Under Test (SUT) is scanning for available
APs at a given location is illustrated in Figure 8. The figure shows the
average number of APs detected at different 2.4 GHz Wi-Fi channels for 6
operating frequencies of the Crazyradio, as well as in case when the radio was
turned off. As visible, the interference from the Crazyradio is significant,
irrespective of its operating frequency.

Hence, our experimentation setup features the (default) possibility of
automatically turning off the Crazyradio while performing a measurement. To
avoid self-interference, the radio is shut down right before the scan starts
and restarted again after the scan has finished. The Crazyflie will go into
position hold mode while the radio is shut down. Small adjustments were made
to the firmware to make working without a radio possible. First, the
`CRTP_TX_QUEUE_SIZE` was increased so that full scan results can be
temporarily stored until the radio comes back online and the results can be
sent to the controlling application. Second, the
`COMMANDER_WDT_TIMEOUT_SHUTDOWN` was increased to 10 sec. This timeout is a
safety measure, if there is no set-point received within this interval, the
Crazyflie will shut down under the assumption that something went wrong. The
default value does not allow to bridge the radio shutdown period.

When the UAV loses its radio connection, it also loses its ability to get new
waypoints (i.e., target locations) from the base station. When no new setpoint
is received for over 500 ms, the UAV will set its attitude angles (pitch, roll
and yaw) to 0 in order to keep itself stabilized. Figure 4 details how the
base station’s custom Python client can forward waypoints to the Commander in
the UAV’s firmware by making use of the CFlib library. In order to make the
UAV hold its position after shutting down the radio connection, an extra
FreeRTOS task was added to the driver of the ESP8266 deck that will feed back
the scanning position every 100 ms to the UAV’s commander during such a scan.
This task gets resumed at the start of the scanning task and suspended at the
end of it so that it does not interfere with regular waypoint activities. This
feedback process results in the UAV having stability and guaranteed lack of
interference while scanning.

##  III Validation

###  III-A Collection of 3D REM-generating Measurements

We demonstrate111Demonstration video: https://youtu.be/fxDkR-Qat6w one
instantiation of the proposed system design for small UAVs-supported
autonomous generation of fine-grained 3D indoor REMs. In particular, the
system was deployed in a living room of an apartment in a large apartment
building in Antwerp, Belgium. The 3D volume for the UAVs to scan is a
rectangular cuboid of 3.74 m long (x-axis), 3.20 m wide (y-axis) and 2.10 m
high (z-axis). At each of the 8 corners of the cuboid, a localization anchor
was placed and manually localized for enabling the UAVs to estimate their
locations within the volume. Once the self-calibration procedure was finished,
using the controlling application we have instructed a fleet of Crazyflie UAVs
to generate an REM of the volume.

We have opted at generating an REM of 2.4 GHz ISM Wi-Fi. This was supported
through AI Thinker ESP-01 modules with Espressif Systems ESP8266 Wi-Fi chips.
The modules were soldered on Crazyflie prototyping decks and enabled to
interact with the UAVs as one of their expansion decks. Using the Crazyflie
2021.06 firmware release as basis, a custom driver for this ESP-01 module was
written. This driver communicates with the ESP-01 module over its UART
interface by sending AT instructions and parsing the output. Since the module
is only used to scan for available access points, it suffices that the driver
supports just the following AT instructions: i) `AT` \- testing AT start-up,
ii) `AT+CW_MODE_CUR` \- setting the current Wi-Fi mode (to put the module in
station mode), iii) `AT+CWLAP` \- listing the available APs and scanning for
Wi-Fi beacons, `AT+CWLAPOPT` \- formatting the output of the `AT+CWLAP`
instruction to $\langle ssid,rssi,mac,channel\rangle$ tuples.

Wi-Fi REM was selected due to self-interference with CrazyRadio, as its
frequency range overlaps with the 2.4 GHz band that the Wi-Fi modules use. In
order to get an idea of how pronounced this expected interference is, the
Crazyradio was run on different frequencies over its range in 25 MHz
increments (2400, 2425, 2450, 2475, 2500 and 2525 MHz). At each of these
frequencies, 3 access point scans were done using the ESP-01 module on the
Crazyflie. To have a baseline for comparison, an additional three scans were
done with the Crazyradio turned off. The interference generated by the
Crazyradio while scanning for access points is illustrated in Figure 8. It
shows for every Wi-Fi channel the average count of access points that were
detected over the 3 runs and shows this for 6 different frequencies of the
Crazyradio as well as the radio turned off. Wi-Fi channels that did not have
any detected APs were left out for clarity. These scans were done in a short
timespan and with the Crazyflie and Crazyradio in a fixed position. Not only
does this data clearly show that the interference of the Crazyradio is
significant, it also demonstrates the benefits of our design decision to turn
off the Crazyradio while performing REM-generating measurements.

To mitigate interference among UAVs, the UAVs are run in a sequence, not
jointly. The LPS is configured to use the Time Difference of Arrival-based
localization procedure. The Crazyflie is advertised as having a flight time of
up to 7 min depending on how it is used. This is, however, without the weight
and power consumed by the LPD and the custom ESP8266 deck. There are also
several other factors that can influence the UAV’s endurance including flight
and scan parameters, the choice between TWR and TDoA or the distance to the
anchors, to name a few. To get a notion of the UAV’s endurance in a baseline
scenario, a UAV was manually flown until it became less responsive and its
motions erratic, considering a fully charged standard battery, eight active
anchors in TWR mode, periodic scanning mode with an interval of 8 sec, with a
beacon scan duration of around 2 sec. The UAV was kept in a steady position
about 1 m above ground level for the duration of the test. The UAV was able to
perform 36 scans over a timespan of 6 min and 12 sec before it experienced
erratic behaviour.

Obviously, the endurance is expected to be lower when the UAVs visit different
locations and scan more frequently. With this constraint in mind, 72 locations
evenly spread over the volume were identified, with each UAV responsible for
scanning 36 of them. The UAVs had 4 sec to fly from a location to another and
3 sec for scanning. Thus, scanning 36 locations was expected to take at least
4 min and 12 sec. Contributing to that the time required for takeoff and
landing, the UAVs were expected to operate at their operating limits.

The UAVs were controlled by a base station, i.e., a laptop running the custom
Python client software. The client was responsible for sending the UAVs to a
waypoint and instructing the scanning. Once the scanning at a waypoint was
finished, the UAV would send the results back to the client for parsing and
storing for further processing. The client was configured to be able to
control multiple UAVs with a matching set of waypoints and parameters such as
radio address, starting position, and yaw. While in the paper we show the
operation of a two-UAV system, the system can be scaled by simply adding sets
of waypoints and above-mentioned parameters.

Using this setup, data was collected for further analysis and processing. A
total of 2696 samples were collected, 1495 by UAV A and 1201 by UAV B. During
data collection, UAV A was active for 5 min and 3 sec and UAV B for precisely
5 min. In the collected samples, there were 73 distinct Medium Access Control
(MAC) addresses and 49 Service Set Identifiers, with the mean RSS of around
-73 dBm. When looking at the samples collected per UAV and scanned location
(cf., Figure 8), we see no issues with the number of samples collected by UAV
B, although this number is generally lower than for UAV A. There are
environmental factors that can play a role however i) the positive x-axis and
negative y-axis point towards the center of the apartment building where we
can expect to see more signals, ii) there is a wall segment that is 40 cm
wider where UAV B’s measurements are taken compared to UAV A, as illustrated
in Figure 8.

We expected to observe a gradual increase in the number of RSS observations
(from more APs) toward the center of the building, irrespective of which UAV
collected the samples. An illustration of this can be seen in Figure 8, which
shows a histogram per axis that groups the x and y values in bins of 0.5 m
with the height representing the number of samples collected by the UAVs in
that bin. We can clearly see that the number of samples collected increases
with an increasing x-coordinate and a decreasing y-coordinate.

###  III-B Generation of Fine-grained 3D REMs

A few pre-processing steps were taken before utilizing the data in an ML
model. Since SSIDs can be shared between devices, they were generally not
used. Instead, RSS readings were grouped based on their MAC addresses. The
timestamps were left out of consideration as well. The time difference between
the first and last collected sample was less than 10 min, hence considered as
irrelevant. MAC addresses with less than 16 samples were dropped, since the
goal was to predict RSS values of APs with sufficient number of measurements.
Finally, MAC and channel features were considered as categorical and one-hot
encoded. This pre-processing results in 2565 retained samples (131 dropped).

Since the data is locational as it represents signal quality in a 3D space, a
k-nearest neighbor regressor was considered. As features, the x, y, z
coordinates were chosen, as well as the one-hot encoded MAC addresses.
Including the one-hot encoded MAC addresses has the advantage that samples
with a different MAC address are considered farther away than similar samples
with the same MAC address. The kNN regressor was configured to use Euclidean
distance by setting $metric$=minkowski and $p$=2, as yielded by the the grid
search considering an exhaustive set of hyperparameters. Similarly, the
$weights$ and $n_{neighbors}$ parameters were tuned using a grid search where
the optimal values were $weights=distance$ and $n_{neighbors}$=3.

We have considered the samples with different MAC addresses to be further
distinguishable by multiplying the one-hot encoded values by the factor of 3
and setting the $n_{neighbors}$ parameter to 16, as yielded by the grid
search. Moreover, as an intuitive alternative to assigning samples with
different MAC addresses a greater distance, we considered a kNN estimator per
MAC address. To achieve that, we kept the hyperparameters of these MAC-based
regressors the same as previously, and took samples with the same MAC address
into account, reducing the feature set to only the x, y, z coordinates.

The last solution we considered was a neural network. In particular, we have
built and tested the network for different configurations, including; i)
multiple hidden layers with a varying amount of nodes, ii) normalized RSS
values, iii) multiple inputs: 1 for the x, y, z coordinates and 1 for the hot-
encoded MAC addresses that get combined into a common hidden layer, and iv)
different activation functions and optimizers. The optimized neural network
had an input layer for the x, y, z coordinates and the one-hot encoded MAC
addresses, sigmoid activation function, hidden layer with 16 fully connected
nodes, linear activation function, output layer with a single node for the
prediction, and Adam optimizer.

The accuracy of the considered estimators was measured based on the Root Mean
Square Error (RMSE) of their predictions. In order to have an unbiased view on
an estimators’ predictive capacity, the preprocessed data was split into a
training ($75\%$) and test ($25\%$) sets. For those estimators that require an
additional validation set for tuning their hyperparameters, the validation set
was taken out of the training set. In order to assess more elaborate
estimators we used a baseline estimator that always returns the mean per MAC
address.

Figure 8 shows a comparison of the RMSEs for the different utilized
predictors. The predictor generally utilizing the mean per MAC address
resulted in an RMSE of $4.8107$ dBm. The k-NN-based regressors generally
yielded slightly better performance than the baseline, with the best
performing one being the k-NN algorithm where the one-hot encoded MAC address
feature was multiplied by the factor of 3 and $n_{neighbors}$=16, resulting in
the RMSE of $4.4186$ dBm. Finally, the optimally tuned neural network with a
single hidden layer of 16 nodes yielded the RMSE of 4.4870 dBm. While this is
slightly better compared to the baseline, it does fall short of the best
k-nearest neighbors solution, as indicated in the figure. This collection of
regressors was used with a relatively small set of collected samples, which
explains their comparable performance. Nonetheless, we believe that this
exercise demonstrates the toolchain-like utilization of the proposed system
with different ML tools.

##  IV Conclusion

We have proposed a system for autonomous collection of 3-Dimensional (3D)
Radio Environmental Maps (REMs)-generating measurements. The main components
of the system include UWB localization-enabled UAVs acting as the carriers of
REM-generating devices. The location-annotated measurements for 3D REM
generation are then streamlined in a Machine Learning (ML) entity for fine-
grained prediction of signal quality indicators at locations not visited by
the UAVs.

Future work will include instantiating the system on a larger number of UAVs,
followed by exhaustive data collection for a set of representative
environments and different wireless technologies. These datasets will be used
for deriving the fundamental limitations on the density of 3D REMs. The UAV
localization is currently utilizing UWB, hence the REM-generating device has
to operate in a different frequency band to mitigate the self-interference
effects. This represents an obvious limitation of our system. To address this
issue, future work will focus on integrating the BitCraze’s infrared system
called Lighthouse for UAV localization, which features comparable precision,
while requiring less anchors and being cheaper. In addition to further self-
interference mitigation, this effort is expected to make the system even
easier to deploy.

## Acknowledgments

Filip Lemic was supported by the EU Marie Curie project ”Scalable
Localization-enabled In-body Terahertz Nanonetwork” (nr. 893760). In addition,
this work received support from the University of Antwerp’s Research Fund
(BOF).

## References

  * [1] H Birkan Yilmaz, Tuna Tugcu, Fatih Alagöz and Suzan Bayhan  “Radio environment map as enabler for practical cognitive radio networks”  In _IEEE Communications Magazine_ 51.12 IEEE, 2013, pp. 162–169
  * [2] Filip Lemic  “Enriched Training Database for improving the WiFi RSSI-based indoor fingerprinting performance”  In _IEEE Consumer Communications & Networking (CCNC)_, 2016, pp. 875–881  IEEE
  * [3] Serena Santi et al.  “Location-Based Vertical Handovers in Wi-Fi Networks With IEEE 802.11 ah”  In _IEEE Access_ 9 IEEE, 2021, pp. 54389–54400
  * [4] Giuseppe Caso  “Vifi: Virtual fingerprinting wifi-based indoor positioning via multi-wall multi-floor propagation model”  In _IEEE Transactions on Mobile Computing_ 19.6 IEEE, 2019, pp. 1478–1491
  * [5] Elena Meshkova  “Experimental spectrum sensor testbed for constructing indoor radio environmental maps”  In _IEEE Dynamic Spectrum Access Networks (DySPAN)_ , 2011, pp. 603–607  IEEE
  * [6] Filip Lemic  “Location-based discovery and vertical handover in heterogeneous low-power wide-area networks”  In _IEEE Internet of Things Journal_ 6.6 IEEE, 2019, pp. 10150–10165
  * [7] Daniel Romero  “Aerial spectrum surveying: Radio map estimation with autonomous uavs”  In _IEEE International Workshop on Machine Learning for Signal Processing (MLSP)_ , 2020, pp. 1–6  IEEE
  * [8] Francesco Nex and Fabio Remondino  “UAV for 3D mapping applications: a review”  In _Applied geomatics_ 6.1 Springer, 2014, pp. 1–15
  * [9] Ayon Chakraborty et al.  “SkyRAN: a self-organizing LTE RAN in the sky”  In _Proceedings of the 14th International Conference on emerging Networking EXperiments and Technologies_ , 2018, pp. 280–292
  * [10] Dimitrios Lymberopoulos  “A realistic evaluation and comparison of indoor location technologies: Experiences and lessons learned”  In _ACM Information Processing in Sensor Networks_ , 2015, pp. 178–189
  * [11] Appala Chekuri and Myounggyu Won  “Automating wifi fingerprinting based on nano-scale unmanned aerial vehicles”  In _2017 IEEE 85th Vehicular Technology Conference (VTC Spring)_ , 2017, pp. 1–5  IEEE
  * [12] Izhak Rubin and Runhe Zhang  “Placement of UAVs as communication relays aiding mobile ad hoc wireless networks”  In _MILCOM 2007-IEEE Military Communications Conference_ , 2007, pp. 1–7  IEEE
  * [13] Kai Daniel and Christian Wietfeld  “Using public network infrastructures for UAV remote sensing in civilian security operations”, 2011
  * [14] Ken Mendes  “Drone-based Autonomous Generation of 3D Maps of Indoor Radio Environments”  University of Antwerpen, https://www.uantwerpen.be/images/uantwerpen/personalpage52077/files/drone-based_indoor_rem_final.pdf, 2020
  * [15] Mark W Mueller  “Fusing UWB range measurements with accelerometers and rate gyroscopes for quadrocopter state estimation”  In _Robotics and Automation (ICRA)_ , 2015, pp. 1730–1736  IEEE

  *[3D]: 3-Dimensional
  *[ML]: Machine Learning
  *[UAVs]: Unmanned Aerial Vehicle
  *[RSS]: Received Signal Strength
  *[REM]: Radio Environmental Map
  *[SNR]: Signal-to-Noise Ratio
  *[CSI]: Channel State Information
  *[GPS]: Global Positioning System
  *[REMs]: Radio Environmental Map
  *[UWB]: Ultra Wide-Band
  *[APs]: Access Point
  *[LPD]: Loco Positioning Deck
  *[LPS]: Loco Positioning System
  *[ISM]: Industrial, Scientific, and Medical
  *[SUT]: System Under Test
  *[MAC]: Medium Access Control
  *[RMSE]: Root Mean Square Error