# Magnetic ground state of supported monatomic Fe chains from first principles

B Nagyfalusi1,2, L Udvardi2,3 and L Szunyogh2,3 1 Wigner Research Centre for
Physics, Institute for Solid State Physics and Optics, H-1525 Budapest,
Hungary  2 Department of Theoretical Physics, Institute of Physics, Budapest
University of Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest,
Hungary  3 MTA-BME Condensed Matter Research Group, Budapest University of
Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary
<EMAIL_ADDRESS>
###### Abstract

A new computational scheme is presented based on a combination of the
conjugate gradient and the Newton-Raphson method to self-consistently minimize
the energy within local spin-density functional theory, thus to identify the
ground state magnetic order of a finite cluster of atoms. The applicability of
the new ab initio optimization method is demonstrated for Fe chains deposited
on different metallic substrates. The optimized magnetic ground states of the
Fe chains on Rh(111) are analyzed in details and a good comparison is found
with those obtained from an extended Heisenberg model containing first
principles based interaction parameters. Moreover, the effect of the different
bilinear spin-spin interactions in the formation of the magnetic ground states
is monitored. In case of Fe chains on Nb(110) spin-spiral configurations with
opposite rotational sense are found as compared to previous spin-model results
which hints on the importance of higher order chiral interactions. The
wavelength of the spin-spiral states of Fe chains on Re(0001) was obtained in
good agreement with scanning tunneling microscopy experiments.

††: J. Phys.: Condens. Matter

Keywords: magnetic nanoclusters, ab initio, magnetic ground state

##  1 Introduction

The rapid development of scanning tunneling microscopy (STM) allows to create
finite clusters of magnetic atoms on solid surfaces with atomic precision [1,
2]. Magnetic chains on top of a superconducting surface has attracted
considerable interest during the past years due to the challenge of finding
Majorana bound states that may open a new prospect in quantum computation [3,
4, 5, 6, 7]. The magnetic structure of the chains has crucial role in the
emergence of topological superconductivity. A recent theoretical and
experimental review on the physics of atomic spin chains deposited on surfaces
is provided in [8].

Complex magnetic systems are often described in terms of a generalized
Heisenberg model, where the interaction between spins $\mathbf{s}_{i}$ and
$\mathbf{s}_{j}$ at sites $i$ and $j$ is ascribed to a tensorial exchange
coupling $\mathbf{J}_{ij}$. This exchange coupling can further be split into
three terms [9]: the isotropic exchange, the symmetric anisotropic exchange
and the antisymmetric anisotropic exchange or the Dzyaloshinky-Moriya
interaction (DMI) [10, 11]. Of particular interest is the DMI, which can lead
to the formation of spin-spiral ground states [12, 13, 14, 15] or, in the
presence of external magnetic field, to the emergence of magnetic skyrmions in
chiral bulk magnets [16, 17, 18] and in ultrathin magnetic films [19, 20, 21,
22]. Although the bilinear Heisenberg model provides a sufficient description
of a large number of magnetic materials, it has been reported that higher
order spin-spin interactions are necessary to account for the complex magnetic
ordering in several systems [23, 24, 25, 26, 27].

In order to avoid using spin models the magnetic ground state of
nanostructures can directly be determined from first principles calculations
based on density functional theory. The real-space linear muffin-tin orbital
method in the atomic sphere approximation (RS-LMTO-ASA) extended to non-
collinear magnetic systems has been succesfully used to obtain the magnetic
configurations of supported nanoclusters[28, 29, 30]. Including the effects of
spin-orbit coupling, this parameter-free method was applied to explore the
complex magnetic ground state of Mn nanowires on Ag(111) and Au(111) [31], as
well as of Cr nano-islands on Pd(111) [32]. Using the constrained local moment
method employed with the multiple scattering Green’s function technique[33,
34] the existence of a canted magnetic state of a Co chain along a Pt(111)
surface step edge [35] was demonstrated from first principles in good
agreement with experiment [36]. Deriving explicit expressions for the torque
acting on the magnetic moments, the formation of a domain wall through a nano-
contact has been studied within an ab-initio framework [37]. Using the same
approach for the torque and by solving the stochastic Landau-Lifshitz-Gilbert
equations, it became possible to perform finite temperature spin-dynamics
simulations for monatomic Co chains on top of Au(001) [38].

In this work a modified version of the conjugate gradient and Newton-Raphson
methods used in [37] is proposed to find the magnetic ground state of
deposited clusters from first principles. First the numerical method will be
described in details, then it will be applied to Fe chains deposited on the
(111) surface of face centered (fcc) Rh. The magnetic ground states derived
from the new ab-initio optimization technique are compared to those obtained
from spin models parametrized from first principles and the effect of the
different kind of spin-spin interactions is discussed. Further applications
are presented for Fe chains deposited on the (110) surface of body-centered
cubic (bcc) Nb and on the (0001) surface of hexagonal closed pack (hcp) Re,
where spin-spiral ground states are found in agreement with previous
theoretical works[25, 39] and experiment[4].

##  2 Methods

###  2.1 Embedded cluster calculations

The electronic structure of the host system, composed of a semi infinite
substrate, a finite number of substrate and empty sphere layers describing the
surface region and a semi infinite vacuum, has been calculated in the terms of
the fully relativistic screened Korringa-Kohn-Rostoker (SKKR) method[40]. Then
the Green’s function embedding technique based on multiple scattering theory
[41] has been applied to determine the electronic and magnetic properties of
the finite atomic clusters. The charge and spin density have been relaxed not
only for the magnetic atoms forming the chain but also in the first shell of
atomic cells around them in the substrate and in the surrounding vacuum. For
the construction of the effective potentials, $V_{\mathrm{eff},i}(\mathbf{r})$
and exchange-correlation fields, $\mathbf{B}_{\mathrm{xc},i}(\mathbf{r})$,
where $i$ labels atomic positions, the atomic sphere approximation (ASA) has
been applied. The exchange-correlation field in the atomic sphere
corresponding to site i is then given by
$\mathbf{B}_{\mathrm{xc},i}(\mathbf{r})=B_{\mathrm{xc},i}(r)\,{\bf s}_{i}$,
where ${\bf s}_{i}$ is a unit vector representing the direction of
$\mathbf{B}_{\mathrm{xc},i}$. An angular momentum cut-off $\ell_{\rm max}=3$
was used for the partial waves in multiple scattering theory. The details of
the geometry used in the calculations will be specified for each system under
consideration.

###  2.2 The modified conjugate gradient method

In order to find the magnetic configuration of a cluster of atoms, a
simultaneous iteration of the direction of the magnetization, of the effective
potential and of the exchange-correlation field works flawlessly in case of a
simple systems. However, in case of complex magnetic structures, e.g. for
frustrated antiferromagnetic or for spin-spiral ground states, the convergence
of the above scheme might be less satisfactory. An estimate of the first and
second order change of the energy of the system with respect to the change of
the orientation of the magnetization can help to find the magnetic ground
state configuration.

