# RTS-GT: Robotic Total Stations Ground Truthing dataset

Maxime Vaidis1, Mohsen Hassanzadeh Shahraji, Effie Daum, William Dubois,
Philippe Giguère, François Pomerleau1

###### Abstract

Numerous datasets and benchmarks exist to assess and compare Simultaneous
Localization and Mapping (SLAM) algorithms. Nevertheless, their precision must
follow the rate at which SLAM algorithms improved in recent years. Moreover,
current datasets fall short of comprehensive data-collection protocol for
reproducibility and the evaluation of the precision or accuracy of the
recorded trajectories. With this objective in mind, we proposed the Robotic
Total Stations Ground Truthing dataset (RTS-GT) dataset to support
localization research with the generation of six- Degrees Of Freedom (DOF)
ground truth trajectories. This novel dataset includes six-DOF ground truth
trajectories generated using a system of three Robotic Total Stations tracking
moving robotic platforms. Furthermore, we compare the performance of the RTS-
based system to a Global Navigation Satellite System (GNSS)-based setup. The
dataset comprises around sixty experiments conducted in various conditions
over a period of $17$ months, and encompasses over $49$ kilometers of
trajectories, making it the most extensive dataset of RTS-based measurements
to date. Additionally, we provide the precision of all poses for each
experiment, a feature not found in the current state-of-the-art datasets. Our
results demonstrate that RTS provide measurements that are $22$ times more
stable than GNSS in various environmental settings, making them a valuable
resource for SLAM benchmark development.

Figure 1: Example of two areas provided in the dataset (bird’s-eye view).
_Left:_ a forest. _Right:_ a park on a campus. Maps in black are based on
lidar scans, while the colored spheres represent the scaled-up uncertainty of
the provided ground truth in millimeters. A Clearpath Warthog platform was
used during these experiments, whose paths are illustrated in green.

††footnotetext: 1Northern Robotics Laboratory, Université Laval, Québec City,
Canada, $\\{$maxime.vaidis, francois.pomerleau$\\}$
@norlab.ulaval.ca††footnotetext: *This research was supported by the Natural
Sciences and Engineering Research Council of Canada (NSERC) through the grant
CRDPJ 527642-18 SNOW (Self-driving Navigation Optimized for Winter).

##  I Introduction

Accurate and precise ground truth trajectories in open-source datasets are
essential to evaluate SLAM algorithms [1]. Motion capture such as Vicon or
OptiTrack systems have become the _de facto_ standard to generate such ground
truth in indoor environments [2]. Nonetheless, they are not suitable to be
deployed outside due to direct sunlight corrupting the readings and the need
to cover significantly larger areas. For outdoor deployments, the majority of
datasets use GNSS in Real Time Kinematics (RTK) mode integrated with Inertial
Navigation System (INS) to compute the reference trajectories within several
centimeters accuracy [3, 4, 5]. Meanwhile, such localization systems are
vulnerable to GNSS outages, as seen on the KITTI dataset [6]. A few datasets
rely on High-Definition (HD) maps from a terrestrial laser scanner to obtain
the reference within centimeter accuracy [7, 8]. However, surveying remains an
open challenge, especially in off-road terrains where reference planes are
scarce, making it challenging to align scans.

In recent years, a limited number of datasets used a single RTS to generate
such references indoors or outdoors with various platforms [9, 10]. RTS is a
surveying instrument capable of measuring the position of a reflective target
(hereafter _prism_) with millimeter precision in two modes: 1) static when the
prism is fixed, and 2) dynamic when the prism is in motion. Moreover, two
types of prism exist, active and passive. An active prism is equipped with
LEDs that emit a unique infrared signature. This active signature enables
multiple RTS to automatically track different prisms within their field of
view. Regardless of their accuracy and robustness, a single RTS can only
provide the reference 3D position of the tracked platform. This limitation
arises since only one prism can be tracked while the robot is in motion. By
using three passive prisms, it is possible to obtain the six-DOF reference
trajectory by collecting manually the data when the robot is static [11],
which makes the data collection procedure cumbersome. To the best of our
knowledge, active prisms were never used in any released datasets providing
six-DOF reference trajectories generated by RTS in dynamic mode. Moreover, new
research done recently enables estimation of the uncertainty for a multi-RTS
setup [12]. This uncertainty on ground truth is not available in state-of-the-
art SLAM datasets, which can lead to unbiased comparisons as many SLAM
algorithms are approaching centimeter-level accuracy [13].

The key motivation for this work is to provide a high-quality dataset that
covers a variety of challenging environments and further motivates ground
truth generation for SLAM algorithms. As shown in our previous research [14,
12], generating a reliable six-DOF reference trajectories of a moving robot
with three RTS is now feasible. As a follow-up, we present RTS-GT, a dataset
providing six-DOF reference trajectories originating from two different types
of setup, three RTS and three GNSS receivers. The RTS-GT dataset was collected
during $17$ months in diverse weather conditions, totaling over $49$
kilometers of trajectories. With this dataset, we provide and compare the
precision and the reproducibility of the setups during multiple weather
conditions and in different environments as shown in Figure 1, demonstrating
that RTS are more reliable than GNSS to generate ground truth trajectories.

##  II Related work

TABLE I: Comparison of public datasets containing ground truth trajectories. The top half of the table shows the most popular datasets used for SLAM evaluation. The bottom half shows public datasets that contain ground truth generated by the most accurate setups. The RTS-GT dataset is the only that allows dynamic six-DOF ground truth generation with RTSs on a large scale. The symbol ∗ indicates the value is estimated by us. The letter _O_ means _Outside_ and the letter _I_ means _Inside_. _Dataset_ | _Area_ |  _Platform_ | _Sensors_ |  _Ground truth_ | _DOF_ |  _Accuracy_ | _Distance_ |  _Weather_ | _Site_
---|---|---|---|---|---|---|---|---|---
Lidar |  Cam. | IMU |  GNSS | Radar
KAIST [3] | Urban |  Car | ✓ |  ✓ | ✓ |  ✓ | - |  INS | 3 |  $<$ $10\text{\,}\mathrm{cm}$∗ | $191\text{\,}\mathrm{km}$ |  - | O
Wild-Places [15] | Forest |  Handheld | ✓ |  ✓ | ✓ |  - | - |  HD map | 6 |  $<$ $30\text{\,}\mathrm{cm}$ | $33\text{\,}\mathrm{km}$ |  - | O
NCLT [5] | Campus |  Segway | ✓ |  ✓ | ✓ |  ✓ | - |  INS | 6 |  $\thickapprox$ $10\text{\,}\mathrm{cm}$∗ | $147\text{\,}\mathrm{km}$ |  - | I/O
KITTI [6] | Urban |  Car | ✓ |  ✓ | ✓ |  ✓ | - |  INS | 6 |  $<$ $10\text{\,}\mathrm{cm}$∗ | $39\text{\,}\mathrm{km}$ |  - | O
nuScenes [16] | Urban |  Car | ✓ |  ✓ | ✓ |  ✓ | ✓ |  HD map | 6 |  $<$ $10\text{\,}\mathrm{cm}$ | $242\text{\,}\mathrm{km}$ |  Rain | O
Oxford [4] | Urban |  Car | ✓ |  ✓ | ✓ |  ✓ | ✓ |  INS | 6 |  $\thickapprox$ $1\text{\,}\mathrm{cm}$ | $280\text{\,}\mathrm{km}$ |  Fog, rain, snow | O
Hilti-Oxford [8] | Urban |  Handheld | ✓ |  ✓ | ✓ |  - | - |  HD map | 6 |  $<$ $1\text{\,}\mathrm{cm}$ |  $<$ $10\text{\,}\mathrm{km}$ |  - | I/O
Hilti SLAM [10] | Urban |  Handheld | ✓ |  ✓ | ✓ |  - | - |  Hilti PLT 300 | 3 |  $\thickapprox$ $3\text{\,}\mathrm{mm}$ |  $<$ $10\text{\,}\mathrm{km}$ |  - | I/O
Euroc [17] | Indoor |  UAV | - |  ✓ | ✓ |  - | - |  Leica MS50 | 3 |  $\thickapprox$ $1\text{\,}\mathrm{mm}$ | $0.9\text{\,}\mathrm{km}$ |  - | I
UZH-FPV [9] | Campus |  UAV | - |  ✓ | ✓ |  - | - |  Leica MS60 | 3 |  $\thickapprox$ $1\text{\,}\mathrm{mm}$ | $10\text{\,}\mathrm{km}$ |  - | I
TUM [2] | Campus |  Handheld | - |  ✓ | ✓ |  - | - |  OptiTrack | 6 |  $\thickapprox$ $1\text{\,}\mathrm{mm}$ |  $<$ $1\text{\,}\mathrm{km}$* |  - | I/O
Challenging dataset [11] | Campus, mountain |  Rig | ✓ |  - | ✓ |  - | - |  Leica TS15 | 6 |  $\thickapprox$ $1\text{\,}\mathrm{mm}$ |  $<$ $1\text{\,}\mathrm{km}$ |  - | I/O
_RTS-GT dataset (Ours)_ | Campus, forest |  UGV | ✓ |  - | ✓ |  ✓ | - |  Trimble S7 | 6 |  $\thickapprox$ $4\text{\,}\mathrm{mm}$ | $49\text{\,}\mathrm{km}$ |  Fog, rain, snow | I/O