In the spirit the magnetic force theorem [42], $V_{\mathrm{eff},i}(r)$ and
$B_{\mathrm{xc},i}(r)$ are kept fixed for small changes of the magnetic
orientations. This implies that the number of the electrons in the finite
cluster is not preserved, therefore, the variation of the energy of the system
at zero temperature is replaced by the variation of the grand potential,
$\Omega=E-\varepsilon_{\mathrm{F}}N$, where the energy $E$ of the system is
approximated by the single particle (band) energy,

$\Delta\Omega=\int\limits_{-\infty}^{\varepsilon_{\mathrm{F}}}\Delta
n(\varepsilon)(\varepsilon-\varepsilon_{\rm
F})\mathrm{d}\varepsilon=-\int\limits_{-\infty}^{\varepsilon_{\mathrm{F}}}\Delta
N(\varepsilon)\mathrm{d}\varepsilon\;,$ (1)

with $\varepsilon_{\mathrm{F}}$ and
$N(\varepsilon)=\int_{-\infty}^{\varepsilon}n(\varepsilon^{\prime})\mathrm{d}\varepsilon^{\prime}$
being the Fermi energy and the integrated density of states (DOS),
respectively. The main advantage of the above formula is that within the
multiple scattering theory it can easily be calculated using Lloyd’s formula
[43, 9].

It is straightforward to derive a formula for the first order change of the
grand potential when the exchange correlation field in an atomic sphere at
site i is rotated around a unit vector $\mathbf{n}_{i}$ with an angle
$\Delta\varphi_{i}$[9]:

$\Delta\Omega^{(1)}_{i}=\mathbf{T}_{i}\Delta{\varphi}_{i}\,,$ (2)

where the torque $\mathbf{T}_{i}$ is defined as

$\mathbf{T}_{i}=\frac{1}{\pi}\int\limits_{-\infty}^{\varepsilon_{\mathrm{F}}}\mathrm{Im\,Tr}\left(\frac{i}{\hbar}\left[\mathbf{J},t^{-1}_{i}\right]\tau_{ii}\right)\mathrm{d}\varepsilon\,,$
(3)

and $\Delta{\varphi}_{i}\equiv\mathbf{n}_{i}\Delta\varphi_{i}$,
$\mathbf{J}=\mathbf{L}+\mathbf{S}$ is the matrix of the total angular momentum
operator, $t_{i}$ is the single-site scattering matrix and $\tau_{ii}$ is the
site-diagonal matrix block of the scattering path operator (see [41]). Note
that these matrices are defined in orbital momentum and spin space.

First, $V_{\mathrm{eff},i}$ and $B_{\mathrm{xc},i}$ are kept fixed and we
search for the directions of the exchange-correlation field ${\bf s}_{i}$ at
each site corresponding to the lowest grand potential. Once the restricted
minimum of $\Omega$ is found, the set of ${\bf s}_{i}$ is kept frozen, while
$V_{\mathrm{eff},i}$ and $B_{\mathrm{xc},i}$ are iterated according to the
local spin-density approximation. The whole procedure is repeated until the
error of the effective potential and the exchange correlation field as well as
the torque on each site become smaller than a predefined value.

During the optimization ${\bf s}_{i}$ is moving on the surface of a unit
sphere. The minimum of the grand potential should be found on the direct
product of $N$ unit spheres, where $N$ is the number of atoms included in the
optimization process, thus the standard optimization methods working in
Euclidean space doesn’t apply. An extension of the conjugate gradient
algorithm to Riemann manifolds has been used [44, 45, 46] as briefly
summarized in the following. Let $\mathbf{t}_{i}^{(k)}$ be the torque vector
which specifies the rotation for site i at the k-th iteration step. At the
beginning of the iteration procedure $\mathbf{t}_{i}^{(0)}$ is set equal to
$\mathbf{T}_{i}$ given by (3). In the first part of an iteration step we look
for the minimum of the grand potential with respect to $\alpha$,

$\Omega(\alpha)=\Omega\left\\{O\left(\mathbf{n}_{i}^{(k)},\alpha
t_{i}^{(k)}\right){\bf s}_{i}\right\\}\;,$ (4)

where $t_{i}^{(k)}=|\mathbf{t}_{i}^{(k)}|$,
$\mathbf{n}_{i}^{(k)}=\mathbf{t}_{i}^{(k)}/{t}_{i}^{(k)}$ and
$O(\mathbf{n},\varphi)$ denotes a rotation around vector $\mathbf{n}$ by an
angle of $\varphi$. The vector $\mathbf{t}_{i}^{(k)}$ plays the same role as
the search direction in a conventional conjugate gradient method. Once the
minimum with respect to $\alpha$ has been found, we calculate the torque using
(3) in the new configuration and define the effective torque for the next
iteration step,

$\mathbf{t}_{i}^{(k+1)}=\mathbf{T}_{i}^{(k+1)}+\beta_{i}\mathbf{t}_{i}^{(k)}\,,$
(5)

with the multipliers $\beta_{i}$ from the Polak-Ribière formula[47],

$\beta_{i}=\frac{\left(\mathbf{T}_{i}^{(k+1)}-\mathbf{T}_{i}^{(k)}\right)\mathbf{T}_{i}^{(k+1)}}{\mathbf{T}_{i}^{(k+1)}\mathbf{T}_{i}^{(k)}}\,.$
(6)

The iterations are repeated until satisfactory convergence is achieved.

In order to further increase the efficiency of the optimizing procedure, we
applied a Newton-Raphson method adopted to the Riemann manifold. The second
order change of the grand potential necessary for the Newton-Raphson method
can also be calculated using the Lloyd’s formula [9],

$\Delta\Omega^{(2)}_{ij}={\frac{1}{2}}\Delta{\varphi}_{i}H_{ij}\Delta{\varphi}_{j}\,,$
(7)

with the $3\times 3$ matrix,

$\displaystyle
H_{ij}=\frac{1}{\pi}\int\limits_{-\infty}^{\varepsilon_{\mathrm{F}}}\frac{1}{\hbar^{2}}\mathrm{Im\,Tr}\left(\left[\mathbf{J},t^{-1}_{i}\right]\tau_{ij}\left[\mathbf{J},t^{-1}_{j}\right]\tau_{ji}\right.$
$\displaystyle\left.-\delta_{ij}\left\\{\left[\mathbf{J},\left[\mathbf{J},t^{-1}_{i}\right]\right]\right\\}\tau_{ii}\right)\mathrm{d}\varepsilon\,.$
(8)

The change of the grand potential is then given up to second order by

$\displaystyle\Delta\Omega$
$\displaystyle=\sum_{i}\Delta\Omega^{(1)}_{i}+\sum_{ij}\Delta\Omega^{(2)}_{ij}=$
(9)
$\displaystyle=\sum_{i}\mathbf{T}_{i}\Delta{\varphi}_{i}+\frac{1}{2}\sum_{ij}\Delta{\varphi}_{i}H_{ij}\Delta{\varphi}_{j}\;.$

The component of the total angular momentum operator parallel to ${\bf s}_{i}$
commutes with the single site scattering matrix $t_{i}$, $[{\bf
s}_{i}\mathbf{J},t_{i}]=0$. Consequently the following relationships apply,

$\displaystyle\mathbf{T}_{i}{\bf s}_{i}=0,$ (10) $\displaystyle H_{ij}{\bf
s}_{j}={\bf 0}\,.$ (11)

This implies that a set of rotation
$\displaystyle\left\\{\Delta{\varphi}_{i}\right\\}$ lowering the grand
potential must be constrained to the manifold perpendicular to
$\displaystyle\left\\{{\bf s}_{i}\right\\}$. For that reason two auxiliary
unit vectors per site, $\mathbf{e}_{1i}$ and $\mathbf{e}_{2i}$, are introduced
which together with ${\bf s}_{i}$ form a right-handed orthonormal system as
depicted in figure 1 and the rotations are chosen as

$\Delta{\varphi}_{i}=\Delta\varphi_{1i}\mathbf{e}_{1i}+\Delta\varphi_{2i}\mathbf{e}_{2i}\,.$
(12)

Figure 1:  Rotation of the magnetization vector ${\bf s}_{i}$ around the
orthogonal directions $\mathbf{e}_{1i}$ and $\mathbf{e}_{2i}$ by angles
$\Delta\varphi_{1i}$ and $\Delta\varphi_{2i}$, respectively. The new
magnetization direction is denoted by ${\bf s}_{i}^{\prime}$.

Varying the change of the grand potential (9) with respect to
$\Delta\varphi_{\gamma i}$ $(\gamma=1,2)$, we obtain the Newton equations,

$\displaystyle T_{\gamma i}+\sum_{\gamma^{\prime}j}\mathcal{H}_{\gamma
i,\gamma^{\prime}j}\Delta\varphi_{\gamma^{\prime}j}\equiv$
$\displaystyle\mathbf{e}_{\gamma
i}\mathbf{T}_{i}+\sum_{\gamma^{\prime}j}\mathbf{e}_{\gamma
i}H_{ij}\mathbf{e}_{\gamma^{\prime}j}\Delta\varphi_{\gamma^{\prime}j}=0\,.$
(13)

Inverting the Hessian, $\mathcal{H}=\left\\{\mathcal{H}_{\gamma
i,\gamma^{\prime}j}\right\\}$, a new set of $\Delta\varphi_{\gamma i}$ is
calculated,

$\Delta\varphi_{\gamma i}=-\sum_{\gamma^{\prime}j}\mathcal{H}^{-1}_{\gamma
i,\gamma^{\prime}j}T_{\gamma^{\prime}j}\,,$ (14)

thus a new spin-configuration is generated. The Newton-Raphson method
converges to the local minimum if the $\mathcal{H}$ is positive definite. In
our implementation, at the end of every line search of the conjugate gradient
procedure the Hessian was checked and, if all the eigenvalues were positive,
the optimization was continued using the Newton-Raphson method. Beyond the
stable spin-moments of the iron atoms, we also included the induced moments of
the non-magnetic host atoms in the minimization of the grand potential.
Nevertheless, if the induced moment on a non-magnetic site was smaller than a
threshold value (chosen as $0.01\mu_{\rm B}$), the direction of the exchange
field at this site was set in the next iteration parallel to the direction of
the magnetization prescribed by the local spin-density approximation.

In all cases the optimization procedure has been done until the relative error
of the effective potential and of the magnitude of the exchange-correlation
field became smaller than $10^{-9}$ and $10^{-7}$, respectively and the
overall torque on the sites, $T=\displaystyle\sqrt{\sum_{i}T_{i}^{2}}$,
decreased below $10^{-2}\;\mathrm{mRyd}$.

###  2.3 Spin model calculations

We compared the results of the first principles optimization process to those
obtained from the following spin model,

$H=-\frac{1}{2}\sum\limits_{ij}\mathbf{s}_{i}\mathbf{J}_{ij}\mathbf{s}_{j}+\sum\limits_{i}\mathbf{s}_{i}\mathbf{K}_{i}\mathbf{s}_{i}\;,$
(15)

where $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ represent the direction of the
magnetic moments at site $i$ and $j$ of the magnetic cluster, $\mathbf{K}_{i}$
is the onsite anisotropy tensor and $\mathbf{J}_{ij}$ is a $3\times 3$ matrix,
called the tensorial exchange coupling. $\mathbf{J}_{ij}$ can be decomposed as
[9],

$\displaystyle\mathbf{J}_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\textrm{Tr}\left(\mathbf{J}_{ij}\right)\mathbf{I}+\frac{1}{2}\left(\mathbf{J}_{ij}-\mathbf{J}^{T}_{ij}\right)$
(16)
$\displaystyle+\frac{1}{2}\left(\mathbf{J}_{ij}+\mathbf{J}^{T}_{ij}-\frac{2}{3}\textrm{Tr}\left(\mathbf{J}_{ij}\right)\mathbf{I}\right),$

where $J_{ij}=\frac{1}{3}\textrm{Tr}\left(\mathbf{J}_{ij}\right)$ is the
isotropic Heisenberg exchange coupling,
$\mathbf{s}_{i}\frac{1}{2}\left(\mathbf{J}_{ij}-\mathbf{J}^{T}_{ij}\right)\mathbf{s}_{j}=\mathbf{D}_{ij}\left(\mathbf{s}_{i}\times\mathbf{s}_{j}\right)$
is the antisymmetric DMI, and the third, symmetric traceless part of
$\mathbf{J}_{ij}$ stands for the pseudo-dipolar anisotropy. These parameters
were obtained using the spin cluster expansion[48, 49, 50, 51], where the
self-consistent effective potentials and exchange-correlation fields of the
magnetic cluster and the neighboring host region were employed.

The ground state of the spin model was determined by consecutive Metropolis
and Landau-Lifshitz-Gilbert (LLG) spin dynamics simulations [38]. For each
system ten simulations were performed from independently chosen random initial
configurations. For the considered systems we found that all the ten
simulations ended in the same spin-configuration (or to its time-reversed
counterpart) within a relative error of $10^{-6}$.

##  3 Results

###  3.1 Fe/Rh(111)