Among the majority of datasets used for benchmarking SLAM algorithms, the
prominent reference trajectory generation method is GNSS-Aided INS. This
approach relies on integrating data from one or more GNSS receivers with high-
frequency measurements obtained from a Inertial Measurement Unit (IMU) within
an INS [18]. This fusion provides the six-DOF pose estimation for the
platform. The KITTI dataset [6] was the first dataset to introduce a SLAM
algorithm benchmark using over $39\text{\,}\mathrm{km}$ of reference
trajectory data, generated by an GNSS-Aided INS. This system was mounted on a
car equipped with lidar and camera sensors for continuous data collection.
Given this success within the scientific community, numerous other benchmarks
emerged, employing sensor-equipped cars and similarly utilizing GNSS-Aided INS
for production of reference trajectories datasets [3, 4, 5]. Nonetheless,
GNSS-Aided INS systems encounter challenges in urban environments, such as the
GNSS canyon effect caused by buildings, issues related to signal multi-path,
and limited satellite visibility. As a result, the system’s accuracy can
fluctuate, typically ranging around $10\text{\,}\mathrm{cm}$ which may be too
significant for specific benchmarking scenarios such as the TUM-VI dataset
[2]. This dataset focused on high-speed drone localization using camera images
and IMU data. Given the drone’s aggressive dynamics, a one-millimeter-accurate
Optitrack system was used to generate the six-DOF reference trajectory. Other
datasets have employed high-definition 3D maps provided by terrestrial
scanners to generate reference trajectories [16]. The obtained accuracy is
under ten centimeters, and the dataset supports object detection and tracking
to train autonomous vehicle algorithms under clear or rainy conditions. More
recently, [8] also utilizes high-resolution mapping of urban locations to
reconstruct a camera, lidar, and IMU setup’s trajectory in six-DOF. Lidar
scans are directly aligned with the high-definition 3D map, resulting in
millimeter-level accuracy. For natural environments, [15] used an SLAM
algorithm to reconstruct their setup’s trajectory in six-DOF and used it as a
ground truth with an $30\text{\,}\mathrm{cm}$ accuracy. All of these methods
are limited in accuracy, as seen with GNSS-Aided INS systems, or are difficult
to apply in outdoor environments, such as high-resolution terrestrial laser
scanning or Optitrack measurement methods. Our RTS-GT dataset overcomes these
limitations, offering a dataset in various indoor and outdoor environments
with centimeter-level accuracy through the use of multiple RTS.

As mentioned earlier, previous datasets have already employed RTS to generate
reference trajectories for their robotic platforms. The EuRoc [17] and UZH-FPV
[9] datasets have employed a Leica RTS-based system to generate reference
trajectories for drones navigating in indoor and outdoor environments.
Similarly, the Hilti SLAM dataset [10] uses a Leica RTS-based system for
dynamic tracking of a setup equipped with lidar, cameras, and an IMU. To
generate reference trajectories using RTS, two strategies exist in the
literature. The first strategy, employed by [11], involves placing three
passive prisms on a robotic platform and manually capturing static position
measurements using an RTS. This method provides the platform’s six-DOF with
millimeter-level accuracy. However, the number of obtained reference poses is
limited due to the time-consuming nature of manual RTS measurements. The
second strategy involves dynamically tracking an active prism with an RTS. The
positional measurement accuracy is in the order of millimeters, and the entire
prism trajectory can be reconstructed. However, the generated reference
trajectory has only three-DOFs as a single prism is tracked in real-time.

The Table I provides an overview of the mentioned datasets along with their
respective characteristics. To the best of our knowledge, there is not any
reference dataset that provides a six-DOF reference trajectory of a dynamic
platform using a multi-RTS-based system. Moreover, there is a lack of
precision information in available datasets regarding the precision of the
dynamic RTS or GNSS measurements. Hence, we present the RTS-GT dataset, the
first dataset focusing on creating six-DOF reference trajectories using three
RTS tracking three distinct active prisms. Additionally, we provide the
precision of measurements and the final poses, a unique feature not previously
provided in other datasets.

##  III Hardware

Two mobile robotic platforms were deployed for the dataset: a Clearpath
Warthog Uncrewed Ground Vehicle (UGV), and a SuperDroid HD2 UGV. The Warthog
robot is equipped with a RoboSense RS-32 lidar and an XSens MTi-10 IMU,
whereas the HD2 robot’s sensor payload consists of a Velodyne-VLP 16 lidar and
an XSens MTi-30 IMU. We generalized our approach to generate reference
trajectories by using and comparing two different types of setup for the two
robots, as shown in Figure 2 and Figure 3 for deployments of the HD2 and
Warthog platforms, respectively.

Figure 2: Setup used in the tunnel sites. _Left_ : RTS setup with the HD2
robot. _Right_ : one deployment done in a $120\text{\,}\mathrm{m}$ tunnel.
Because the floor was slippery, heavy weights were added to stabilize the RTS
tripods.

The first ground truth setup is composed of three Trimble S7 RTS. Each RTS
tracks a single Trimble MultiTrack Active Target MT1000 prism at a maximum
achievable measurement rate of $2.5\text{\,}\mathrm{Hz}$. For this active
prism, the nominal range at the tracking mode is $800\text{\,}\mathrm{m}$ and
the nominal position measurement accuracy is $4\text{\,}\mathrm{mm}$. To
collect the data from each RTS, a Raspberry Pi 4 is used as a client through a
USB connection. The data are then sent to a master Raspberry Pi 4 through LoRa
Shield radio modules from Dragino operating with a radio frequency of
$905\text{\,}\mathrm{MHz}$. The chosen modulation allows the Raspberry Pi to
reliably send data over a distance of up to $800\text{\,}\mathrm{m}$ at
$366\text{\,}\mathrm{Bps}$ in open space. The second ground truth setup is
composed of a set of four GNSS receivers used for outdoor experiences. In this
setup, three GNSS receivers are mounted on the Warthog UGV and the fourth GNSS
receiver serves as a static base station. For the sake of comparison, two
different models of GNSS receivers are used in this dataset, _Reach RS+_ and
_Trimble R10-2_. When operating in RTK rover/base mode, the _Trimble R10-2_
receiver has a vertical accuracy of $21.2\text{\,}\mathrm{mm}$ and horizontal
accuracy of $11.3\text{\,}\mathrm{mm}$, while the _Reach RS+_ has a vertical
accuracy of $19.8\text{\,}\mathrm{mm}$ and horizontal accuracy of
$9.9\text{\,}\mathrm{mm}$. All three prisms and GNSS receivers were mounted on
the Warthog robot for outdoor experiments. The HD2 UGV was used to collect
data in the university tunnels. In such environments, GNSS receivers are not
functional, thus, only the prisms can be used.

Figure 3: Setup used on the campus with the Warthog UGV. A GNSS fixed base
station sends corrections to the three GNSS rovers on the robot. Three prisms
are tracked by three RTS. Data are collected by three Raspberry Pi clients
connected by USB to the RTS. A LoRa communication protocol is used to send
data to a Raspberry Pi master located on the UGV. The lidar and IMU are on the
front part of the Warthog.

##  IV Data collection

Our main contribution revolves around delivering a dataset featuring openly
available ground truth trajectories obtained using diverse setups. The RTS-GT
dataset was gathered across three distinct environmental settings,
encompassing various times of day and spanning around diverse weather
conditions, summarized in the Table II.

TABLE II: Table of deployments in the RTS-GT dataset. The weather legend is as follows: C=Clear, FR=Freezing rain, S=Snow, R=Rain. The robot legend means W=Warthog and H=HD2. | _Month_ |  _Exp._ | _Length_ |  _Weather_ | _Robot_ |  _Setup_
---|---|---|---|---|---|---
| Feb. 22 |  1 | $1.58\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS
| Mar. 22 |  6 | $7.45\text{\,}\mathrm{km}$ |  FR, S, C | W |  RTS/GNSS
| May 22 |  8 | $14.43\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS
Campus | Jun. 22 |  4 | $3.42\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS
| Jul. 22 |  6 | $2.67\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS
| Nov. 22 |  14 | $3.91\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS
| Dec. 22 |  4 | $2.1\text{\,}\mathrm{km}$ |  C, R | W |  RTS/GNSS
| Jul. 23 |  2 | $2\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS
| May 22 |  6 | $1.55\text{\,}\mathrm{km}$ |  C | H |  RTS
Tunnel | Jul. 22 |  5 | $1.58\text{\,}\mathrm{km}$ |  C | H |  RTS
| Sep. 22 |  9 | $0.85\text{\,}\mathrm{km}$ |  C | H |  RTS
Forest | Nov. 22 |  4 | $3.16\text{\,}\mathrm{km}$ |  C | W |  RTS/GNSS

###  IV-A Environments

The first deployment site is the campus of Université Laval in Quebec City.
This campus comprises buildings, open spaces, and wooded areas. Multiple views
of the campus are depicted in Figure 1-_right_ and Figure 3. Overall, $45$
different experiments through $20$ distinct deployments were conducted with
the Warthog UGV, totaling $37.6\text{\,}\mathrm{km}$ of trajectories on the
campus. During these data collection campaigns, various weather conditions
were encountered, including clear weather, rain, and even a snowstorm. Data
from both RTS and GNSS setups are available for each of the conducted
experiments, enabling a comparison of their respective accuracies in
generating reference trajectories.