Rh crystallizes in fcc structure with a lattice constant of
$a_{\mathrm{Rh}}=3.803\,$Å implying an inter-layer distance of 2.196 Å between
the Rh(111) planes. For the vertical position of the Fe chains with respect to
the uppermost Rh monolayer we used the interlayer distance between an Fe
monolayer and the top Rh monolayer of the Rh(111) substrate, $d_{\mathrm{Fe-
Rh}}=2.057\,$Å (corresponding to an inward relaxation of -7.3 %), reported by
Lehnert _et al._ [52]. The relaxation of the interlayer distances between the
Rh layers found to be 1 % or less are neglected in our calculations for
simplicity, see figure 2. It should be noted that in [52] the hexagonal close-
packed (hcp) stacking position of the Fe monolayer has been found slightly
lower in energy than the fcc stacking. In [53] it was demonstrated that the
magnetocrystalline anisotropy energy (MAE) is remarkably sensitive to the
stacking position of Fe and Co adatoms on Rh(111) and Pd(111) surfaces. Since
the MAE and, in general, relativistic effects are of great importance in
determining the magnetic state of supported nanoclusters, we investigated the
magnetic ground state of monatomic Fe chains occupying both fcc and hcp
stacking positions above the topmost Rh monolayer.

Figure 2:  Sketch of the simplified layer geometry for embedded Fe chains on
Rh(111) substrate. The Rh and Fe atoms, as well as the empty spheres in the
vacuum region are shown with blue, red and green spheres, respectively.

The electronic structure of the host Rh(111) surface was calculated using the
SKKR method[40]. The Rh(111) surface was modeled by eight monolayers of Rh and
four monolayers of empty spheres (vacuum) between semi-infinite bulk Rh and
semi-infinite vacuum regions. The vertical geometry of the layers was
identical with that shown in figure 2, the Fe monolayer being replaced by a
layer of empty spheres. The iron chains were then embedded into this vacuum
layer by means of the embedding technique based on the KKR Green’s function
method[41] and self-consistent calculations of the magnetic ground states have
been performed as described in Section II. We studied close-packed chains of
4, 5, 6 and 7 Fe atoms deposited in the $[1\overline{1}0]$ direction, where
the Fe atoms occupied either the fcc or the hcp hollow positions above the
Rh(111) surface.

Figure 3:  Magnetic ground states of monatomic chains of 4, 5, 6 and 7 Fe
atoms on a Rh(111) surface with fcc stacking obtained from ab-initio
optimization. The $x$, $y$ and $z$ axes correspond to the $[1\overline{1}0]$,
$[1\overline{1}\overline{2}]$ and $[111]$ directions, respectively Red spheres
represent the Fe atoms, while the blue ones represent atoms in the underlying
Rh layer. The red vectors show the direction of the spin magnetic moments of
the Fe atoms, while their magnitude in $\mu_{\mathrm{B}}$ is written inside
the red spheres. The panels on the right show the normal-to-chain views of the
spin-configurations of the four considered chains in the same order as in the
left panels.

In figure 3 we present the ground state magnetic configuration for the
considered Fe chains with fcc stacking. The magnitudes of the Fe magnetic
moments inside the chain show slight oscillations in the range of 3.17-3.18
$\mu_{\mathrm{B}}$, while the edge atoms exhibit a modestly enhanced moment of
3.26-3.27 $\mu_{\mathrm{B}}$ which can be attributed to the reduced
coordination of the edge Fe atoms. The magnitude of the induced magnetic
moments of the Rh atoms nearest the Fe atoms (not displayed in figure 3) is
around $0.2\,\mu_{\mathrm{B}}$.

Before discussing the ground-state spin-configurations displayed in figure 3
we note that the point group of the chains contains two operations: the
identity and a mirror plane perpendicular to the chain. In addition, time
reversal that turns around the spin vectors is also a symmetry operation of
the system. Accordingly, there exist two kinds of ground-state spin-
configurations being invariant either under the mirror transformation or under
simultaneous action of mirror transformation and time reversal. It then
follows that the spin vectors can transform either as axial or polar vectors
under the mirror plane operation.

As can be seen from figure 3a, either of the above symmetry of the magnetic
configurations can be obtained from the ab-initio optimization depending on
the length of the chains: the 4- and 7-atom-long chains show polar vector
symmetry, while the magnetic moments in the 5- and 6-atom-long chains behave
as axial vectors. Inferring the ground-state spin-configurations, it is
tempting that the building blocks of these configurations are pairs of Fe
spins being aligned closely parallel to each other ($\uparrow\uparrow$), while
these dimers seem to be coupled antiparallel to each other
($\uparrow\uparrow\downarrow\downarrow$). This structure is particularly
obvious for the chains containing even number of Fe atoms, as what follows
referred to as even chains. In case of chains containing odd number of Fe
atoms (odd chains) the alternation of $\uparrow\uparrow$ and
$\downarrow\downarrow$ dimers is necessarily broken, which gives rise to
strongly noncollinear spin-configurations. This might be associated with the
fact that the central spin in an odd chain must be aligned either normal to
the chain (polar vector symmetry) or parallel to the chain (axial vector
symmetry). Apparently, noncollinearity of the spins is present in the even
chains too. Furthermore, all the ground state spin-configurations are coplanar
with the planes tilted to the $z$ axis by roughly the same angle.

As what follows we make an attempt to analyze the main features of the ground-
state spin configurations for the different chains in terms of the bilinear
spin model (15). The leading terms in the spin model are the isotropic
exchange interactions. In Table 1 we present the nearest neighbor (NN) and
next nearest neighbor (NNN) isotropic interactions between the Fe atoms.
Apparently, the NN and NNN isotropic interactions are ferromagnetic (FM) and
antiferromagnetic (AFM), respectively. These interactions are consistent with
the alternating $\uparrow\uparrow\downarrow\downarrow$ structure seen as the
dominant feature of the even chains. For $n=6$ the large FM coupling between
the inner atoms, $J_{34}$, further stabilizes this spin structure. In the case
of odd number of atoms this double-pair-wise antiferromagnetic configuration
is obviously broken. It should be noted that a
$\uparrow\uparrow\downarrow\downarrow$ double-row-wise AFM structure has been
theoretically predicted [54] and experimentally observed [55] in an Fe
monolayer on Rh(111). While the emergence of such a spin-structure in a
monolayer requires the presence of higher-order (in particular, three-site
four-spin) interactions [55], for a chain of finite length it can be
stabilized by the inhomogeneous bilinear isotropic interactions only.

Table 1:  Nearest neighbor and next nearest neighbor isotropic exchange
interactions $J_{ij}$ (in mRy) between sites $i$ and $j$ in fcc-stacked
monatomic chains of $n$ Fe atoms. Label 1 denotes atoms at the left edge of
the chain. Only the independent interactions are shown.

  * $i-j$ | $n=4$ | $n=5$ | $n=6$ | $n=7$
---|---|---|---|---
$1-2$ |  -0.604 |  -0.572 |  -0.686 |  -0.605
$1-3$ | -0.562 | -0.624 | -0.614 | -0.599
$2-3$ |  -0.114 |  -0.331 |  -0.558 |  -0.322
$2-4$ |  | -0.441 | -0.456 | -0.406
$3-4$ |  |  |  -1.181 |  -0.227
$3-5$ |  |  |  | -0.420