For the second deployment site, underground tunnels located beneath the
university campus were selected. Four deployments were conducted in these
tunnels to gather data from RTS setup during $20$ different experiments,
covering a total of $4\text{\,}\mathrm{km}$ of trajectories with the HD2
platform. These tunnels have lengths of several hundred meters, which might
cause particular challenges for SLAM algorithms. Unlike other setups based on
motion capture or ultra-wideband, our RTS setup can generate reference
trajectories at long ranges without the need to alter the environment. The
Figure 2 shows an example of a tunnel as a long-range environment.

Finally, the last site is located in the Montmorency Forest, which belongs to
Université Laval. This site contains numerous paths for snowmobiles and cross-
country ski trails in a dense forest. Two deployments were carried out with
the Warthog robot, along with RTS and GNSS systems. In total,
$3.2\text{\,}\mathrm{km}$ of trajectories were acquired during four different
experiments. An example of the RTS system in the forest is depicted in Figure
1-_left_. The lidar was available during the majority of deployments to
collect data for SLAM algorithms.

###  IV-B Ground truth protocol

In this section, we introduce the standardized protocols that we employed for
each of the setups during all conducted deployments, split between RTS and
GNSS. Starting with RTS protocol, the field deployment of the three RTS is
carried out in several steps:

  1. 1.

($30\text{\,}\mathrm{min}$) RTS units are acclimated to the ambient
temperature before data collection can begin. This step is necessary to
prevent any condensation effects that could bias the measurements, especially
in winter;

  2. 2.

($5$-$60\text{\,}\mathrm{min}$) Tripods for the RTS are set up, while the RTS
units adjust to the ambient temperature;

  3. 3.

($10$-$15\text{\,}\mathrm{min}$) RTS units are mounted on tripods, and we
roughly level RTS units visually. To achieve a finer leveling for RTS units,
calibrated electronic sensors are utilized;

  4. 4.

($3$-$5\text{\,}\mathrm{min}$) Raspberry Pi devices are powered on and
connected via USB to the RTS units to retrieve their data and send it to the
master unit for recording. These Raspberry Pi devices also assign a unique
prism number to each RTS unit to be tracked;

  5. 5.

($1$-$2\text{\,}\mathrm{min}$) Active prisms are mounted on the UGV and
powered on with their unique ID. To facilitate data processing of multiple
experiments conducted by the robots, the same prism ID and positions were
adopted during all data acquisition operations;

  6. 6.

($1$-$10\text{\,}\mathrm{min}$) Verification of proper prism tracking is
performed, followed by measurements of static prism positions to enable post-
processing extrinsic calibration of RTS units.

The complete setup process takes between $50\text{\,}\mathrm{min}$ and
$120\text{\,}\mathrm{min}$, which depends on the number of available operators
and weather conditions. At the end of the data collection procedure, prisms
have to remain on the UGV to immediately perform an extrinsic sensor
calibration, discussed later.

As for the GNSS protocol, we need a minimum of three GNSS receivers as rover
plus one GNSS receiver as a base station to generate a six-DOF ground truth
with only GNSS-based setup. This rover/base configuration, known as RTK, can
achieve centimeter-level accuracy by correcting errors of the GNSS receivers
used as rovers. The base station, consisting of a GNSS receiver, provides
corrections for the rover GNSS observations by simultaneously monitoring the
same satellites as the rover receivers. The base station can be fixed at a
predetermined, known, and stationary location (e.g., a geodesic pillar stock
or a geodetic survey marker) or an unknown position. If the position of the
base station is known in advance, GNSS receivers have just to be powered on
and once the radio connection between the base and rovers is established the
system is operational and the collection of data can start. On the other hand,
if the position of the base station is unknown, a waiting time of at least
$15\text{\,}\mathrm{min}$ is necessary before data collection to let the rover
or the base receivers have enough time to boot, average their positions, and
achieve optimal reading of the satellite constellations in the sky. The
positional corrections are then transmitted through messages via a radio link
from the base station to the rover receivers, where they are employed to
correct the real-time positions of the rover. Moreover, the internal radio of
the GNSS receiver has a maximum range of $2\text{\,}\mathrm{km}$. To achieve a
greater range, we have to use the external radio, which can theoretically
transmit up to $10\text{\,}\mathrm{km}$ under optimal conditions.

###  IV-C Calibration and time synchronization

This section specifies how the systems were calibrated and synchronized, along
with the data format given by the dataset. Attention is given to 1) extrinsic
calibration and 2) time synchronization. First, the extrinsic calibration
process involves obtaining the precise pose of all the sensors on the robots.
Calibration is performed using one RTS at the end of each deployment to
closely match the conditions of the experiments as temperature, pressure, and
humidity can affect measurements, especially in winter. Retro-reflective
targets, as shown in Figure 3, are stuck atop the lidar and GNSS receivers.
The prisms are left in active mode for this calibration. Ten repetitions are
performed for each millimeter-precise position to increase precision and to
determine their uncertainties. This method provides the relative positions of
the prisms, GNSS receivers, and lidar to each other. Since the IMU is too
small, its extrinsic calibration with the lidar is based on the work of [19],
who uses lidar and IMU data to perform calibration within a tenth of a degree
using a modified four-DOF SLAM algorithm. Translation between the lidar and
IMU is determined using the Computer-Aided Design (CAD) model of their
support. These calibrations were performed after each deployment to ensure
precise referenced measurements. As for time synchronization, a Network Time
Protocol (NTP) daemon was used to synchronize the clocks between the two on-
board computers on the robot (i.e., the main computer of the UGV and the
Raspberry Pi master connected to the robot network). Ten minutes was allowed
after booting up both of the computers so that the client clock could adjust.
Measurements from the wheel encoders, lidar, and IMU were timestamped using
the data logging computer’s clock. Synchronization between the master and
client Raspberry Pi is achieved using a modified NTP protocol designed for
LoRa communication, as described in [20]. Initial synchronization is performed
at the start of data collection, followed by repetitions every five minutes
for each client. This method ensures clock precision at the level of one to
two milliseconds. As the GNSS devices are not connected to the setup, temporal
synchronization to the RTS data is established using a classical maximum
likelihood state estimator, similar to the approach taken by [17].

###  IV-D Data format

For each deployment and experiment, ROS 2 rosbag files are provided for data
gathered by the Warthog and HD2 UGV. These files include lidar scans, IMU
measurements, motors and encoders data, as well as rigid transformation
between all sensor frames. Unified Robot Description Format (URDF) files are
provided along with these rosbag files for each UGV. Note that no INS were
used to process the IMU measurements. All IMU raw data is available in the
different rosbag files. With the RTS setup, prism positions, and time
synchronizations are provided for rosbag files in ROS 1 and ROS 2. The static
prism positions computed for the RTS extrinsic calibration are given in a text
file GCP.txt, and the sensor extrinsic calibration results are given in
another text file calibration_raw.txt. Finally, NMEA and UBX files are
provided for all GNSS receivers used during experiments. The RTS-GT dataset is
available at https://github.com/norlab-ulaval/RTS_project, along with a
toolbox code to process the data.

##  V Discussion and challenges

To assess the disparities among the setups, we conducted an analysis of $27$
outdoor experiments, covering a total distance of $20.6\text{\,}\mathrm{km}$,
during which we tracked the trajectories of the robot using GNSS and RTS-based
setups when available. We evaluated the precision of each system by employing
_inter-distances_ between prisms and GNSSs. These distances are calculated
between every synchronized triplet of prism positions or GNSS receiver
positions recorded during each experiment, respectively referred as inter-
prism distances and inter-GNSS distances. Subsequently, each of these distance
triplets is compared to their corresponding calibrated distance (evaluated at
step six of our protocol) and taken as a reference to obtain the errors in
precision. Furthermore, we employed an inter-precision distance to quantify
variation in precision between different experiments conducted on the same
site at different times. The purpose of this distance is to highlight the
reproducibility of the different setups to generate ground truth trajectories.
Sets of closest positions in the ground truth trajectories recorded during
separate experiments are computed by a nearest neighbors algorithm.
Subsequently, the inter-precision distance is computed for each set by
subtracting each average inter-prism or inter-GNSS distances of each position.

###  V-A Precision and reproducibility

Results of inter-prism distances, inter-GNSS distances, and precision on the
final translation and rotation of the robot are shown in Table III for $15$
deployments. The precision on the final pose is estimated by the same method
used in our previous work [12]. It can be seen that RTS precision is stable in
multiple environments, while the GNSS can have a variation of more than
$300\text{\,}\mathrm{m}\mathrm{m}$ in the forest environment, as well as in
open space. Moreover, higher distances between a prism and its RTS affect the
precision, leading to higher uncertainties.

TABLE III: Median values for different metrics obtained over $15$ RTS and GNSS setups. All deployments happened during 2022 and the date format is day/month. For each row, light colors indicate low values, whereas dark means high values. | _Environment_ |  |  |  |  |  |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1pt]2-Z | _Open space_ |  |  |  |  |  |  |  |  |  _Building_ |  |  _Forest_ |  |  |
|  24/02 |  07/03 |  14/03 |  16/03 |  22/06 |  30/06 |  11/07 |  29/11 |  05/12 |  12/03 |  16/11 |  31/03 |  09/11 |  10/11 |  24/11
[0.5pt] Inter-prism distance [mm] | 6.5 | $2.4$ |  4.1 | 4.1 |  9.7 | 5.7 |  22.9 |  4.6 |  3.7 |  4.3 |  4.5 |  6.6 |  18.0 |  5.0 |  3.2
Inter-GNSS distance [mm] | 5.5 | $5.2$ |  4.6 | 6.4 |  37.2 | 368.0 |  111.0 |  6.9 |  20.5 |  230.0 |  17.7 |  288.0 |  149.0 |  394.0 |  423.0
Translation error [mm] | 15.2 | $28.7$ |  25.3 | 10.6 |  12.3 | 12.8 |  35.2 |  22.9 |  12.7 |  23.7 |  12.0 |  19.0 |  12.8 |  6.7 |  21.0
Rotation error [deg] | 1.51 | $1.6$ |  2.05 | 1.4 |  0.9 | 0.94 |  4.0 |  1.94 |  1.33 |  1.93 |  0.85 |  1.06 |  0.97 |  0.53 |  1.59
RTS Range [m] | 39.0 | $20.0$ |  18.0 | 19.0 |  90.0 | 32.0 |  87.0 |  39.0 |  32.0 |  71.0 |  42.0 |  38.0 |  42.0 |  59.0 |  37.4

The results depicted in Figure 4-_Left_ offer general insights into the errors
associated with the inter-prism and inter-GNSS distances. These findings
reveal that the RTS acquisition system consistently achieves a median
precision of approximately $4.5\text{\,}\mathrm{mm}$, whereas the GNSS system
exhibits a $22$ times lower median precision, being around
$118.1\text{\,}\mathrm{mm}$. It is important to emphasize that the inter-
distance errors highlight the highest precision of the RTS acquisition system
compared to the GNSS system. This discrepancy can be attributed to the
relatively low error inherent to the RTS acquisition system, in contrast to
the absolute error associated with GNSS.

Figure 4:  Distribution of errors arising for the two setups. _Left_ : inter-
prism and inter-GNSS distances. _Right_ : inter-precision distances. Results
obtained from the RTS are denoted in blue, whereas those from the GNSS are
shown in orange. The median error values are in the center of each box, and
the IQR is indicated alongside for reference.

To evaluate the reproducibility between experiments, we employ a nearest
neighbor distance calculation within three meters, with data expressed in the
GNSS frame. As illustrated in Figure 4-_Right_ , the RTS setup consistently
exhibits high reproducibility, with a median error of
$4.2\text{\,}\mathrm{mm}$. The GNSS have the same level of reproducibility,
with a median of $6.2\text{\,}\mathrm{mm}$. However, the GNSS IQR, with a
value of $133\text{\,}\mathrm{mm}$, is $14$ times more important than the RTS.
This observation underscores the stability of precision across all experiments
for the RTS compared to the GNSS.

###  V-B Challenges encountered

The first encountered challenge was the leveling of the RTS in winter. It is
sometimes necessary to remove snow from the ground for step two of our
protocol to prevent tripods from gradually sinking into the snow, which would
impact the ground truth. In forested areas where the snow depth can reach
several meters, as shown in Figure 5-_Left_ , surface snow is compacted to
provide stability. Secondly, as depicted in Figure 1, obstacles between prisms
and RTS can disrupt measurements. To address this issue, RTS placement
locations are pre-selected based on the planned trajectory of the UGV, and the
heights of both RTS and prisms are varied to minimize occlusion risks. Another
challenge is the presence of dust or dirt on RTS lenses, as shown in Figure
5-_Middle_ , which can degrade RTS performance. Therefore, lenses should be
wiped regularly. Furthermore, since prisms and GNSS units are elevated on the
UGV, they are susceptible to vibrations during motions. These vibrations can
lead to positioning errors of up to $1\text{\,}\mathrm{cm}$. To deal with this
problem, metal supports have been added to dampen vibrations, as depicted in
Figure 5-_Right_. Finally, unlike GNSS data, RTS measurement timestamps are
not globally valid. To address this issue and obtain the six-DOF pose, data
interpolation is performed. This interpolation can reduce the accuracy of the
final estimated pose, especially for high UGV dynamic and low-rate measurement
setups.

Figure 5: Issues encountered while collecting the dataset. _Left_ :
difficulties of leveling on deep snow. _Middle_ : dust on a lens which was
interfering with the tracking mode of the RTS. _Right_ : supports added on the
Warthog to reduce vibrations on the top prism.

##  VI Conclusion

In this paper, we have introduced a novel dataset, the RTS-GT dataset,
designed to compare various ground truthing setups. The RTS-GT dataset stands
out as a unique dataset for providing six-DOF reference trajectories of a
moving robotic platform generated by three RTS. The dataset encompasses GNSS
and RTS data as well, along with information from lidar, IMU, and encoder
sensors. It facilitates the application of SLAM algorithms and the assessment
of results to the ground truth data. Furthermore, the RTS-GT dataset covers
diverse environments and weather conditions to assess the quality of the
different ground truth setups. Additionally, tools for quantifying their
precision are provided, a feature not present in previous datasets. An
extensive analysis of the precision and reproducibility of different
trajectory generation setups was conducted. The results indicate that RTS
systems deliver more precise and reproducible data compared to GNSS solutions,
even when used in RTK mode. These results demonstrate the challenges posed by
generating six-DOF ground truth trajectories in outdoor and indoor
environments. Moreover, it shows that multiple RTS can be used to benchmark
six-DOF SLAM algorithm comparisons.

## References

  * [1] Li Qingqing, Yu Xianjia, Jorge Peña Queralta and Tomi Westerlund  “Multi-Modal Lidar Dataset for Benchmarking General-Purpose Localization and Mapping Algorithms”  In _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 3837–3844  DOI: 10.1109/IROS47612.2022.9981078
  * [2] David Schubert et al.  “The TUM VI Benchmark for Evaluating Visual-Inertial Odometry”  In _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 1680–1687  DOI: 10.1109/IROS.2018.8593419
  * [3] Yukyung Choi et al.  “KAIST Multi-Spectral Day/Night Data Set for Autonomous and Assisted Driving”  In _IEEE Transactions on Intelligent Transportation Systems_ 19.3, 2018, pp. 934–948  DOI: 10.1109/TITS.2018.2791533
  * [4] Dan Barnes et al.  “The Oxford Radar RobotCar Dataset: A Radar Extension to the Oxford RobotCar Dataset”  In _IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 6433–6438  DOI: 10.1109/ICRA40945.2020.9196884
  * [5] Nicholas Carlevaris-Bianco, Arash K. Ushani and Ryan M. Eustice  “University of Michigan North Campus long-term vision and lidar dataset”  In _International Journal of Robotics Research_ 35.9, 2015, pp. 1023–1035
  * [6] Andreas Geiger, Philip Lenz and Raquel Urtasun  “Are we ready for autonomous driving? The KITTI vision benchmark suite”  In _IEEE Conference on Computer Vision and Pattern Recognition_ , 2012, pp. 3354–3361  DOI: 10.1109/CVPR.2012.6248074
  * [7] Kamak Ebadi et al.  “Present and Future of SLAM in Extreme Environments: The DARPA SubT Challenge”  In _arXiv preprint arXiv:2208.01787_ , 2022  arXiv: http://arxiv.org/abs/2208.01787
  * [8] Lintong Zhang et al.  “Hilti-Oxford Dataset: A Millimeter-Accurate Benchmark for Simultaneous Localization and Mapping”  In _IEEE Robotics and Automation Letters_ 8.1, 2023, pp. 408–415  DOI: 10.1109/LRA.2022.3226077
  * [9] Jeffrey Delmerico et al.  “Are We Ready for Autonomous Drone Racing? The UZH-FPV Drone Racing Dataset”  In _International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 6713–6719  DOI: 10.1109/ICRA.2019.8793887
  * [10] Michael Helmberger et al.  “The Hilti SLAM Challenge Dataset”  In _IEEE Robotics and Automation Letters_ 7.3, 2022, pp. 7518–7525  DOI: 10.1109/LRA.2022.3183759
  * [11] François Pomerleau, Ming Liu, Francis Colas and Roland Siegwart  “Challenging data sets for point cloud registration algorithms”  In _The International Journal of Robotics Research_ 31.14, 2012, pp. 1705–1711  DOI: 10.1177/0278364912458814
  * [12] Maxime Vaidis et al.  “Uncertainty analysis for accurate ground truth trajectories with robotic total stations”  In _IEEE International Conference on Intelligent Robots and Systems (IROS)_ , 2023
  * [13] Carlos Campos, Jose M.M. Montiel and Juan D. Tardos  “Inertial-Only Optimization for Visual-Inertial Initialization”  In _IEEE International Conference on Robotics and Automation (ICRA)_ IEEE, 2020  DOI: 10.1109/icra40945.2020.9197334
  * [14] Maxime Vaidis et al.  “Extrinsic calibration for highly accurate trajectories reconstruction”  In _IEEE International Conference on Robotics and Automation (ICRA)_ , 2023, pp. 4185–4192  DOI: 10.1109/ICRA48891.2023.10160505
  * [15] Joshua Knights et al.  “Wild-Places: A Large-Scale Dataset for Lidar Place Recognition in Unstructured Natural Environments”  In _IEEE International Conference on Robotics and Automation (ICRA)_ , 2023, pp. 11322–11328  DOI: 10.1109/ICRA48891.2023.10160432
  * [16] Holger Caesar et al.  “nuScenes: A Multimodal Dataset for Autonomous Driving”  In _IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 11618–11628  DOI: 10.1109/CVPR42600.2020.01164
  * [17] Michael Burri et al.  “The EuRoC micro aerial vehicle datasets”  In _The International Journal of Robotics Research_ 35.10, 2016, pp. 1157–1163  DOI: 10.1177/0278364915620033
  * [18] Mohinder S Grewal, Lawrence R Weill and Angus P Andrews  “Global positioning systems, inertial navigation, and integration”  John Wiley & Sons, 2007
  * [19] Vladimír Kubelka, Maxime Vaidis and François Pomerleau  “Gravity-constrained point cloud registration”  In _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 4873–4879  DOI: 10.1109/IROS47612.2022.9981916
  * [20] Maxime Vaidis, Philippe Giguère, François Pomerleau and Vladimír Kubelka  “Accurate outdoor ground truth based on total stations”  In _18th Conference on Robots and Vision (CRV)_ , 2021, pp. 1–8  DOI: 10.1109/CRV52889.2021.00012

  *[SLAM]: Simultaneous Localization and Mapping
  *[RTS-GT]: Robotic Total Stations Ground Truthing dataset
  *[DOF]: Degrees Of Freedom
  *[RTS]: Robotic Total Station
  *[GNSS]: Global Navigation Satellite System
  *[RTK]: Real Time Kinematics
  *[INS]: Inertial Navigation System
  *[HD]: High-Definition
  *[IMU]: Inertial Measurement Unit
  *[DOFs]: Degrees Of Freedom
  *[UGV]: Uncrewed Ground Vehicle
  *[CAD]: Computer-Aided Design
  *[NTP]: Network Time Protocol
  *[URDF]: Unified Robot Description Format
  *[GNSSs]: Global Navigation Satellite System
  *[IQR]: Interquartile Range