Figure 4:  Ground-state spin configurations of a 4-atom-long Fe chain (left)
and a 5-atom-long Fe chain (right) in fcc-stacked positions on Rh(111) surface
based on spin models as follows: (a) isotropic exchange interactions and on-
site anisotropy, (b) the previous spin model extended with DMI, and (c) full
tensorial interactions.

For the cases of the 4-atom-long and 5-atom-long chains we demonstrate the
effect of different kinds of exchange interactions in forming the magnetic
ground state. For this reason, we consecutively considered three sets of spin-
model parameters: (a) isotropic interactions and on-site anisotropy matrices,
(b) the previous spin-model extended by DMI and (c) the spin-model with full
tensorial interactions (15). The spin-configurations obtained by spin-dynamic
simulations using these spin-models are presented in figure 4.

In case of the 4-atom-long chain, see figure 4 left panel, when only isotropic
exchange interactions and on-site anisotropies are present, the simulations
resulted in to a nearly collinear double-pair-wise AFM configuration along the
chain direction due to an easy $x$-axis anisotropy. This configuration clearly
respects a polar-vector symmetry of the spins. The $y$-components of the DM
vectors introduce non-collinearity of the spin-structure in the $x$-$z$ plane
by keeping the polar-vector symmetry. In addition, the plane of the spins is
slightly rotated to the $y$ axis, which shows the preference of the $y$ axis
with respect to the $z$ axis by on-site anisotropy (for a discussion of the
tilting of spin-spiral states see [25]). Switching on the symmetric exchange
anisotropy doesn’t considerably affect the spin-configuration.

As expected, in case of the 5-atom-long chain, see figure 4 right panel, the
frustration of the isotropic exchange interactions causes a strongly non-
collinear spin arrangement, with the middle spin oriented along the chain,
again due to easy $x$-axis anisotropy. The ground state exhibits, therefore,
an axial-vector symmetry. The $y-z$ view of the spin configuration shows that
the plane of the spins is slightly tilted away from the $x-y$ plane, which is
mainly due to the nonvanishing $xz$ components of the on-site anisotropy
matrices. Switching on the DMI further stabilizes this configuration by just
slightly changing the relative angles between the spins. On the other hand,
the plane of the spins is tilted close to the $x-z$ plane which can be
attributed to large $y$ components of the DM vectors [25]. As for the 4-atom-
long chain, including the symmetric exchange anisotropy has negligible effect
on the magnetic ground state.

Table 2:  The mean deviation $\Delta$, defined by (17) between the magnetic
ground states obtained from ab initio optimization and from the spin model
(15).

  * $n=4$ | $n=5$ | $n=6$ | $n=7$
---|---|---|---
0.332 | 0.417 | 0.297 | 0.601

The ground state configurations of the Fe chains on Rh(111) obtained from ab
initio optimization and based on the bilinear tensorial spin model compare
remarkably well. However, in particular, for the longer chains larger
differences occur mostly at the ends of the chains. In order to quantify the
differences between the magnetic configurations obtained from these two
methods we use the average mean deviation of the spin vectors $\Delta$, the
square of which is defined as

$\Delta^{2}=\frac{1}{n}{\sum_{i=1}^{n}(\mathbf{s}^{\mathrm{m}}_{i}-\mathbf{s}^{\mathrm{a}}_{i})^{2}}\;,$
(17)

with $\mathbf{s}^{\mathrm{m}}_{i}$ and $\mathbf{s}^{\mathrm{a}}_{i}$ being the
spin vectors at site $i$ obtained from the spin model and the ab initio
optimization, respectively. The mean deviations for the different chains are
summarized in Table 2. It is apparent that the agreement between the magnetic
ground states based on ab initio optimization and on the spin model is
significantly better for even chains which might be attributed to the missing
geometrical frustration, thus to reduced non-collinearity in these systems.

Figure 5:  Magnetic ground state of monatomic chains of 4, 5, 6 and 7 Fe atoms
on Rh(111) surface with hcp stacking obtained from ab-initio optimization. For
the details of the figure see the caption of figure 3. The fourth entry from
the top shows the magnetic ground state of the 7-atom-long chain with axial
vector symmetry, while the last entry depicts a metastable state of this chain
with polar vector symmetry.

As we mentioned before we also applied the ab initio optimization technique to
find the magnetic ground state of the Fe chains in hcp stacking. Our results
are shown in figure 5. The magnitudes of the spin moments are very similar to
those obtained for the chains with fcc stacking. Concerning the ground-state
spin configurations, the double pair-wise AFM structure keep dominating, which
indicates that the isotropic interactions are similar for both kinds of
stacking geometry of the Fe chains. Seemingly, the even chains remain almost
collinear. This behavior could be attributed to relatively weak DM
interactions compared to strong isotropic interactions. However, a spin-model
study as in case of the fcc-stacked chains yields a ground-state configuration
for the 6-atom with enhanced non-collinearity, implying that the bilinear spin
model (15) is less suitable to recover the magnetic ground state obtained from
the ab initio optimization.

As in case of fcc stacking the odd chains with hcp stacking show a largely
non-collinear magnetic ground-state, while the plane of the magnetic moments
is almost normal to the surface. For the 7-atom-long chain the magnetic ground
state shows axial-vector symmetry, see the fourth entry in figure 5. In
addition, we also found a metastable state with polar-vector symmetry depicted
in the last entry of figure 5, which is by $16.3\,\mathrm{meV}$ higher in
energy than the ground state. These two configurations are related to each
other by about 90∘ global rotation around the normal of the plane of the
spins, while the metastable configuration is slightly tilted from the $z$
axis. Note that both configurations have right-handed chirality due to
positive $y$ components of the DM vectors.

###  3.2 Fe chains on Nb(110) and Re(0001)

In this section we present our main results on the magnetic structure of Fe
chains on top of Nb(110) and Re(0001) surfaces. These systems are of
particular interest, since they are supposed to host Majorana edge states[4,
5, 6]. Bulk Nb has a bcc structure with lattice constant of
$a_{\mathrm{Nb}}=3.3004$ Å. Using the Vienna Ab-initio Simulation Package
(VASP)[56] Laszloffy et al.[39] found that by putting an Fe adatom on the
(110) surface of Nb, the average vertical distance of the atoms in the two
uppermost Nb layers decreased to $2.2542$ Å from the bulk value of $2.3337$ Å
, while the vertical distance between the Fe adatom and the closest Nb atoms
decreased to $1.9868$ Å. In the calculations of atomic clusters of Fe these
values were used for the interlayer distances between the uppermost two Nb
layers, as well as between the top Nb layer and the first two vacuum layer, in
which the Fe atoms were embedded[39], respectively. We adopted this geometry
in our present calculations of the Fe chains in terms of the embedding Green’s
function technique. The host surface was modelled by eight Nb layers and four
empty sphere layers sandwiched between semi-infinite bulk Nb and vacuum
regions. The top view of an Fe chain on the Nb(110) surface is shown in figure
6.