# Universal patterns of skyrmion magnetizations unveiled by defect
implantation

Imara Lima Fernandes  Samir Lounis

###### Abstract

Skyrmions are spin-swirling textures hosting wonderful properties with
potential implications in information technology. Such magnetic particles
carry a magnetization, whose amplitude is crucial to establish them as robust
magnetic bits, while their topological nature gives rise to a plethora of
exquisite features such as topological protection, the skyrmion and
topological Hall effects as well as the topological orbital moment. These
effects are all induced by an emergent magnetic field directly proportional to
the three-spin scalar chirality,
$\chi=(\mathbf{S}_{i}\times\mathbf{S}_{j})\cdot\mathbf{S}_{k}$, and shaped by
the peculiar spatial dependence of the magnetization. Here, we demonstrate the
existence of novel chiral magnetizations emerging from the interplay of spin-
orbit interaction and either $\chi$ or the two-spin vector chirality
$\boldsymbol{\kappa}=\mathbf{S}_{i}\times\mathbf{S}_{j}$. By scrutinizing
correlations among the spin, orbital (trivial and chiral) magnetizations, we
unveil from ab-initio universal patterns, quantify the rich set of
magnetizations carried by single skyrmions generated in PdFe bilayer on
Ir(111) surface and demonstrate the ability to engineer their magnitude via
controlled implantation of impurities. We anticipate that our findings can
guide the design of disruptive storage devices based on skyrmionic bits by
encoding the desired magnetization with strategic seeding of defects.

## Introduction

The identification of universal patterns characterizing properties of
materials is crucial for various fields of science. This enables to design
materials with the right combination of elements to achieve desired
characteristics. A fundamental example is the Mendeleev periodic table that
consists of a graphical representation of chemical elements, which depending
on their atomic numbers and location in the table share trends with periodic
patterns, such as electron negativity, affinity, metallicity and cohesion. In
the context of magnetism, Hund’s rules provide generic dependencies of the
spin and orbital moments as function of the filling of electronic states. The
spin moment follows an inverse parabolic behavior with a maximum expected at
half-filling of the $d$-states, while the orbital moment exhibits an S-shape
passing through zero at half-filling with early (late) transition elements
having an orbital moment antiferromagnetically (ferromagnetically) aligned to
the spin moment. The Slater-Pauling rule [1, 2] illustrates another primary
general prescription in magnetism, which uses the electronic filling as a
fundamental descriptor dictating how the total magnetization of an alloy is
shaped by the amount of valence electrons injected to or removed upon
alteration of the concentration of the atomic building-blocks. Instead of the
electronic filling, the radius of the 3$d$ electron shell acts as the
descriptor defining the celebrated Bethe-Slater curve [3, 4, 5], a common
textbook example, which is fundamental to understand the magnetic properties
and magnetic ordering of transition metal materials.

In this work, we unveil a universal behavior characterizing the magnetization
carried by complex spin-textures. We exemplify our findings by addressing
magnetic skyrmions [6, 7], which are topological objects confined in two-
dimension [8] and heavily prospected as a mean of encoding magnetic bits for
information technology [9, 10]. Here the magnetic bit is defined by the
magnetization carried by the skyrmion, which is different from the surrounding
background. Clearly the ability to engineer the amplitude of the magnetization
is paramount for potential device applications. Moreover, it is the peculiar
shape of the magnetization that defines the emergent magnetic field giving
rise to the topological properties characterizing skyrmions, such as the
topological charge as well as the topological and skyrmion Hall effects [8].
Currently three types of magnetization are discussed in the literature (see
Fig. 1d). Besides the spin and usual orbital magnetization (UOM) induced
respectively by the exchange splitting and spin-orbit interaction (SOI), which
lift spin and orbital degeneracy, the chiral orbital magnetization (COM) was
recently introduced as a quantity emerging independently from SOI due to
magnetic non-coplanarity [11, 12, 13, 14, 15, 16, 17]. Non-vanishing ground-
state currents [18], giving rise to an orbital magnetization, develop
proportionally to the three-spin scalar chirality
$\chi_{ijk}=(\mathbf{S}_{i}\times\mathbf{S}_{j})\cdot\mathbf{S}_{k}$, i.e. the
solid angle defined by a plaquette of three magnetic moments with
$\mathbf{S}_{i}$ being the unit vector of the magnetic moment carried by atom
$i$. $\chi$ induces the emergent magnetic field [19, 20] responsible for
various Hall effects [19, 21, 22], triggers electromagnetic induction [23, 24]
and is distinct from the two-spin chirality vector
$\boldsymbol{{\kappa}}_{ij}=\mathbf{S}_{i}\times\mathbf{S}_{j}$ associated
with the Dzyaloshinskii-Moriya interaction [25, 26]. Ultimately, the total
three-spin scalar chirality summed up over the whole topological spin-texture
with weakly rotating spin moments is quantized by being proportional to the
topological charge, which enforces a topological behavior on key transport
phenomena and promotes the chiral orbital magnetization to a topological one
[12]. We refer to the aforementioned chiral orbital magnetization as the
three-spin chiral orbital magnetization (3-spin COM).

Utilizing fundamental concepts in multiple-scattering theory, we demonstrate
the emergence of new chiral orbital magnetizations triggered by SOI and non-
collinearity (Fig. 1d) via mechanisms involving multi-spin plaquettes. The
first one is directly proportional to the two-spin chirality
$\boldsymbol{\kappa}$, called then 2-spin COM, and therefore changes sign upon
chirality switching, also identified semi-classically [16]. The second one
looks like an interference effect produced by both the three-spin chirality
$\chi$ and SOI and emerges from a four-spin plaquette. It is therefore refered
to as the 4-spin COM. We quantify from first-principles the various
magnetizations carried by a realistic few-nanometers wide Néel-type skyrmion
known to emerge in fcc-PdFe bilayer deposited on Ir(111) surface [27, 28, 29,
30, 31] and explore the magnetization responses to single atomic defects
implanted in the Pd overlayer (see Figure 1a and Methods). We demonstrate that
the impurities can inject large magnetizations into the skyrmion. As defects
we consider 3$d$ (Ti, V, Cr, Mn, Fe, Co, Ni) that can carry significant spin
moments (Figure 1c) and 4$d$ (Zr, Nb, Mo, Tc, Ru, Rh) transition metal atoms
with much weaker spin moments as well as Cu and Ag atoms, which are positioned
close to the core of the skyrmion (Figure 1a-b). We identify common patterns
in the spin and orbital magnetization as function of the atomic number of the
defects, i.e. the impurities electronic filling, which are surprisingly also
executed by the SO-independent chiral orbital magnetization. The unveiled
patterns are however distinct for 3$d$ and 4$d$ impurities and find their
origin in the magnetic interactions dictating the deformation of the
skyrmionic texture. We propose the three-spin scalar chirality as a common
descriptor for the spin, orbital and chiral magnetizations, which can be
probed by the imposed impurities. While the conventional SO-induced orbital
magnetization is expected to be linear with the spin magnetization, the three-
spin scalar chirality is a cubic function of the magnetization, which explains
the unveiled dependencies and defines the generic trends.