Figure 6:  Top view of an array of Fe atoms on top of a Nb(110) surface
studied in the present work. The Fe and Nb atoms are marked by red circles,
the blue circles stand for the Nb atoms. The $x$, $y$ and $z$ axes correspond
to the $[1\overline{1}0]$, [001] and $[110]$ directions, respectively. The
chains are deposited along the $[1\overline{1}0]$ direction.

We optimized the magnetic ground state of 5-, 10- and 15-atom long Fe chains
deposited along the $[1\overline{1}0]$ ($x$) direction. Similar to the spin-
model study in [39], we found that the ground state of the investigated chains
is a spin-spiral. In figure 7 we present two optimized magnetic configurations
for the 10-atom long chain. Remarkably, the spin-moments of the Fe atoms is by
about 1 $\mu_{\textrm{B}}$ less than in the chains on the Rh(111) surface
which is a typical volume effect. Due to the fact that the hard axis of the
system is $y$, the spins rotate in the $x$-$z$ plane. Based on a spin model
with exchange interactions obtained from first principles, in [39] a spin
spiral wavelength of $\lambda$=3.39 $a_{x}$ was found, where
$a_{x}=\sqrt{2}\,a_{\mathrm{Nb}}$ is the distance between the subsequent iron
atoms in the chain. For the 10-atom-long chain, we performed a least square
fit of the rotation angles of the spins with respect to $x$ axis to a linear
function and we obtained a slope of 2.17 rad/site which translates into a
spin-spiral wavelength of $\lambda$=2.90 $a_{x}$, which means a decrease by
about 15% as compared to [39].

It is obvious that the metastable state shown in figure 7a and the ground
state shown in figure 7b have opposite rotational sense. The ground state is
lower in energy by $1.7\,\mathrm{meV}$ than the metastable state, converted
into about $0.08\,\mathrm{meV}$/Fe, which is indeed in the energy range of the
DMI of these chains[39]. It is pointed out in [39] that the frustrated
isotropic interactions and the relatively strong anisotropy can stabilize spin
spirals with opposite chirality for longer chains. It should be noted,
however, that the ground state spin-spiral in figure 7b has a left-handed
(anti-clockwise) rotational sense as opposed to that obtained from the spin
model[39] which has a right-handed (clockwise) rotational sense. This
discrepancy might be attributed to higher-order chiral spin-spin
interactions[25, 57] missing in the spin-model used in [39].

Figure 7:  (a) Metastable and (b) ground state magnetic configurations
obtained by ab initio optimization for a 10-atom long Fe chain on Nb(110). Red
and blue circles represent the Fe and Nb atoms, respectively. The red
(unit)vectors show the direction of the magnetization vector, while its
magnitude is written inside the red spheres in $\mu_{\mathrm{B}}$ units. The
$x$, $y$ and $z$ directions correspond in order to the $[1\overline{1}0]$,
[001] and $[110]$ directions. The pictures on the right show the $y$-$z$ side-
view of the spin vectors. Note that the rotational sense of spin-spiral (a) is
right-handed (clockwise) around the positive $y$ axis, while spin-spiral (b)
has a left-handed (anti-clockwise) rotational direction.

Finally, we turn to Fe chains deposited on top of the (0001) surface of hcp
Rhenium. In our calculations we used the bulk in-plane lattice constant of
$a_{\mathrm{2D}}=2.761$ Å and the interlayer relaxations found by Lászlóffy et
al. [25] in terms of the VASP code[56]. As compared to the interlayer distance
$2.281$ Å of the (0001) planes in bulk Re, the distance between the uppermost
Re layer and the Fe adatoms decreased to $2.228\,$Å, while the distance
between the surface and subsurface Re layers reduced to $2.16\,$Å. The
Re(0001) surface was modeled as an interface region between a semi-infinite
bulk Re and vacuum consisting of 8 monolayers of Re and four monolayers of
empty spheres (vacuum). The iron chains were embedded into the lowermost
vacuum layer along the nearest neighbor direction ($x$).

The ground state magnetic configuration of a 15-atom-long chain of Fe atoms
obtained via ab initio optimization is shown in figure 8. This ground state is
obviously a similar spin-spiral state as found by Laszloffy et al.[25] also by
ab initio spin-dynamics using a cut-off $\ell_{\mathrm{max}}=2$ for the
multiple scattering calculations (see figure 5(e) of [25]). A linear fit of
the rotational angles of the spin vectors in figure 8 yields a spin-spiral
wavelength of about $3.4\,a_{\mathrm{2D}}$. Calculations on chains of 11, 13,
17 and 19 Fe atoms resulted in similar wavelengths. This wavelength is
considerably lower than the value of about $5\,a_{\mathrm{2D}}$ obtained in
[25]. Nevertheless, both values are close to $\lambda=4\,a_{\mathrm{2D}}$
found experimentally using spin polarized scanning tunneling microscopy by Kim
et al.[4] for a 40-atom long chain.

Figure 8:  Ground-state spin-configuration of a 15-atom long Fe chain on
Re(0001) substrate. The red spheres correspond to the Fe atoms, while the blue
ones belong the underlying Re layer. The $x$, $y$ and $z$ directions
correspond to the $[1000]$, $[1\overline{2}00]$ and $[0001]$ directions of the
hcp crystal structure, respectively. Red vectors show the directions of the
atomic magnetic moments, while their magnitudes are written inside the spheres
in $\mu_{\mathrm{B}}$ units.

##  4 Summary

We performed fully converged unconstrained self-consistent calculations of the
magnetic ground state of Fe chains deposited on different metallic surfaces by
using a newly develepod method based on a combination of the conjugate
gradient and Newton-Raphson methods within the multiple scattering Green’s
function embedding technique. In case of Fe chains in fcc stacking positions
on top of Rh(111) we found that the building block of the obtained magnetic
configurations is a double pair-wise antiferromagnetic structure that subsists
for chains of even number of Fe atoms, while for odd number of Fe atoms the
geometrical frustration leads to strongly non-collinear magnetic ground
states. We established that, as the length of the chain is varied, spin-
configurations with either axial or polar vector symmetry occur. Simulations
based on a tensorial spin-model showed that the isotropic exchange couplings
and the onsite anisotropy sufficiently explain the main features of the
magnetic ground state, but the Dzyaloshinky-Moriya interactions play an
important role in the formation of the final, complex spin-configuration. A
numerical comparison of the ground state spin-configurations from the two
methods showed larger deviations for higher degree of complexity (non-
collinearity) of the magnetic states indicating that higher order spin-spin
interactions should be included in spin-model studies. Our calculations
revealed that the Fe chains in hcp stacking positions have similar spin-
structures as in fcc stacking. For the 7-atom long chain the magnetic ground
state had axial vector symmetry, but we also found a metastable state with
polar vector symmetry stabilized by DM interactions.

In case of Fe chains deposited on Nb(110) surface along the second nearest
neighbor direction, we found almost degenerate spin-spiral states with
opposite chiralities in line with a recent spin-model study by Laszloffy et
al.[39], who demonstrated that the quasi-degeneracy of these spin-spiral
states is the consequence of the frustration of isotropic interactions and the
DM interactions, while spin-spirals not preferred by the DMI are stabilized by
strong magnetic anisotropy. For Fe chains on Re(0001) surface, our first
principles optimization process led to spin spiral states with a wavelength
comparable with the experimental finding[4].

In conclusion, we demonstrated that the optimization method introduced in this
work is an efficient tool to explore the complex magnetic ground states of
supported atomic clusters. This might be particularly helpful for a material
specific design of new magnetic nanostructures for quantum information
technology based on low-dimensional topological superconductivity.

The authors are thank András Lászlóffy for the help in the spin model
calculations and Levente Rózsa for fruitful discussions. This research was
funded by the Ministry of Innovation and the National Research, Development,
and Innovation Office of Hungary under the project No. K131938, project No.
FK124100 and under Grant TKP2021-NVA-02.

## References

## References

  * [1] Eigler D and Schweizer E 1990 Nature 344(18) 524–526 URL https://doi.org/10.1038/344524a0
  * [2] Meier F, Zhou L, Wiebe J and Wiesendanger R 2008 Science 320 82–86 ISSN 0036-8075 (Preprint https://science.sciencemag.org/content/320/5872/82.full.pdf) URL https://science.sciencemag.org/content/320/5872/82
  * [3] Andolina G M and Simon P 2017 Physical Review B 96 URL https://doi.org/10.1103/physrevb.96.235411
  * [4] Kim H, Palacio-Morales A, Posske T, Rózsa L, Palotás K, Szunyogh L, Thorwart M and Wiesendanger R 2018 Science Advances 4 (Preprint https://advances.sciencemag.org/content/4/5/eaar5251.full.pdf) URL https://advances.sciencemag.org/content/4/5/eaar5251
  * [5] Beck P, Schneider L, Rózsa L, Palotás K, Lászlóffy A, Szunyogh L, Wiebe J and Wiesendanger R 2021 Nature Communications 12 1–9
  * [6] Crawford D, Mascot E, Shimizu M, Schneider L, Beck P, Wiebe J, Wiesendanger R, Jeschke H O, Morr D K and Rachel S 2021 Majorana modes with side features in magnet-superconductor hybrid systems (Preprint 2109.06894)
  * [7] Küster F, Brinker S, Hess R, Loss D, Parkin S, Klinovaja J, Lounis S and Sessi P 2021 Non-majorana zero energy modes in diluted spin chains proximitized to a superconductor (Preprint 2112.05708)
  * [8] Choi D J, Lorente N, Wiebe J, von Bergmann K, Otte A F and Heinrich A J 2019 Reviews of Modern Physics 91 URL https://doi.org/10.1103/revmodphys.91.041001
  * [9] Udvardi L, Szunyogh L, Palotás K and Weinberger P 2003 Phys. Rev. B 68(10) 104436 URL https://link.aps.org/doi/10.1103/PhysRevB.68.104436
  * [10] Dzyaloshinsky I 1958 Journal of Physics and Chemistry of Solids 4 241 – 255 ISSN 0022-3697 URL http://www.sciencedirect.com/science/article/pii/0022369758900763
  * [11] Moriya T 1960 Phys. Rev. 120(1) 91–98 URL https://link.aps.org/doi/10.1103/PhysRev.120.91
  * [12] Bode M, Heide M, Von Bergmann K, Ferriani P, Heinze S, Bihlmayer G, Kubetzka A, Pietzsch O, Blügel S and Wiesendanger R 2007 Nature 447 190–193
  * [13] Udvardi L, Antal A, Szunyogh L, Buruzs Á and Weinberger P 2008 Physica B: Condensed Matter 403 402–404 ISSN 0921-4526 URL https://www.sciencedirect.com/science/article/pii/S092145260700600X
  * [14] Menzel M, Mokrousov Y, Wieser R, Bickel J E, Vedmedenko E, Blügel S, Heinze S, von Bergmann K, Kubetzka A and Wiesendanger R 2012 Phys. Rev. Lett. 108(19) 197204 URL https://link.aps.org/doi/10.1103/PhysRevLett.108.197204
  * [15] Schweflinghaus B, Zimmermann B, Heide M, Bihlmayer G and Blügel S 2016 Phys. Rev. B 94(2) 024403 URL https://link.aps.org/doi/10.1103/PhysRevB.94.024403
  * [16] Mühlbauer S, Binz B, Jonietz F, Pfleiderer C, Rosch A, Neubauer A, Georgii R and Böni P 2009 Science 323 915–919 URL https://doi.org/10.1126/science.1166767
  * [17] Münzer W, Neubauer A, Adams T, Mühlbauer S, Franz C, Jonietz F, Georgii R, Böni P, Pedersen B, Schmidt M, Rosch A and Pfleiderer C 2010 Phys. Rev. B 81(4) 041203 URL https://link.aps.org/doi/10.1103/PhysRevB.81.041203
  * [18] Yu X Z, Onose Y, Kanazawa N, Park J H, Han J H, Matsui Y, Nagaosa N and Tokura Y 2010 Nature 465 901–904 URL https://doi.org/10.1038/nature09124
  * [19] Yu X Z, Kanazawa N, Onose Y, Kimoto K, Zhang W Z, Ishiwata S, Matsui Y and Tokura Y 2010 Nature Materials 10 106–109 URL https://doi.org/10.1038/nmat2916
  * [20] Kiselev N, Bogdanov A, Schäfer R and Rößler U 2011 Journal of Physics D: Applied Physics 44 392001
  * [21] Sampaio J, Cros V, Rohart S, Thiaville A and Fert A 2013 Nature nanotechnology 8 839–844
  * [22] Romming N, Hanneken C, Menzel M, Bickel J E, Wolter B, von Bergmann K, Kubetzka A and Wiesendanger R 2013 Science 341 636–639
  * [23] Kurz P, Bihlmayer G, Hirai K and Blügel S 2001 Phys. Rev. Lett. 86(6) 1106–1109 URL https://link.aps.org/doi/10.1103/PhysRevLett.86.1106
  * [24] Heinze S, von Bergmann K, Menzel M, Brede J, Kubetzka A, Wiesendanger R, Bihlmayer G and Blügel S 2011 Nature Physics 7 713–718 URL https://doi.org/10.1038/nphys2045
  * [25] Lászlóffy A, Rózsa L, Palotás K, Udvardi L and Szunyogh L 2019 Phys. Rev. B 99(18) 184430 URL https://link.aps.org/doi/10.1103/PhysRevB.99.184430
  * [26] Hoffmann M and Blügel S 2020 Phys. Rev. B 101(2) 024418 URL https://link.aps.org/doi/10.1103/PhysRevB.101.024418
  * [27] Krönlein A, Schmitt M, Hoffmann M, Kemmer J, Seubert N, Vogt M, Küspert J, Böhme M, Alonazi B, Kügel J, Albrithen H A, Bode M, Bihlmayer G and Blügel S 2018 Physical Review Letters 120 URL https://doi.org/10.1103/physrevlett.120.207202
  * [28] Bergman A, Nordström L, Burlamaqui Klautau A, Frota-Pessôa S and Eriksson O 2006 Phys. Rev. B 73(17) 174434 URL https://link.aps.org/doi/10.1103/PhysRevB.73.174434
  * [29] Bergman A, Nordström L, Klautau A B, Frota-Pessôa S and Eriksson O 2006 Surface Science 600 4838–4842 ISSN 0039-6028 URL https://www.sciencedirect.com/science/article/pii/S0039602806008491
  * [30] Bergman A, Nordström L, Burlamaqui Klautau A, Frota-Pessôa S and Eriksson O 2007 Phys. Rev. B 75(22) 224425 URL https://link.aps.org/doi/10.1103/PhysRevB.75.224425
  * [31] Cardias R, Bezerra-Neto M M, Ribeiro M S, Bergman A, Szilva A, Eriksson O and Klautau A B 2016 Phys. Rev. B 93(1) 014438 URL https://link.aps.org/doi/10.1103/PhysRevB.93.014438
  * [32] de Melo Rodrigues D C, Pereiro M, Bergman A, Eriksson O and Klautau A B 2016 Journal of Physics: Condensed Matter 29 025807 URL https://doi.org/10.1088/0953-8984/29/2/025807
  * [33] Stocks G M, Ujfalussy B, Wang X, Nicholson D M C, Shelton W A, Wang Y, Canning A and Györffy B L 1998 Philosophical Magazine B 78 665–673 (Preprint https://doi.org/10.1080/13642819808206775) URL https://doi.org/10.1080/13642819808206775
  * [34] Újfalussy B, Wang X D, Nicholson D M C, Shelton W A, Stocks G M, Wang Y and Gyorffy B L 1999 Journal of Applied Physics 85 4824–4826 (Preprint https://doi.org/10.1063/1.370494) URL https://doi.org/10.1063/1.370494
  * [35] Újfalussy B, Lazarovits B, Szunyogh L, Stocks G M and Weinberger P 2004 Phys. Rev. B 70(10) 100404 URL https://link.aps.org/doi/10.1103/PhysRevB.70.100404
  * [36] Gambardella P, Dallmeyer A, Maiti K, Malagoli M, Eberhardt W, Kern K and Carbone C 2002 Nature 416 301–304
  * [37] Balogh L, Palotás K, Udvardi L, Szunyogh L and Nowak U 2012 Phys. Rev. B 86(2) 024406 URL https://link.aps.org/doi/10.1103/PhysRevB.86.024406
  * [38] Rózsa L, Udvardi L and Szunyogh L 2014 Journal of Physics: Condensed Matter 26 216003
  * [39] Lászlóffy A, Palotás K, Rózsa L and Szunyogh L 2021 Nanomaterials 11 ISSN 2079-4991 URL https://www.mdpi.com/2079-4991/11/8/1933
  * [40] Szunyogh L, Újfalussy B and Weinberger P 1995 Phys. Rev. B 51(15) 9552–9559 URL https://link.aps.org/doi/10.1103/PhysRevB.51.9552
  * [41] Lazarovits B, Szunyogh L and Weinberger P 2002 Phys. Rev. B 65(10) 104441 URL https://link.aps.org/doi/10.1103/PhysRevB.65.104441
  * [42] Jansen H J F 1999 Phys. Rev. B 59(7) 4699–4707 URL https://link.aps.org/doi/10.1103/PhysRevB.59.4699
  * [43] Lloyd P 1967 Proceedings of the Physical Society 90 207–216
  * [44] Honkela A, Tornio M, Raiko T and Karhunen J 2008 Natural conjugate gradient in variational inference Neural Information Processing ed Ishikawa M, Doya K, Miyamoto H and Yamakawa T (Berlin, Heidelberg: Springer Berlin Heidelberg) pp 305–314 ISBN 978-3-540-69162-4
  * [45] Smith S T 2013 Geometric optimization methods for adaptive filtering (Preprint 1305.1886)
  * [46] Edelman A, Arias T and Smith S 1998 SIAM Journal on Matrix Analysis and Applications 20
  * [47] Polak E and Ribiere G 1969 Rev. Francaise Informat Recherche Operationelle 3(16) 35–43
  * [48] Drautz R and Fähnle M 2004 Phys. Rev. B 69(10) 104404 URL https://link.aps.org/doi/10.1103/PhysRevB.69.104404
  * [49] Szunyogh L, Udvardi L, Jackson J, Nowak U and Chantrell R 2011 Phys. Rev. B 83(2) 024401 URL https://link.aps.org/doi/10.1103/PhysRevB.83.024401
  * [50] Deák A, Szunyogh L and Ujfalussy B 2011 Phys. Rev. B 84(22) 224413 URL https://link.aps.org/doi/10.1103/PhysRevB.84.224413
  * [51] Lászlóffy A, Udvardi L and Szunyogh L 2017 Phys. Rev. B 95(18) 184406 URL https://link.aps.org/doi/10.1103/PhysRevB.95.184406
  * [52] Lehnert A, Dennler S, Błoński P, Rusponi S, Etzkorn M, Moulas G, Bencok P, Gambardella P, Brune H and Hafner J 2010 Phys. Rev. B 82(9) 094409 URL https://link.aps.org/doi/10.1103/PhysRevB.82.094409
  * [53] Błoński P, Lehnert A, Dennler S, Rusponi S, Etzkorn M, Moulas G, Bencok P, Gambardella P, Brune H and Hafner J 2010 Phys. Rev. B 81(10) 104426 URL https://link.aps.org/doi/10.1103/PhysRevB.81.104426
  * [54] Hardrat B, Al-Zubi A, Ferriani P, Blügel S, Bihlmayer G and Heinze S 2009 Phys. Rev. B 79(9) 094411 URL https://link.aps.org/doi/10.1103/PhysRevB.79.094411
  * [55] Krönlein A, Schmitt M, Hoffmann M, Kemmer J, Seubert N, Vogt M, Küspert J, Böhme M, Alonazi B, Kügel J, Albrithen H A, Bode M, Bihlmayer G and Blügel S 2018 Phys. Rev. Lett. 120(20) 207202 URL https://link.aps.org/doi/10.1103/PhysRevLett.120.207202
  * [56] Kresse G and Furthmüller J 1996 Phys. Rev. B 54(16) 11169–11186 URL https://link.aps.org/doi/10.1103/PhysRevB.54.11169
  * [57] Brinker S, dos Santos Dias M and Lounis S 2019 New Journal of Physics 21 083015 URL https://doi.org/10.1088/1367-2630/ab35c9