Figure 1: Various orbital magnetizations emerging in a magnetic skyrmion. a
Schematic drawing of the cross-section of PdFe/Ir(111) surface, which
incorporates an atomic defect (green sphere) implanted in the Pd overlayer. b
Spin structure of a magnetic skyrmion interacting with a Ti defect. The green
arrow and sphere represent the defect. c Spin moments of the series of 3$d$
and 4$d$ atomic defects when interacting with the ferromagnetic PdFe bilayer.
Positive (negative) values indicate a ferromagnetic (antiferromagnetic)
alignement of the impurity spin moment with respect to the substrate
magnetization. d Schematic representation of various mechanisms affecting the
electronic level splitting from which emerge distinct contributions to the
orbital moment. In a solid, the crystal field quenches the orbital moment,
while spin-orbit coupling lifts the orbital degeneracy. Even without spin-
orbit coupling, an orbital moment (3-spin chiral orbital moment) emerges from
the non-coplanarity of three magnetic moments. With spin-orbit coupling an
additional chiral orbital moment is enabled by the two-spin vector chirality
(2-spin chiral orbital moment). A third chiral orbital moment (not shown) can
be induced by a plaquette of four spins, which can be interpreted as an
interference term depending on spin-orbit interaction and the three-spin
scalar chirality.

## Results

###  0.1 Chirality-induced magnetizations.

The shape of a magnetic skyrmion is modified as soon as a defect is implanted
at its vicinity. This naturally modifies the total magnetic moment carried by
the non-trivial spin-texture. As demonstrated in Supplementary Note 1 and
illustrated in Supplementary Figure 1, we expect new orbital moments induced
by the non-coplanarity and/or chirality of a spin-texture besides the three
types of magnetic moments discussed in the literature: (i) the spin, (ii) the
orbital and (iii) chiral (topological) orbital magnetic moments. Utilizing
multiple-scattering theory, we demonstrate that the orbital moment
$\mathbf{m}^{\text{orb}}$ at a given atomic site $i$ can be approximately
decomposed into four contributions categorized in terms of multi-spin inducing
processes with and without SOI:

$\mathbf{m}^{\mathrm{orb}}=\mathbf{m}^{\mathrm{orb}}_{\mathrm{1-spin}}+\mathbf{m}^{\mathrm{orb}}_{\mathrm{2-spin}}+\mathbf{m}^{\mathrm{orb}}_{\mathrm{3-spin}}+\mathbf{m}^{\mathrm{orb}}_{\mathrm{4-spin}}.$
(1)

We start with the most basic contribution,
$\mathbf{m}^{\mathrm{orb}}_{\mathrm{1-spin}}$, which is the usual orbital
moment (UOM) induced by SOI and is collinear with the atomic spin moment. It
is thus an orbital moment that primarily arises from a single spin term,
$\mathbf{S}_{i}$, in contrast to the two-spin induced orbital moment (2-spin
COM), $\mathbf{m}^{\mathrm{orb}}_{\mathrm{2-spin}}$, which is chiral since it
is proportional to the two-spin chirality $\boldsymbol{\kappa}$, i.e.
$\propto\left(\mathbf{S}_{i}\times\mathbf{S}_{j}\right)$. Owing to the chiral
nature of this term, rotating the moments clockwise or counter-clockwise
modifies the orbital moment in an opposite fashion. It is linear with SOI and
is finite when inversion symmetry is broken. The three-spin induced
$\mathbf{m}^{\mathrm{orb}}_{\mathrm{3-spin}}\propto\chi$ is the chiral orbital
moment (3-spin COM) triggered by the three-spin scalar chirality $\chi$
($\propto(\mathbf{S}_{i}\cdot\left(\mathbf{S}_{j}\times\mathbf{S}_{k}\right)$),
while the four-spin induced orbital moment (4-spin COM),
$\mathbf{m}^{\mathrm{orb}}_{\mathrm{4-spin}}$, is a sort of interference
moment enabled by SOI and proportional to the three-spin scalar chirality
($\mathbf{S}_{i}\chi$). Interestingly, among the four contributing terms, only
$\mathbf{m}^{\mathrm{orb}}_{\mathrm{3-spin}}$ does not request SOI to be
finite. We note that a five-spin chirality term can induce an orbital moment
[15] even without SOI, but it is found negligible in the current study.
Furthermore, for a magnetically collinear state, all chiral orbital
magnetizations cancel out.

###  0.2 Unveiling skyrmion’s magnetization patterns by defect implantation.

Following the method of extraction described in Supplementary Note 1, each of
these orbital moments together with the spin moment summed up over a whole
magnetic skyrmion are presented in Fig. 2, which gives an insight on the
sensitivity of the total skyrmion’s magnetization to the electronic filling of
atomic defects implanted at the vicinity of the skyrmmion’s core. We note that
the skyrmion magnetization is defined as the difference in the total magnetic
moment with and without the skyrmion, which is here of Néel type and has a
diameter of about 2.2n. Overall the net skyrmion magnetization points along
the z-direction.

The 3-spin COM carried by the whole skymion can be dramatically modified by
the presence of atomic defects (see Fig. 2a). Co or Cr impurities can enhance
it by about 6% while a strong reduction can be achieved by Ti. Intriguingly,
one can clearly see that on the one hand 3$d$ atomic defects induce 3-spin
COMs exhibiting a W-shape as function of the impurities atomic number. The
curve shows two minima located at Cr and Co separated by a maximum found at
Fe. On the other hand, 4$d$ defects trigger 3-spin COMs following a V-shape,
with a single minimum at Rh. Whilst the late transition elements (say starting
from Tc) induce moments similar to those of the 3$d$ defects, early elements
(Zr, Nb and Mo) strongly reduce the skyrmion’s chiral orbital moment.

Surprisingly the identified trends, W- and V-shapes for respectively 3$d$ and
4$d$ impurities as schematically illustrated in Fig. 2d, are exactly followed
by the UOM conventionally induced by SOI and the spin magnetization of the
skyrmion (Fig. 2b-c). These results indicate some general pattern behavior
connecting the three types of magnetizations and a common hidden link to the
three-spin scalar chirality. In contrast to the 3-spin COM, the UOM is much
larger and can change by almost 1$\mu_{B}$ across the different impurities
while the spin magnetization exhibits alterations reaching about 25$\mu_{B}$.
This highlights the potential of using atomic defects to engineer the
magnetization carried by a single skyrmion. The newly proposed COM, enabled by
SO and the two-spin chirality (Fig. 2e) or the three-spin chirality (Fig. 2f),
are one order of magnitude smaller than the SO-independent chiral orbital
moment. The shapes of the 2-spin and 4-spin COM for the 4$d$ impurities seem
similar to that obtained for the rest of the magnetization, with a minimum
observed for Mo (2-spin COM) and Pd (4-spin COM). The 2-spin COM induced by
the 3$d$ atomic defects exhibits the W-shape, with a minimum obtained for Cr,
while the second minimum has been pushed to the limit of the series of
transition elements (Cu). For the 4-spin COM, the first minimum induced by the
early 3$d$ transition elements in the rest of magnetizations disappears
leading to a V-like shape with a single minimum for Ni. As aforementioned, the
4-spin COM is a sort of interference term, which is in the particularly
investigated skyrmion weak in comparison to the other studied orbital
magnetizations. In the following, we focus our discussion on the spin
magnetization, the UOM and 3-spin COM, which are significantly larger in
magnitude than the 2-spin and 4-spin COM. The two-dimensional maps of the spin
and orbital moments across the skyrmion with and without impurities are
illustrated in Supplementary Figures 3-8.

Figure 2: Universal trends of distinct skyrmion-carried magnetizations
unveiled by implanted atomic defects. a Total 3-spin chiral orbital
magnetization (COM), b usual orbital magnetization (UOM), c spin
magnetization, e 2-spin COM and f 4-spin COM. d Schematic illustration of the
magnetizations, chiral or not, pertaining to a single skyrmion as function of
the electronic filling across the transition metal series.

###  0.3 Origin of the spin-magnetization patterns.

To understand the origin of the magnetization patterns and the dissimilarity
between those pertaining to 3$d$ versus 4$d$ impurities, we explore how the
latter impact the spin magnetization of the underlying Fe layer, which hosts
the non-trivial spin-textures. Two mechanisms are at play: (i) the competition
between the defect-substrate magnetic interaction with the intralayer Fe-Fe
magnetic interaction, and (ii) the alteration of the latter by the presence of
impurities. In contrast to 4$d$ impurities, 3$d$ atomic defects carry a large
spin magnetic moment (see Fig. 1c) and impose a strong magnetic exchange
interaction to the substrate. As shown in Supplementary Note 2 for a simple
model, the z-magnetization is expected to be proportional to the ratio between
the impurity-substrate magnetic interaction and the nearest neighboring
interaction within the Fe layer, $|\mathrm{J}_{\mathrm{Fe-
imp}}/\mathrm{J}_{\mathrm{Fe-Fe}}|$ [32]. This is calculated for the
investigated skyrmion and plotted in Fig. 3a, which reveals the aforementioned
W-shape. On the same figure, we deliberately show what the Heisenberg model
predicts for the 4$d$ impurities, although the model is not appropriate since
the magnetic moments are much smaller than the 3$d$ ones and are more
sensitive to both the surrounding and spin-rotation. The latter ingredients
prevent the use of the Heisenberg model, which interestingly predicts a
behavior similar to that induced by the 3$d$ atomic defects.

The origin of the V-shape found for the 4$d$ impurities lies in the way the
effective magnetic interaction among the substrate Fe atoms is affected by the
presence of the defects. In the defect-free case, the overlayer Pd atoms
enhance the nearest neighboring Fe-Fe magnetic interaction in a linear fashion
[28]. Indeed the Pd atoms are characterized by a ferromagnetic Stoner
susceptibility, which at the vicinity of the magnetic substrate carry an
induced moment facilitated by the Stoner criterion for magnetism. The bare
nearest neighboring Fe-Fe magnetic interaction is then enhanced by 16% in
average [28]. A similar enhancement is expected for impurities that tend to
couple ferromagnetically to the substrate, which is the case for example of Rh
(see Fig. 3d). For instance, when replacing the Pd atom by the Rh, the
magnetic moment locally increases from $\qty{0.3}\mu_{\text{B}}$ to
$\qty{0.9}\mu_{\text{B}}$. This stiffens the skyrmion and enlarges its core as
it can be grasped from Fig. 3c, which leads to a negative contribution to the
z-component of the magnetization and reduces it. Mo impurity, however, couples
antiferromagnetically to the Fe substrate, this affects the surrounding Pd
atoms by reducing their magnetization as shown in Fig. 3e, which can be
compared to the Rh-case in Fig. 3f. This diminishes the magnitude of the
effective Fe-Fe magnetic interaction in comparison to the defect-free scenario
as illustrated in Fig. 3d. In this case, the outer-region of the skyrmion,
contributing positively to the z-magnetization, wins even more over the core
region (Fig. 3b) and the total magnetization is enhanced. For completeness,
two-dimensional maps of the spin moments of Pd overlayer atoms are presented
in Supplementary Figures 9-11.

Figure 3: Impact of substrate-defect interactions on the trends of skyrmion’s
magnetization. a Ratio between the next neighboring impurity-substrate and Fe-
Fe magnetic exchange interactions for the case of 3$d$ and 4$d$ atomic defects
as a function of the impurities atomic number. Change of the substrate Fe
magnetic moments due to the presence of a b Mo and c Rh defects when comparing
with a defect-free skyrmion. d Local change of the first nearest neighbour
exchange interaction within the Fe-layer due to Mo- and Rh-defects. The red
(blue) color of the bond between the atoms indicates the exchange interaction
decreased (increased) as compared to the defect-free substrate. Modification
of the overlayer Pd atoms induced by e Mo and f Rh defects as obtained in the
collinear state.

###  0.4 Correlation between spin, orbital and chiral magnetizations of
skyrmions.

The previous subsection deals with the origin of the patterns characterizing
the spin magnetization carried by a skyrmion. Here, we explain how the
unveiled trends are also obeyed by both the UOM and 3-spin COM. To get first
insight, we present in Fig. 4 various correlation plots systematically
collected for the different defects implanted close to the core of the
magnetic skyrmion. For instance, Fig. 4a shows that the 3-spin COM carried by
the skyrmion is linear with the total three-spin scalar chirality, while Fig.
4b recovers the expected linear dependence between the UOM and the spin
magnetization. An intriguing cubic dependence is, however, found between the
three-spin chiral orbital and spin magnetizations (Fig. 4c), which can be
fitted by the following equation

$\text{3-spin
COM}=\mathrm{A}\,\mathrm{M}_{\text{sky}}^{\text{spin}}+\mathrm{B}\,{\mathrm{M}_{\text{sky}}^{\text{spin}}}^{2}+\mathrm{C}\,{\mathrm{M}_{\text{sky}}^{\text{spin}}}^{3},$
(2)

where $\mathrm{A}=-1.6410^{-2}$, $\mathrm{B}=-1.3910^{-3}\mu_{B}^{-1}$ and
$\mathrm{C}=1.8710^{-5}\mu_{B}^{-2}$ are material-dependent parameters. We
show analytically in Supplementary Note 3 that
$\mathrm{A}=\frac{\sqrt{3}}{2}$, $\mathrm{B}=0$ and
$\mathrm{C}=-\frac{\sqrt{3}}{18}$ for a non-collinear trimer in a Néel state
having initially its three magnetic moments lying in-plane, which
simultaneously interact with an atomic magnetic defect. A more involved
analytical evaluation comfort the cubic dependence presented in Eq. 2 and
conjecture therefore that it settles a universal correlation between the total
3-spin COM and the spin magnetization as well as the UOM.

While the trends dictating the spin-magnetization are induced by the
modification of magnetic interactions across the transition atomic defect via
electronic filling of the impurities, it is the SOI and the cubic dependence
of the three-spin scalar chirality that enforce the same trends on
respectively the conventional orbital and chiral magnetization.

Figure 4: Universal correlations between skyrmion’s magnetizations and device
concept. a The skyrmion’s 3-spin COM as function of the total 3-spin scalar
chirality collected from every magnetic profile induced by defects of the 3$d$
and 4$d$ series. b The skyrmion’s UOM and c 3-spin COM as function of the
skyrmion’s spin magnetization. While the 3-spin COM and UOM exhibit a linear
behavior, as expected, with respectively the 3-spin scalar chirality and spin
magnetization, a cubic function describes the relation between the 3-spin COM
and spin magnetization. d Based on the universal dependence of the skyrmion’s
magnetization with respect to the atomic number of an impurity, a device is
envisaged with a controlled implantation of defects to engineer the amplitude
of the magnetization. The correlation among the different types of skyrmion’s
magnetizations enables to monitor, design and quantify the chiral orbital
contributions once the spin magnetization is measured.

## Discussion

In this study, we delve into the intricate realm of the magnetization carried
by skyrmions. Hinging on multiple-scattering mechanisms, we unveil chiral
orbital magnetizations, sparked by multi-spin plaquettes intertwined with
spin-orbit interaction, non-collinearity, two-spin and three-spin chiralities.
Employing rigorous first-principles quantification, we explore the magnetic
landscape of a realistic Néel-type skyrmion residing in an fcc-PdFe bilayer
atop an Ir(111) surface. Our inquiry extends to the realm of single atomic
defects in the Pd overlayer, revealing the remarkable ability of impurities,
to infuse substantial magnetizations into the skyrmion.

We discern compelling patterns universally followed by the spin and orbital
magnetizations (usual and chiral) intricately tied to the atomic
characteristics of defects and electronic filling. These patterns showcase
nuances between 3$d$ and 4$d$ impurities, stemming from magnetic interactions
that sculpt the deformation of the skyrmionic texture. Dissecting the
correlations between the different skyrmion’s magnetizations influenced by the
implanted defects reveals a cubic relation between the spin and chiral orbital
magnetization instead of the linear dependence characterizing the spin and
orbital magnetizations.

By elucidating the universal dependencies governing skyrmion magnetization, we
contribute to the foundational principles that govern magnetism, paving the
way for establishing overarching trends in the intricate domain of chiral
magnetism. Our work anticipates the prospect of quantifying chiral
magnetization, an aspect yet unmeasured, by precisely determining spin
magnetization, strategically leveraging defects that amplify its magnitude.

We envision the transformative potential of strategically implanting defects
into surfaces to engineer and manipulate the magnetization carried by
skyrmions. Notably, our findings indicate that 3$d$ impurities yield to a more
important amplifying effect of the skyrmion’s magnetization compared to their
4$d$ counterparts, presenting a promising avenue for advancing skyrmionic bits
in the development of disruptive storage devices.

###  0.5 Computational details.

The simulations are based on density functional theory within the local spin
density approximation as implemented in the all-electron full-potential
scalar-relativistic Korringa-Kohn-Rostoker (KKR) Green function method with
spin-orbit interaction added self-consistently. [33, 34] The method is ideal
to treat within an embedding scheme a single confined magnetic skyrmion
interacting with atomic defects all embedded on a perfect magnetic substrate.
The skyrmion is hosted by the Fe layer and lives within 37 Fe atoms. Nearest
neighbors are included, which lead to an embedding cluster of 124 atoms. The
spin-texture characterizing the skyrmion is obtained in a self-consistent
fashion. The Pd/Fe/Ir(111) slab consists of a fcc-stacked PdFe bilayer
deposited on 34 layers of Ir with atomic positions obtained from ab initio. We
assume an angular momentum cutoff at lmax=3 for the orbital expansion of the
Green function and the energy contour contained 42 grid points in the upper
complex plane with seven Matsubara poles and a Brilloun zone mesh of
30$\times$30 k-points for the self-consistent description of skyrmion-free
substrate properties (see more technical details in Refs. [28, 29, 35, 36,
37]).

###  0.6 Analytical derivation of the new chiral magnetizations.

The concept used and presented in Supplementary Note 1 to establish the
existence of orbital magnetizations emerging from multiple-scattering events
bears some similarities to that utilized to unveil new chiral multi-site
magnetic interactions [38, 39, 40, 41], hyperfine interaction [42], and
magnetoresistances [31].

###  0.7 Data availability

The data that support the findings of this study are present in the paper and
the Supplementary Information.

###  0.8 Code availability

The KKR code is a rather complex ab-initio DFT-based code, which is in general
impossible to use without proper training on the theory behind it and on the
practical utilization of the code. We are happy to provide the latter code
upon request.

We thank Moritz Wineterott for fruitful discussions. This work is supported by
the Deutsche Forschungsgemeinschaft (DFG) through SPP 2137 “Skyrmionics”
(Project LO 1659/8-1). We acknowledge the computing time granted by the JARA-
HPC Vergabegremium and VSR commission on the supercomputers JURECA at
Forschungszentrum Jülich and at the supercomputing centre of RWTH Aachen
University.

S.L. initiated, designed, and supervised the project. I.L.F. performed the
first-principles calculations and post-processed the data. S.L. provided the
analytical derivations of the various orbital magnetizations. I.L.F. and S.L.
discussed the results and S.L. wrote the manuscript.

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to S.L. (email:
s.lounis@fz-juelich.de).

## References

## References

  * [1] Slater, J. C.  The ferromagnetism of nickel. ii. temperature effects.  _Phys. Rev._ 49, 931–937 (1936).  URL https://link.aps.org/doi/10.1103/PhysRev.49.931.
  * [2] Pauling, L.  The nature of the interatomic forces in metals.  _Phys. Rev._ 54, 899–904 (1938).  URL https://link.aps.org/doi/10.1103/PhysRev.54.899.
  * [3] Slater, J. C.  Cohesion in monovalent metals.  _Phys. Rev._ 35, 509–529 (1930).  URL https://link.aps.org/doi/10.1103/PhysRev.35.509.
  * [4] Slater, J. C.  Atomic shielding constants.  _Phys. Rev._ 36, 57–64 (1930).  URL https://link.aps.org/doi/10.1103/PhysRev.36.57.
  * [5] Sommerfeld, A. & Bethe, H.  _Elektronentheorie der Metalle_ , 333–622 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1933).
  * [6] Bogdanov, A. N. & Yablonskii, D.  Thermodynamically stable “vortices” in magnetically ordered crystals. the mixed state of magnets.  _J. Exp. Theor. Phys._ 95, 178 (1989).
  * [7] Rössler, U. K., Bogdanov, A. N. & Pfleiderer, C.  Spontaneous skyrmion ground states in magnetic metals.  _Nature_ 442, 797–801 (2006).
  * [8] Nagaosa, N. & Tokura, Y.  Topological properties and dynamics of magnetic skyrmions.  _Nat. Nanotech._ 8, 899–911 (2013).
  * [9] Fert, A., Cros, V. & Sampaio, J.  Skyrmions on the track.  _Nat. Nanotech._ 8, 152–156 (2013).
  * [10] Fert, A., Reyren, N. & Cros, V.  Magnetic skyrmions: advances in physics and potential applications.  _Nat. Rev. Mater._ 2, 17031 (2017).
  * [11] Shindou, R. & Nagaosa, N.  Orbital ferromagnetism and anomalous hall effect in antiferromagnets on the distorted fcc lattice.  _Phys. Rev. Lett._ 87, 116801 (2001).  URL https://link.aps.org/doi/10.1103/PhysRevLett.87.116801.
  * [12] dos Santos Dias, M., Bouaziz, J., Bouhassoune, M., Blügel, S. & Lounis, S.  Chirality-driven orbital magnetic moments as a new probe for topological magnetic structures.  _Nature Communications_ 7, 13613 (2016).
  * [13] Hanke, J.-P. _et al._ Role of berry phase theory for describing orbital magnetism: From magnetic heterostructures to topological orbital ferromagnets.  _Phys. Rev. B_ 94, 121114 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.94.121114.
  * [14] dos Santos Dias, M. & Lounis, S.  Insights into the orbital magnetism of noncollinear magnetic systems.  _Spintronics X_ 10357, 103572A (2017).  URL https://doi.org/10.1117/12.2275305.
  * [15] Bouaziz, J., Dias, M. d. S., Guimarães, F. S. M., Blügel, S. & Lounis, S.  Impurity-induced orbital magnetization in a rashba electron gas.  _Phys. Rev. B_ 98, 125420 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.125420.
  * [16] Lux, F. R., Freimuth, F., Blügel, S. & Mokrousov, Y.  Engineering chiral and topological orbital magnetism of domain walls and skyrmions.  _Communication Physics_ 1, 60 (2018).
  * [17] Brinker, S., dos Santos Dias, M. & Lounis, S.  Spin, atomic, and interatomic orbital magnetism induced by $3d$ nanostructures deposited on transition metal surfaces.  _Phys. Rev. Mater._ 4, 024404 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevMaterials.4.024404.
  * [18] Tatara, G. & Garcia, N.  Quantum toys for quantum computing: Persistent currents controlled by the spin josephson effect.  _Phys. Rev. Lett._ 91, 076806 (2003).  URL https://link.aps.org/doi/10.1103/PhysRevLett.91.076806.
  * [19] Bruno, P., Dugaev, V. K. & Taillefumier, M.  Topological hall effect and berry phase in magnetic nanostructures.  _Phys. Rev. Lett._ 93, 096806 (2004).  URL https://link.aps.org/doi/10.1103/PhysRevLett.93.096806.
  * [20] Everschor-Sitte, K. & Sitte, M.  Real-space berry phases: Skyrmion soccer (invited).  _Journal of Applied Physics_ 115, 172602 (2014).  URL https://doi.org/10.1063/1.4870695.
  * [21] Lux, F. R., Freimuth, F., Blügel, S. & Mokrousov, Y.  Chiral hall effect in noncollinear magnets from a cyclic cohomology approach.  _Phys. Rev. Lett._ 124, 096602 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.124.096602.
  * [22] Bouaziz, J., Ishida, H., Lounis, S. & Blügel, S.  Transverse transport in two-dimensional relativistic systems with nontrivial spin textures.  _Phys. Rev. Lett._ 126, 147203 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.126.147203.
  * [23] Nagaosa, N.  Emergent inductor by spiral magnets.  _Japanese Journal of Applied Physics_ 58, 120909 (2019).
  * [24] Yokouchi, T. _et al._ Emergent electromagnetic induction in a helical-spin magnet.  _Nature_ 586, 232 (2020).
  * [25] Dzyaloshinsky, I.  A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics.  _J. Phys. Chem. Sol._ 4, 241 – 255 (1958).
  * [26] Moriya, T.  Anisotropic superexchange interaction and weak ferromagnetism.  _Phys. Rev._ 120, 91–98 (1960).
  * [27] Romming, N. _et al._ Writing and deleting single magnetic skyrmions.  _Science_ 341, 636–639 (2013).
  * [28] Lima Fernandes, I., Bouaziz, J., Blügel, S. & Lounis, S.  Universality of defect-skyrmion interaction profiles.  _Nat. Commun._ 9, 4395 (2018).
  * [29] Fernandes, I. L., Bouhassoune, M. & Lounis, S.  Defect-implantation for the all-electrical detection of non-collinear spin-textures.  _Nature Commun._ 11, 1–9 (2020).
  * [30] Fernandes, I. L., Chico, J. & Lounis, S.  Impurity-dependent gyrotropic motion, deflection and pinning of current-driven ultrasmall skyrmions in PdFe/Ir(111) surface.  _Journal of Physics: Condensed Matter_ (2020).  URL http://iopscience.iop.org/10.1088/1361-648X/ab9cf0.
  * [31] Fernandes, I. L., Blügel, S. & Lounis, S.  Spin-orbit enabled all-electrical readout of chiral spin-textures.  _Nature Communications_ 13, 1576 (2022).  URL https://doi.org/10.1038/s41467-022-29237-0.
  * [32] Lounis, S., Mavropoulos, P., Zeller, R., Dederichs, P. H. & Blügel, S.  Noncollinear magnetism of cr and mn nanoclusters on ni(111): Changing the magnetic configuration atom by atom.  _Phys. Rev. B_ 75, 174436 (2007).  URL https://link.aps.org/doi/10.1103/PhysRevB.75.174436.
  * [33] Bauer, D. S. G.  Development of a relativistic full-potential first-principles multiple scattering green function method applied to complex magnetic textures of nano structures at surfaces.  _PhD dissertation at the RWTH-Aachen_ (2013).
  * [34] Crum, D. M. _et al._ Perpendicular reading of single confined magnetic skyrmions.  _Nat. Commun._ 6, 8541 (2015).
  * [35] Arjana, I. G., Lima Fernandes, I., Chico, J. & Lounis, S.  Sub-nanoscale atom-by-atom crafting of skyrmion-defect interaction profiles.  _Scientific Reports_ 10, 14655 (2020).  URL https://doi.org/10.1038/s41598-020-71232-2.
  * [36] Bouhassoune, M. & Lounis, S.  Friedel oscillations induced by magnetic skyrmions: From scattering properties to all-electrical detection.  _Nanomaterials_ 11 (2021).  URL https://www.mdpi.com/2079-4991/11/1/194.
  * [37] Aldarawsheh, A. _et al._ Emergence of zero-field non-synthetic single and interchained antiferromagnetic skyrmions in thin films.  _Nature Communications_ 13, 7369 (2022).  URL https://doi.org/10.1038/s41467-022-35102-x.
  * [38] Brinker, S., Dias, M. d. S. & Lounis, S.  The chiral biquadratic pair interaction.  _New Journal of Physics_ 21, 083015 (2019).
  * [39] Lounis, S.  Multiple-scattering approach for multi-spin chiral magnetic interactions: application to the one- and two-dimensional rashba electron gas.  _New Journal of Physics_ 22, 103003 (2020).  URL https://dx.doi.org/10.1088/1367-2630/abb514.
  * [40] Grytsiuk, S. _et al._ Topological chiral magnetic interactions driven by emergent orbital magnetism.  _Nature Communications_ 11, 1–7 (2020).
  * [41] Brinker, S., dos Santos Dias, M. & Lounis, S.  Prospecting chiral multisite interactions in prototypical magnetic systems.  _Phys. Rev. Research_ 2, 033240 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033240.
  * [42] Shehada, S., dos Santos Dias, M., Abusaa, M. & Lounis, S.  Interplay of magnetic states and hyperfine fields of iron dimers on mgo (001).  _Journal of Physics: Condensed Matter_ 34, 385802 (2022).  URL https://dx.doi.org/10.1088/1361-648X/ac8135.