Spirit: Multifunctional Framework for Atomistic Spin Simulations

TL;DR

Spirit is a comprehensive, user-friendly framework for atomistic spin simulations that supports various interactions, methods, and visualizations, enabling detailed analysis of magnetic systems at the atomic scale.

Contribution

It introduces a versatile, GPU-accelerated simulation framework with advanced methods for studying static and dynamic magnetic phenomena at the atomic level.

Findings

Validated methods on systems like vortices and skyrmions

Supports detailed thermodynamical and dynamical simulations

Provides efficient CPU and GPU parallelization

Abstract

The \textit{Spirit} framework is designed for atomic scale spin simulations of magnetic systems of arbitrary geometry and magnetic structure, providing a graphical user interface with powerful visualizations and an easy to use scripting interface. An extended Heisenberg type spin-lattice Hamiltonian including competing exchange interactions between neighbors at arbitrary distance, higher-order exchange, Dzyaloshinskii-Moriya and dipole-dipole interactions is used to describe the energetics of a system of classical spins localised at atom positions. A variety of common simulations methods are implemented including Monte Carlo and various time evolution algorithms based on the Landau-Lifshitz-Gilbert equation of motion, which can be used to determine static ground state and metastable spin configurations, sample equilibrium and finite temperature thermodynamical properties of magnetic materials and nanostructures or calculate dynamical trajectories including spin torques induced by stochastic temperature or electric current. Methods for finding the mechanism and rate of thermally assisted transitions include the geodesic nudged elastic band method, which can be applied when both initial and final states are specified, and the minimum mode following method when only the initial state is given. The lifetime of magnetic states and rate of transitions can be evaluated within the harmonic approximation of transition-state theory. The framework offers performant CPU and GPU parallelizations. All methods are verified and applications to several systems, such as vortices, domain walls, skyrmions and bobbers are described.

Tables

Table 1. Table 1: Parameters of the Gaussians in the energy surface of the single-spin system shown in Fig. 7 in the main text.

a	$\sigma$	$c_x$	$c_y$	$c_z$
$-1.10$	$0.06$	$-0.20$	$0.00$	$-0.90$
$0.80$	$0.15$	$-1.00$	$0.20$	$-0.20$
$-0.90$	$0.10$	$1.00$	$-0.20$	$-0.10$
$0.09$	$0.03$	$0.80$	$0.50$	$-0.80$
$0.15$	$0.07$	$0.80$	$-0.50$	$-0.70$
$-0.90$	$0.10$	$0.50$	$1.20$	$-0.40$
$-0.90$	$0.10$	$0.20$	$-0.90$	$-0.40$

Equations

$$\mathcal{H}=$$

$$-\sum _{⟨ij⟩}{J}_{ij}{\mathbf{n}}_i\cdot {\mathbf{n}}_j-\sum _{⟨ij⟩}{\mathbf{D}}_{ij}\cdot ({\mathbf{n}}_i\times {\mathbf{n}}_j)$$

$$+\frac{1}{2}\frac{{\mu }_0}{4\pi }\sum _{\begin{array}{c}i,j\\ i\ne j\end{array}}{\mu }_i{\mu }_j\frac{({\mathbf{n}}_i\cdot {\hat{\mathbf{r}}}_{ij})({\mathbf{n}}_j\cdot {\hat{\mathbf{r}}}_{ij})-{\mathbf{n}}_i{\mathbf{n}}_j}{{r_{ij}}^3},$$

$$E_{\mathrm{Quad}}=-\sum _{ijkl}{K}_{ijkl}\left({\mathbf{n}}_i\cdot {\mathbf{n}}_j\right)\left({\mathbf{n}}_k\cdot {\mathbf{n}}_l\right).$$

$$\nu =\frac{\left|{\oint }_C\mathbf{n}\cdot \mathrm{d}\mathbf{s}\right|}{{\oint }_C|\mathrm{d}\mathbf{s}|}.$$

$$Q(\mathbf{n})=\frac{1}{4\pi }{\int }_{{\mathbb{R}}^2}\mathbf{n}\cdot ({\partial }_x\mathbf{n}\times {\partial }_y\mathbf{n})\mathrm{d}\mathbf{r}.$$

$$v=\frac{1}{2\pi }{\oint }_{\Gamma }\frac{n_x\mathrm{d}{n}_y-n_y\mathrm{d}{n}_x}{n_x^2+n_y^2}$$

$$P=e^{-\Delta E/k_{\mathrm{B}}T},$$

$$\chi =\frac{1}{k_{\mathrm{B}}T}\left(⟨m^2⟩-{⟨m⟩}^2\right),$$

$$C_{\mathrm{V}}=\frac{1}{k_{\mathrm{B}}{T}^2}\left(⟨E^2⟩-{⟨E⟩}^2\right),$$

$$U_4=1-\frac{⟨m^4⟩}{3{⟨m^2⟩}^2}.$$

$$\frac{\partial {\mathbf{n}}_i}{\partial t}=$$

$$-\frac{\gamma \alpha }{(1+{\alpha }^2){\mu }_i}{\mathbf{n}}_i\times ({\mathbf{n}}_i\times {\mathbf{B}}_i^{\mathrm{eff}})$$

$$-\frac{\alpha -\beta }{(1+{\alpha }^2)}u{\mathbf{n}}_i\times ({\mathbf{j}}_e\cdot {\nabla }_{\mathbf{r}}){\mathbf{n}}_i$$

$$+\frac{1+\beta \alpha }{(1+{\alpha }^2)}u{\mathbf{n}}_i\times ({\mathbf{n}}_i\times ({\mathbf{j}}_e\cdot {\nabla }_{\mathbf{r}}){\mathbf{n}}_i),$$

$${\mathbf{B}}_i^{\mathrm{th}}(t,T)=\sqrt{2D_i(T)}{\mathbf{\eta }}_i(t)=\sqrt{2\alpha k_{\mathrm{B}}T\frac{{\mu }_i}{\gamma }}{\mathbf{\eta }}_i(t),$$

$$n_z(t)$$

$$\phi (t)$$

$$n_x(t)$$

$$\mathcal{H}=-\sum _i{K}_1{n}_{ix}^2+K_2{n}_{iy}^2-J\sum _{⟨ij⟩}{\mathbf{n}}_i\cdot {\mathbf{S}}_j.$$

$$F_{\nu }^{\mathrm{tot}}=F_{\nu }^{\mathrm{S}}+F_{\nu }^{\mathrm{E}},$$

$$F_{\nu }^{\mathrm{S}}=(l_{\nu -1,\nu }-l_{\nu ,\nu +1}){\tau }_{\nu },$$

$$F_{\nu }^{\mathrm{E}}=-\nabla E_{\nu }+(\nabla E_{\nu }\cdot {\tau }_{\nu }){\tau }_{\nu },$$

$$F_{\nu }^{\mathrm{S},\mathrm{CI}}=0,F_{\nu }^{\mathrm{E},\mathrm{CI}}=-\nabla E_{\nu }+2(\nabla E_{\nu }\cdot {\tau }_{\nu }){\tau }_{\nu }.$$

$${\Gamma }^{\mathrm{HTST}}=\frac{v}{2\pi }{\Omega }_0{e}^{-\Delta E/k_{\mathrm{B}}T},$$

$${\Omega }_0$$

$$v$$

$$\mathcal{V}{|}^{2N}={\Lambda }^T{T}^T\mathcal{V}{|}^{3N}T\Lambda ,$$

$$F^{\mathrm{eff}}=F-2(F\cdot \hat{\lambda })\hat{\lambda },$$

$$H_{ij}=T_i^T{\bar{H}}_{ij}{T}_j-T_i^{T}I(x^j\cdot {\nabla }^j\bar{\mathcal{H}})T_j,$$

$$K=4J,\mathbf{D}=(0,0,1J),$$

$$Q=\frac{1}{4\pi }\sum _l{A}_l,$$

$$\cos \left(\frac{A_l}{2}\right)=\frac{1+{\mathbf{n}}_i\cdot {\mathbf{n}}_j+{\mathbf{n}}_i\cdot {\mathbf{n}}_k+{\mathbf{n}}_j\cdot {\mathbf{n}}_k}{\sqrt{2\left(1+{\mathbf{n}}_i{\mathbf{n}}_j\right)\left(1+{\mathbf{n}}_j{\mathbf{n}}_k\right)\left(1+{\mathbf{n}}_k{\mathbf{n}}_i\right)}}$$

$${\mathbf{n}}_i\cdot \left({\mathbf{n}}_j\times {\mathbf{n}}_k\right)=0,\text{and}\left|{\mathbf{n}}_i+{\mathbf{n}}_j+{\mathbf{n}}_k\right|\le 1.$$

Full text

I Introduction

Multiscale materials’ simulations have emerged as one of the most powerful and widespread assets in the quest for novel materials with optimal or target properties, functionalities and performance. Simulations are employed to narrow down the design continuum of devices, to decrease the effort required to design novel materials, to substitute experiments that seem unfeasible, to suggest and analyse experiments and provide understanding of the underlying physics on scales ranging from Ångström to millimeters and from femtoseconds to decades. In this context, spintronics is a very active field where multiscale simulations hoffmann_antiskyrmions_2017 play an important role for the conceptualization and development of the next-generation data devices, zutic_spintronics_2004 which includes nanoscale magnetic objects like domain walls or nontrivial magnetic textures such as solitons with a time dilemma of 16 orders of magnitude between writing and saving information. Relating the requested properties to the development and choice of magnetic materials, the simulation approach is highly useful, and a large variety of potential applications exist. bader_colloquium_2006

Since quantum mechanics is the key to understand magnetism, at the quantum mechanical level, ab initio methods, such as density functional theory, can be used to calculate various properties of and interactions between atoms. Due to the computational complexity of such calculations, they can currently only be applied to magnetic structures in crystals with length scales in the order of $1$ nm and cannot be used for time-dependent dynamics simulations on time-scales relevant for spintronics. From ab initio methods one may extract parameters for more approximate, atomistic spin models, such as Heisenberg type spin-lattice Hamiltonians. There, detailed information about the electronic structure is integrated out to effective parameters describing the interaction between pairs of classical spins, so that simulations of magnetization dynamics can be extended over the timescale of nanoseconds for systems of hundreds of nanometers. The third level of the multiscale approach in spintronics is the well-known micromagnetic approximation brown_micromagnetics_1963 based on the assumption of a continuous magnetization vector field, defined at any point of the magnetic sample, is valid when changes of the magnetization are much larger in space than the underlying atomic lattice. In contrast, the atomistic spin-lattice model covers the technologically increasingly important length scale from a few to several tens of nanometers.

Here, we present a general purpose, open source, i.e. publicly available, framework for atomistic spin simulations called Spirit. spirit There are actually a number of computational tools available for the simulation of the time- and space-dependent magnetization evolution. Among the software packages for micromagnetic simulations, two of the most impactful and widely known ones are OOMMF oommf and mumax3. mumax This software definitely revolutionized the simulation of magnetic properties of materials and the temporal behavior of devices described by the Landau Lifshitz Gilbert (LLG) dynamics. Based on the micromagnetic approach these methods have well-known limitations, e.g.  the description of antiferromagnets, frustrated magnets, higher order non-pairwise interactions (e.g. three-spin or four-spin interaction), stochastic spin dynamics and Monte Carlo simulations, etc. Moreover, most micromagnetic software is not interactive or provides quite limited in situ access to the parameters of the modeled system. Among the atomistic simulation programs, UppASD skubic_method_2008 and VAMPIRE evans_atomistic_2014 are examples of well-tested tools that provide important functionalities beyond LLG simulations.

The functionality of the aforementioned softwares can be greatly extended by adding an interactive graphical user interface (GUI) that can be used to control calculations in real time – to not only change parameters, but also interact with the spin texture as demonstrated for example in Ref. rybakov_chiral_2018, . Together with such a GUI, Spirit unifies various computational methods that are commonly applied to atomistic (and in part also to micromagnetic) simulations: Monte Carlo and Landau Lifshitz Gilbert (LLG) dynamics, nowak_thermally_2001 the geodesic nudged elastic band (GNEB) method bessarab_method_2015 , minimum mode following mueller_duplication_2018 (MMF), harmonic transition-state theory bessarab_harmonic_2012 (HTST), and the visualization of eigenmodes. All of these methods are quite distinct, but complementary in nature. braun_topological_2012 For example, LLG dynamics can be used to simulate the time evolution of a magnetic system on a short time scale, while GNEB and/or MMF can be used to find first-order saddle points of the energy landscape – corresponding to transition states for thermally assisted transitions. These calculations can provide important information, such as the energy barrier for the transition and can be used in HTST to calculate the lifetimes of metastable magnetic configurations over a long time scale. The integration of these methods into a single, uniform framework can lead to significant increase in productivity.

The following section will introduce the structure of the Spirit software and subsequent sections will detail the aforementioned methods, ordered by their complexity, which we rank according to the derivatives of the energy required for the implementation of these methods. Provided examples are mainly related to magnetic skyrmions, which represents one of the most rapidly developing fields in modern nanomagnetism.

II The Framework

The framework consists of modular components, as illustrated in Fig. 1: a core library for calculations and input/output; an application programming interface (API) layer to abstract the interaction with the code and provide a uniform interface across various programming languages, e.g. C/C$++$ and Python; a set of user interfaces to enable fast and easy interaction with simulations, powerful real-time visualization and post-processing features, for instance visualization of 2D and 3D magnetization vector fields with corresponding isosurfaces and visualization of eigenmodes.

The visualization of Spirit is available as a standalone library called VFRendering, vfrendering which utilizes advanced features of modern OpenGL, e.g. shaders, available since version 3.2. Note that the images of spin systems in Figs. 3 and 8 have been generated using the graphical user interface of Spirit. Other examples of the visualization features of Spirit can be found in Refs. liu_binding_2018, ; zheng_experimental_2018, ; du_interaction_2018, ; mueller_duplication_2018, . Spirit has also been used for numerical calculations in Refs. hoffmann_antiskyrmions_2017, ; hagemeister_controlled_2018, ; mueller_duplication_2018, ; redies_distinct_2018, .

As the API layer is written in the C programming language, many other languages can be used to call the corresponding functions and the core library can, therefore, be used in many different contexts. An illustration of this flexibility is the implementation of an additional, web based user interface, juspin using JavaScript to call Spirit and WebGL to display the simulated system. The desktop GUI can be used to control parameters in the calculation, such as system size or interaction parameters – useful for fast testing and setup – as well as for direct interaction with the spin textures. The latter is highly useful, for example to set up complex initial states rybakov_chiral_2018 or rectify calculations, such as GNEB paths that have diverged from their intended transition In order to increase productivity in repetitive or long-timescale calculations, Spirit can be used in Python scripts. Such scripts allow to reproduce all steps which can be taken in the GUI and therefore flexible and effective use of clusters and remote machines. Example Python scripts can be found in the code repository spirit . Note that the ability to use Spirit from Python also enables a straightforward integration into multiscale simulations and workflow automation frameworks, such as ASE ASE or AiiDA. AiiDA Documentation of input, features and the APIs, as well as examples of usage can be found online. spirit

The spin-lattice Hamiltonian, as well as all implemented methods and solvers have been abstracted from the specifics of numerical operations, allowing a generic backend, which can optionally use OpenMP for CPU parallelization or CUDA for GPU parallelization. The performance of a simple LLG simulation over different system sizes, including dipolar interactions, is shown in Fig. 2. The performance gain of the parallelization over the single-threaded case is obvious as $10$ cores give almost an order of magnitude across a broad range of system sizes and the GPU can give another order of magnitude at larger system sizes. As expected, the speed drops with the system size. Note that when dipolar interactions are included, due to the usage of FFTs, iterations can be slowed down if a side length of the system is not a power of two.

III Model and Methods

III.1 Hamiltonian

In Spirit, we implemented an extended Heisenberg Hamiltonian aharoni_introduction_2000 ; rado_magnetism_1963 of classical spins $\bm{n}_{i}$ of unit length located at lattice sites $i$ giving rise to the magnetic moment $\bm{m}_{i}=\mu_{i}\bm{n}_{i}$. The general form

[TABLE]

includes (i) the Zeeman term describing the interaction of the spins with the external magnetic field $\bm{B}$, (ii) the single-ion magnetic anisotropy, where $\hat{\bm{K}}_{j}$ are the axes of the uniaxial anisotropies of the basis cell with the anisotropy strength $K_{j}$, (iii) the symmetric exchange interaction given by $J_{ij}$ and the antisymmetric exchange, also called Dzyaloshinskii-Moriya interaction, given by vectors $\bm{D}_{ij}$, where $\braket{ij}$ denotes the unique pairs of interacting spins $i$ and $j$, (iv) the dipolar interactions, where $\hat{\bm{r}}_{ij}$ denotes the unit vector of the bond connecting two spins.

Quite often, the number of pairs for the exchange interactions is limited to nearest or next-nearest neighbors only. In Spirit the implementation of the Hamiltonian (1) does not assume any limitation on the number of or distance between such pairs, meaning that long-ranged and non-isotropic interactions can be considered.

Additionally, higher-order multi-spin-multi-site interactions hoffmann_systematic_2018 are implemented in Spirit as quadruplets of the form

[TABLE]

These can represent the 4-spin-2-site szilva_interatomic_2013 (also called biquadratic), the 4-spin-3-site, kroenlein_magnetic_2018 and the 4-spin-4-site heinze_spontaneous_2011 (also called "4-spin") interactions.

Both the system geometry and the underlying lattice symmetry can be chosen arbitrarily by setting the Bravais vectors and basis cells with any given number of atoms. Spirit also allows the pinning of individual spins or a set of spins, for instance belonging to the boundary layers. One can also introduce defects, such as vacancies and atoms of different types.

Dipolar interactions. The dipole-dipole interaction, due to its long-ranged nature, represents the most complex contribution to the Hamiltonian (1). Direct summation over all interacting spins is of complexity $\mathcal{O}(N^{2})$, where $N$ is the number of spins, resulting in dramatic decrease of performance of the simulations. By making use of fast Fourier transforms (FFTs) and the convolution theorem, the computational complexity can be reduced to $\mathcal{O}(N\log N)$. This convolution method is well-known in micromagnetic simulations, hayashi_ddi_1996 based on a finite difference scheme. To treat arbitrary spin lattices with any given number of atoms in the basis cells, we use an adapted version of this method. In particular, we consider sublattices composed of atoms with the same index in the basis cell. One FFT is performed on each of these sublattices and additional convolutions are required to describe the interactions between the sublattices. An efficient implementation of this scheme is achieved using high performance, robust FFT libraries. fftw ; cufft

To verify our implementation of dipolar interactions, we compared it to the direct evaluation of the sum for random configurations with spatially non-symmetric basis cells and checked the convergence of the stray field of a homogeneously magnetized monolayer against the analytically known result. Here, we show that Spirit correctly reproduces the solution of typical problems, e.g. Ref. hubert_systematic_1999, , by calculating the stray field-induced helicity of a ferromagnetic cube. The helicity is defined as the absolute value of the line integral over the curve $C$ which is composed of the upper edges of the cube:

[TABLE]

In the atomistic case this is discretized into the appropriate sums.

The energy minimization was performed using a Verlet-like velocity projection method (see Appendix D). The results are shown in Fig. 3. The squared helicity is expected to approach the critical field linearly, so that a line can be fitted to extract the precise result from the calculations. We note that the resulting critical field converges to the expected value of $h_{\mathrm{c}}=0.158$ with increasing resolution of the grid, where a cube with $30^{3}$ lattice sites already gives an agreement within $2\%$ and the shown example with $50^{3}$ sites a discretization error of less than $1\%$ with respect to the continuum solution.

Topological charge. In the past years, we witnessed the characterization of smooth magnetization fields $\bm{m}(\bm{r})=\mu\bm{n}(\bm{r})$ in terms of topological classes. In this case, the $\mathbb{S}^{2}$ winding number defined as

[TABLE]

It should be noted that a secondary topological charge $v$ can be considered as the defining index to distinguish between skyrmion ($v=1$) and antiskyrmion ($v=-1$) independently of the background state. It is defined as

[TABLE]

and is also referred to as the $\mathbb{S}^{1}$ winding number in comparison to the conventional topological charge $Q$. Here, $\Gamma$ denotes an oriented Jordan curve around the center of the skyrmion.

The notion of skyrmion/antiskyrmion refers to a local energy minimizer of (1) within the topological class characterized by the relative skyrmion number $N=\pm 1$. For sufficiently regular fields $\bm{n}$ decaying to the background state $n_{z}(\infty)=\pm 1$, the index $N$ is defined relative to this background state, $N(\bm{n})=-n_{z}(\infty)Q(\bm{n})$.hoffmann_antiskyrmions_2017 In a typical situation where the horizontal magnetization field vanishes at a single point (skyrmion center) the relative skyrmion number $N$ agrees with the index of the horizontal magnetization field at the skyrmion center. It is customary to fix the background state $n_{z}(\infty)=1$, which leads to the characterization $Q=-1$ for skyrmions and $Q=+1$ for antiskyrmions.

The evaluation of expression (4) for the spindensity $\bm{n}$ on a discrete lattice, as implemented in Spirit, is outlined in Appendix A.

III.2 Monte Carlo

The Monte Carlo method is well-known in Physics and has a broad range of applications. binder_monte_1981 We have implemented a basic Metropolis algorithm with a cone angle for the displacement of the spins. hinzke_monte_1999 ; nowak_thermally_2001 This requires only the calculation of the energy, making it the most straightforward of the methods implemented in Spirit. While it is a useful tool to calculate equilibrium properties, the drawback is that it cannot resolve time-dependent processes.

One iteration of the Metropolis algorithm will sequentially – but in random order – pick each spin in the system once and perform a trial step. Trial steps are preformed by defining a relative basis in which the current spin is the $z$-axis and choosing a new spin direction by uniformly distributed random variables $\varphi\in[0,2\pi]$ and $\cos(\theta)\in[0,\cos(\theta_{\mathrm{cone}})]$, where $\theta_{\mathrm{cone}}$ is the opening angle of the cone. The trial step is accepted with a probability

[TABLE]

where $\Delta E$ is the energy difference between the previous spin configuration and the trial step. The cone angle can be set by an adaptive feedback algorithm according to a desired acception-rejection ratio.

Using this method, one can, for example, calculate the critical temperature of a spin system. It is known that the isothermal susceptibility is related to the magnetization fluctuations landau_guide_2009

[TABLE]

where $m=\frac{1}{N}|\sum_{i}\bm{n}_{i}|$ is the average magnetization of the sample, while the specific heat relates to fluctuations of the energy

[TABLE]

where both should diverge at the critical temperature $T_{\mathrm{c}}$ for a phase transition, e.g. to the paramagnetic phase. The $4^{\mathrm{th}}$ order Binder cumulant binder_finite_1981 , which is often used to avoid finite size scaling effects, is defined as

[TABLE]

Fig. 4 shows these quantities as results of a Monte Carlo calculation of a cube of ${30\times 30\times 30}$ lattice sites for $J=1$ meV. For a simple cubic ferromagnet, from the high-temperature expansion method, baker_high-temperature_1967 the critical temperature is known to be rocio_modeling_2011 $T_{\mathrm{c}}=1.44~{}J/k_{\mathrm{B}}=16.71$ K. The results shown in Fig. 4 demonstrate the validity of the implementation, as the expected critical temperature is matched with an error of less than $1\%$.

Note that in Monte Carlo methods, the parallel tempering algorithm has proven to be an effective tool.  swendsen_replica_1986 ; hukushima_exchange_1996 ; bottcher_b-t_2017 The usage of Python and a MPI package would enable one to quite easily reproduce this algorithm in a Python script using Spirit.

III.3 Landau-Lifshitz-Gilbert Dynamics

The Landau-Lifshitz-Gilbert (LLG) equation landau_on_1935 ; gilbert_phenomenological_2004 is the well-established equation of motion for the dynamical propagation of classical spins. In its explicit form and including spin torque and temperature contributions, brown_thermal_1963 ; schieback_numerical_2007 it can be written

[TABLE]

in which the terms correspond to (i) precession, (ii) damping, (iii) precession-like current-induced spin torque, and (iv) damping-like current-induced spin torque. $\gamma$ is the electron gyromagnetic ratio, $\alpha$ is the damping parameter, $\bm{B}^{\mathrm{eff}}_{i}$ is the effective field, $\beta$ is a non-adiabatic parameter, $\bm{j}_{e}$ is the electron current vector, and $\nabla_{\bm{r}}=\partial/\partial\bm{r}$ is the spatial gradient acting on the spin orientation. The effective field always contains a component related to the energy gradient $\bm{B}^{\mathrm{eff}}_{i}=-\partial\mathcal{H}/\partial\bm{n}_{i}$, but in this notation for the LLG equation, the effective field may contain also a stochastic thermal field, i.e. $\bm{B}_{i}\to\bm{B}_{i}+\bm{B}_{i}^{\mathrm{th}}$, given by

[TABLE]

where the magnitude is given by the fluctuation-dissipation theorem and $\bm{\eta}_{i}$ is white noise, such that $\braket{\bm{B}_{i\alpha}^{\mathrm{th}}(t,T)}_{t}=0$ and $\braket{\bm{B}_{i\alpha}^{\mathrm{th}}(t,T)\bm{B}_{j\beta}^{\mathrm{th}}(0,T)}_{t}=2D_{i}(T)\delta_{ij}\delta_{\alpha\beta}\delta(t)$. To achieve these properties in an implementation, the vectors $\bm{\eta}_{i}(t)$ can each be created from three independent standard normally distributed random values at every time step. Note also that in time-integration schemes, to fulfill the fluctuation-dissipation relation, the thermal field needs to be normalized by the time step with a factor $1/\sqrt{\delta t}$. For more details on the integration of the stochastic LLG equation see for example references mentink_stable_2010, ; bauer_thermally_2011, ; levente_langevin_2014, and references therein.

Sampling of the stochastic LLG for the same parameters as shown in Fig. 4 is presented in Appendix E, verifying the implementation and the equivalence of the stochastic LLG and Monte Carlo methods.

In order to evolve a spin system in time according to this equation, quite a few well-known solvers can be applied. In Spirit, currently Heun’s method, nowak_thermally_2001 a $4^{\mathrm{th}}$ order Runge-Kutta solver, Depondt’s Heun-like method depondt_spin_2009 , and Mentink’s semi-implicit method B (SIB) mentink_stable_2010 are implemented (see Appendices B and C for details). These methods can also be used for energy minimization by considering only the damping part of the LLG equation. However, experience has shown that a Verlet-like velocity projection solver bessarab_method_2015 can greatly improve convergence to the closest energy minimum, as it carries a fictive momentum (see Appendix D for details).

An easy test for the validity of the implemented dynamics solvers is the Larmor precession and the damping of a single spin in an external magnetic field, as shown in Fig. 5. The analytical equations with which the results can be compared are

[TABLE]

The errors of the Depondt solver, shown in Fig. 5, match those of an equivalent calculation given in Ref. evans_atomistic_2014, .

In order to verify our implementation of spin current induced torques, the results from Ref. schieback_numerical_2007, on the velocity of a domain wall in a head-to-head spin chain were reproduced for various non-adiabatic parameters $\beta$. The chain is oriented along the $x$-axis and the first and the last spin are fixed in $+x$ and $-x$ direction, respectively. As a subset of the general Hamiltonian (1), the Hamiltonian for this example can be written as follows:

[TABLE]

The reference provides analytical equations against which the numerical results were checked. In Fig. 6 we show the data for the average domain wall velocity $\braket{v}$ over applied current $u$ in normalised units.

The approximate prediction thiaville_domain_2004 $\langle v\rangle=\sqrt{u^{2}-u_{c}^{2}}/(1+\alpha^{2})$ fits the results well, as shown in Fig. 6. As expected, we observe the Walker breakdown schryer_motion_1974 and a critical effective velocity of $u_{c}\approx 0.0414$, which is in close agreement with the reported value of $u_{c}\approx 0.0416$. Note, for $\beta=0.1$ and currents larger than $u_{W}$ and for $\beta=0$ and currents larger than $u_{c}$, the wall starts rotating around the $x$-axis.

III.4 Geodesic Nudged Elastic Band Method

When determining the rates of some rare transition events or the lifetimes of metastable magnetic states, LLG dynamics simulations typically are typically unfeasible due to the disparity between the time scales of the simulation and the transition events. An approach to this problem is given by a set of rate theory methods, namely the geodesic nudged elastic band bessarab_method_2015 (GNEB) and minimum mode following mueller_duplication_2018 (MMF) methods together with harmonic transition-state theory bessarab_harmonic_2012 (HTST). The latter two are higher order methods, requiring knowledge of the second derivatives of the energy – the Hessian matrix – and will be described in the following sections.

The GNEB method is a way of calculating minimum energy transition paths between two pre-determined configurations. The path is discretized by a number of spin configurations, in the following called images. In order to converge from an initial guess to a stable, energy-minimized path, spring forces are applied along the path tangents, while energy gradient forces are applied orthogonal to the path tangents. The total force therefore reads

[TABLE]

where $\nu$ is the image index along the chain, $F^{\mathrm{S}}$ is a spring force, and $F^{\mathrm{E}}_{\nu}$ is an energy gradient force. The forces in this section are $3N$-dimensional vectors. A simple definition of the spring force, which gives an equidistant distribution of images in phase space, is given by

[TABLE]

where $l_{\nu,\mu}$ is a measure of distance between images $\nu$ and $\mu$ and $\tau_{\nu}$ is the (normalized) path tangent at image $\nu$. The $F_{\nu}^{\mathrm{E}}$ should pull each image towards the minimum energy path, while leaving the distance to other images unchanged. They can be defined to be orthogonal to the path by orthogonalizing with respect to the tangents

[TABLE]

where $\nabla_{i}=\partial/\partial\bm{n}_{i}$. The path tangents can be easily approximated by finite differences between spin configurations, but in order to avoid the formation of kinks in the path the definitions given in Ref. henkelman_improved_2000, should be used.

In order to precisely find the point of highest energy along the minimum energy path, a first order saddle point of the energy landscape, one can use a so-called climbing image. henkelman_climbing_2000 Convergence onto the saddle point is achieved through the deactivation of the spring force for that image, while inverting the energy gradient force along the path:

[TABLE]

This will cause it to minimize all degrees of freedom, except the tangent to the path, which is instead maximized.

So far, the definitions match those of the regular NEB method. In order to use the NEB method for spin systems, it is necessary to consider the constraint of constant spin length and treat tangents and force vectors accordingly. bessarab_method_2015 For more details see Appendix F and H.

In order to verify and illustrate the GNEB method, we show the example of a single spin in a set of Gaussian potentials (see Appendix G). Fig. 7 shows the initial guess, made by homogeneous interpolation between the initial and final configuration, as well as a relaxed chain of images and a chain with two climbing and one falling image. The climbing images converge onto the saddle points and the falling image onto an additional local minimum, so that the energy barriers are known exactly.

The implementation of the GNEB method can be further tested using a conceptually simple process, which has enough degrees of freedom to pose a challenge for convergence: the destruction of a skyrmion tube in a chiral magnetic thin film. The parameter set is chosen in accordance with a calculation presented in Ref. rybakov_new_2015, , where a novel particle-like state is shown to emerge along the minimum energy path – the chiral bobber. The nucleation of a pair of Bloch points, cutting the skyrmion tube in half, is reported, resulting in the formation of one chiral bobber at each surface of the film. In fact, as we show in Fig. 8, also a single Bloch point can be nucleated at one of the films free surfaces. For these calculations the specific parameters are $J=1$ and $D=0.45~{}J$, meaning that the incommensurate spin spiral has a period of $L_{\mathrm{D}}=13.96~{}a$. We note that the conical phase background – corresponding to the ground state of the system – introduces additional modes with little energy cost associated and this can slow the rate of convergence to the minimum energy path. The climbing-image method henkelman_climbing_2000 was used to converge nearby images onto the maxima along the path and – analogous to what is suggested in the reference – the spring force was modulated to distribute images evenly along the energy curve. The latter improves the convergence onto the maxima, as the resolution for the finite-difference calculation of the tangents at the saddle points is increased. As it is common to calculate cubic polynomials to interpolate between the discrete points, the segment length of these polynomials can be used for the spring forces between the images. In Spirit, an additional parameter is implemented, with which one can set the weighting of energy versus reaction coordinate. Without the climbing image method, energy barrier calculations may be quite imprecise, especially when the resolution near the maximum is low. This is illustrated by the fact that we observe a ratio of the energy barriers between the collapse of the bobber and the Bloch point nucleation of only $3.3$, while Ref. rybakov_new_2015, – not using climbing images – reports a ratio of $4.3$.

The GNEB calculations reveal a crossover between the two Bloch point nucleation mechanisms, where at increasing field it becomes favorable to nucleate just one Bloch point at the surface. It can further be seen that the energy barrier for the collapse of the bobber goes to zero right below the critical field $H_{\mathrm{D}}$, meaning that – in the frame of this model – it can only be stabilised in the conical phase. In order to give additional quantitative reference results for this parameter set, the dependence of the energy barrier on the external magnetic field is also presented in Fig. 9.

III.5 Harmonic Transition-State Theory

As certain processes may be too rare or the desired time scale, which is to be simulated, too large to allow for dynamical simulations, other approaches are essential in estimating stability and the calculation of lifetimes of metastable magnetic states. One can employ the well-known transition-state theory, wigner_transition_1938 which has been used extensively, e.g. in chemical reaction and diffusion calculations. truhlar_current_1996 The rate of transitions can be estimated from the probability of finding the system in the most restrictive and least likely region separating the initial state from possible final states – the transition state, sometimes also called dividing surface. Within the harmonic approximation to transition-state theory bessarab_harmonic_2012 (HTST), one can make simplifications allowing the analytical calculation of the rate, which is then given by an Arrhenius-type law with an exponential dependence on the inverse temperature $T$ and the energy barrier of the transition $\Delta E$:

[TABLE]

where

[TABLE]

where the M and S superscripts indicate the minimum and first order saddle point of the transition. The $\lambda_{i}$ are eigenvalues of the Hessian matrix (see Appendix H), $V$ are the phase space volumes of zero modes (if present, otherwise $V=1$), $N_{0}$ are the number of zero modes – modes with zero eigenvalue – and $a_{i}$ are coefficients in the expansion of the velocity along the unstable mode. The primes next to determinants, products, and sums denote that only positive eigenvalues are taken into account.

The factors $a_{i}$ are in fact velocities: the first row of the dynamical matrix $\mathcal{V}$ transformed into the eigenbasis of the Hessian according to

[TABLE]

where $T$ is a ${3N\times 2N}$ basis matrix of the tangent space and $\Lambda_{ij}$ denote the matrix of the Hessians eigenvectors (in $2N$-representation, i.e. in the basis $T$). See Appendix H for more information.

The implementation has been verified against UppASD, skubic_method_2008 and we additionally present an example for the calculation of the lifetime $\tau$ of an isolated skyrmion in a two-dimensional system. As parameters we chose $J=1$ and $D=0.6~{}J$ and only the radial collapse mechanism is considered, making for a simple structure of the dependence on external field and temperature. Note that this example is purely illustrative and while larger skyrmions would exhibit longer lifetimes, the parameters are chosen to produce a small skyrmion in order to reduce the computational effort. Fig. 10 shows the results for an external field varied between $3.5$ T and $5$ T and temperature between $2$ K and $5$ K. HTST as well as Langers theory, langer_statistical_1969 which is closely related, have recently both been used to calculate skyrmion lifetimes, bessarab_lifetime_2018 ; desplat_thermal_2018 ; von_malottki_skyrmion_2018 showing that energy barriers are in general not enough to estimate the stability of metastable magnetic states.

There are two translational zero modes at the initial state minimum, while – due to the lattice discretisation and the defect-like shape of the skyrmion at the saddle point – there are no zero modes at the saddle point. Consequently, the transition rate prefactor has a linear temperature dependence.

III.6 Minimum Mode Following Method

To find the first order saddle points on the energy surface, without prior knowledge of the possible final states, the minimum mode following method mueller_duplication_2018 can be used. The effective force acting on a spin configuration is defined as

[TABLE]

where $F=-\nabla\mathcal{H}$ is the negative gradient of the energy and $\hat{\lambda}$ is the normalized eigenvector corresponding to the lowest, negative eigenvalue of the Hessian matrix of second derivatives. Note that these vectors and the dot product are $3N$-dimensional for a system with $N$ spins.

The calculation of second derivatives requires further attention, as the requirement of constant length effectively constrains the spins to a sub-manifold $\mathcal{M}^{\mathrm{phys}}\subset\mathcal{E}$ of an embedding space $\mathcal{E}=\mathbb{R}^{3N}$. As is shown in Ref. mueller_duplication_2018, , the covariant second derivatives, valid at all points of the phase space, can be calculated using a projector-based approach. absil_extrinsic_2013 The corresponding ${2N\times 2N}$ Hessian matrix can be represented as

[TABLE]

where $i$ and $j$ are spin indices, $\bar{\mathcal{H}}$ is the smooth continuation of the Hamiltonian to the embedding space, $\bar{H}_{ij}=\partial^{2}\bar{\mathcal{H}}$, $I$ is the ${3\times 3}$ unit matrix and $T_{i}$ is a ${3\times 2}$ matrix that transforms into a tangent space basis of spin $i$. As the Hessian matrix (23) is represented in the $2N$-dimensional tangent basis, the evaluation of an eigenmode in the $3N$-representation of the embedding space $\mathcal{E}$ requires a transformation back, i.e. $\lambda|^{3N}=T\lambda|^{2N}$. Further details on the above mathematical concepts and notations can be found in Appendix H.

For a single spin, the energy landscape and force vectors can be visualized easily as the phase space is two-dimensional. An illustration of the method is shown in Fig. 11 for a system consisting of one movable spin interacting with a second, pinned spin. The parameters of the Hamiltonian are, relative to the exchange constant $J$,

[TABLE]

where the anisotropy $K$ is used to reduce the symmetry of the energy landscape. The figure illustrates how the minimum mode can be used to invert the right part of the gradient force in order to obtain a force that directs the system to a first order saddle point.

The test of a larger and far more complex system has been given in Ref. mueller_duplication_2018, , where the minimum mode following method revealed the existence of a skyrmion duplication mechanism. By defining the force field in the above way, previously unknown transition mechanisms can be found and subsequently used in the calculation of lifetimes. Applying this saddle point search method to three-dimensional systems will likely identify an even larger variety of mechanisms, as the additional dimension can significantly increase the amount of possible transitions.

IV Conclusions

The functionality of a comprehensive simulation framework, Spirit, for studies of atomic scale magnetic systems is presented and various example applications described. It is an open source software written in the C$++$ programming language and is available for free under the so-called MIT license (see Ref. spirit, ). Spirit is a very flexible, high-performance, and interactive tool, able to simulate for example ferromagnets, antiferromagnets, synthetic antiferromagnets, ferrimagnets, noncollinear magnetic structures, vortices or skyrmions. Arbitrary geometries and interactions can be described, such as bulk systems, thin films, exchange bias, multilayers, nanotubes or core-shell nanoparticles. The computational domain can be treated by open and periodic boundary conditions and can be subjected to external magnetic fields, temperature and spin-current induced torques. Due to the fact that it can be used with the Python programming language, Spirit can integrate perfectly into multiscale simulations and workflow automation frameworks, such as ASE ASE or AiiDA. AiiDA It can be used on most common architectures, such as desktop and laptop computers, clusters or supercomputers and even current day mobile devices. The calculations can be parallelized both on CPUs and GPUs.

Various simulation methods have been implemented, including Monte Carlo, Landau-Lifshitz-Gilbert dynamics, Langevin dynamics, geodesic nudged elastic band and minimum mode following methods as well as the calculations of transition rates and lifetimes within the harmonic approximation to transition-state theory. The basic algorithms of these methods have been outlined, their implementation verified and applications to several systems, such as vortices, domain walls, skyrmions and boobers are described. The parameters of the simulation can be set and modified in real time through a graphical user interface and the output of the simulations can be visualized easily.

We note that a micromagnetic description of the energetics could easily be implemented in Spirit and the micromagnetic calculations would then be able to make use of the various simulation methods and visualization features.

Acknowledgements.

The authors acknowledge helpful discussions with Pavel Bessarab, Stephan von Malottki, Jan Müller, Jonathan Chico, Filipp N. Rybakov, and Florian Rhiem. G.P.M. acknowledges funding by the Icelandic Research Fund (grants 185405-051 & 184949-051) and M.H. and S.B. from MAGicSky Horizon 2020 European Research FET Open project (#665095) and from the DARPA TEE program through grant MIPR (#HR0011831554) from DOI.

Appendix A Determination of topological charge for spin density on a lattice

For the proper definition of the topological charge of a discrete lattice of spins $\bm{n}(x_{i},y_{i})$, where $i$ runs over all lattice sites, we follow the definition given by Berg and Lüscher,Berg81 and arrive at the following expression:

[TABLE]

with

[TABLE]

where $l$ runs over all elementary triangles of the hexagonal lattice, and $A_{l}$ is the solid angle, i.e. the area of the spherical triangle with vertices $\bm{n}_{i}$, $\bm{n}_{j}$, $\bm{n}_{k}$, see Fig. 12. The sign of $A_{l}$ is determined as $\operatorname{sign}\left(A_{l}\right)=\operatorname{sign}\left[\bm{n}_{i}\cdot\left(\bm{n}_{j}\times\bm{n}_{k}\right)\right]$.

The sites $i$, $j$, $k$ of each elementary triangle are numbered in a counter-clockwise sense relative to the surface normal vector $\hat{\bm{r}}_{\perp}$ pointing in positive direction of the $z$-axis. The latter means that the numbering should satisfy the condition $\hat{\bm{r}}_{\perp}\cdot(\bm{r}_{ij}\times\bm{r}_{ik})>0$, where $\bm{r}_{ij}$ is a connection vector directed from lattice site $i$ to $j$.

The parameter $A_{l}$ can be thought of as local topological charge, which takes values in the range of $-2\pi<A_{l}<+2\pi$. According to Berg and Lüscher,Berg81 there is a set of exceptional spin configurations for which $Q$ is not defined but still measurable as $A_{l}$ in (26) is defined for all possible spin configuration. The exceptional spin configurations correspond to the case when a spherical triangle degenerates to a great circle $A_{l}=2\pi$. In this case the orientation of $A_{l}$ becomes ambiguous and the position of these elementary triangles $l^{*}$ are considered as exceptional configurations or topological defects of a two-dimensional magnetic structure.

These topological defects satisfy the following condition:

[TABLE]

The elementary triangle $l^{*}$, for which the condition (27) is satisfied, can be considered as the position at which the localization of a topological defect takes place. It is important to note that the definition of the topological charge given above remains correct only for spatially extended two-dimensional systems. This means that a topological analysis of the spin structure on a finite size domain is only defined if periodical boundary conditions are present. In the case of open boundary conditions, strictly speaking, the topological charge is not defined.

Appendix B Heun’s solver

To simplify the following discussion, we write the LLG equation (10) as

[TABLE]

where $\{\bm{n}_{j}\}$ is the set of all spins and we keep the explicit time-dependence of $\bm{A}_{i}$, as the Hamiltonian can be time-dependent, for example when an AC magnetic field is used. Heun’s method is a common and illustrative way to solve ordinary differential equations (ODEs) by first calculating an intermediate prediction step and then “averaging” to obtain the final approximation. Denoting the time step $\delta t$, for an ODE of the form

[TABLE]

the predicted value $y^{p}$ is first calculated as

[TABLE]

and then the approximation for the next step as

[TABLE]

When applied to the LLG equation, where $f\widehat{=}\bm{n}\times\bm{A}$, this integration scheme obviously does not intrinsically preserve the spin length, requiring the re-normalization of the vectors $\bm{n}_{i}$ after a given number of iterations, depending on the required precision. Note that Heun’s method falls into the category of Runge-Kutta methods, which function analogously and therefore all have this property.

In order to improve on this, Ref. depondt_spin_2009, proposes to make use of the fact that the spins are only allowed to rotate, by writing an appropriate rotation matrix $R_{i}$, which is calculated directly from the field $\bm{A}_{i}$. Applied to Heun’s method, the prediction step (30) reads

[TABLE]

To perform the correction step (31), one needs the correction field $\bm{A}^{c}$, which is calculated from the average of the initial and predicted fields:

[TABLE]

From this, in turn, the rotation matrix for the correction step $R_{i}^{c}\left(\bm{A}_{i}^{c}\right)$ is obtained and the final step of the scheme reads

[TABLE]

Higher order Runge-Kutta schemes could apply this approach analogously.

Appendix C Semi-implicit midpoint solver

Instead of using a Runge-Kutta type scheme, as described in Appendix B, Ref. mentink_stable_2010, takes a different approach, using the implicit midpoint (IMP) structure to preserve the spin length. The corresponding prediction step in its implicit form reads

[TABLE]

from which the corrector step can be analogously calculated as

[TABLE]

The solutions to these equations can be obtained analytically by rewriting them in a skew matrix form and applying Cramer’s rule (see Appendix C).

The implicit midpoint method, which the SIB method bases on, solves differential equations of the form $y^{\prime}(t)=f(t,y(t))$,   $y(t_{0})=y_{0}$ (see Eq. (8) of the main text) and an iteration step is defined as

[TABLE]

For the LLG equation (28) and a time step $\delta t$ this leads us to

[TABLE]

The semi-implicit scheme B (SIB) mentink_stable_2010 uses a predictor $\bm{n}^{p}_{i}$ to reduce the implicitness of the equation above by removing the dependence of $\bm{A}_{i}$ on $\bm{n}_{j}(t+\delta t)$. To preserve the spin length the predictor is obtained with the IMP structure.

[TABLE]

Eq. (39) can be rewritten as:

[TABLE]

with the matrix

[TABLE]

The right hand side of Eq. (40) can be easily calculated as:

[TABLE]

To solve Eq. (40) we use Cramer’s rule. The components $n^{p}_{i,\alpha}$ with $\alpha=x,y,z$ of $\bm{n}^{p}_{i}$ are calculated with

[TABLE]

where $M^{\alpha}$ is the same matrix as $M$ but column $\alpha$ is replaced with the vector $\bm{a}$, for example

[TABLE]

We now use the predictor $\bm{n}^{p}_{i}$ in the IMP step (38) to calculate $\bm{n}_{i}(t+\delta t)$:

[TABLE]

The correction step is analogous to the prediction step (compare eqs. (39) and (45)), meaning that the scheme (43) can be applied to obtain $\bm{n}_{i}(t+\delta t)$, too.

Appendix D Velocity projection solver

This description is derived from Ref. bessarab_method_2015, . Verlet-like methods generally find application in solving second order differential equations of the form $\ddot{x}(t)=F(t,x(t))$, $x(t_{0})=x_{0}$, $\dot{x}(t_{0})=v_{0}$, such as Newtons equation of motion. One formulation of this method is to increment both the position and the velocity at each time step

[TABLE]

[TABLE]

The velocity projection is used to accelerate convergence towards local minima and to avoid overstepping due to momentum. The velocity at each time step is damped by projecting it on the force

[TABLE]

Note that the dot product and norm in this equation denote those of $3N$-dimensional vectors.

To apply this scheme to the energy minimization of a spin system, we therefore no longer solve the LLG equation, but instead pretend that the spins are massive particles moving on the surfaces of spheres. The force is then simply

[TABLE]

As the method does not conserve the length of the spins, they should be renormalized after each iteration

[TABLE]

Note that this scheme, too, would most likely benefit from the usage of rotations instead of displacements.

Appendix E Stochastic LLG

Instead of Monte Carlo, one can also sample the stochastic LLG equation over time. We present here the results of such sampling for the same system and parameters, as the example shown in Fig. 4. Recall the expected critical temperature $T_{c}=1.44~{}J/k_{\mathrm{B}}\approx 16.71$ K. Fig. 13 shows the results.

The results shown in Fig. 13 demonstrate the validity of the implementation, as the expected critical temperature of $T_{\mathrm{C}}\approx 16.71$ K is matched with an error of only $1\%$. Note, however, the higher number of samples (compared to Monte Carlo) required to obtain this result: at each temperature $200$k thermalisation steps were made before taking $1$M samples.

Appendix F GNEB tangents and forces

For spin systems, special care has to be taken due to the fact that the phase space is curved (the spins are restricted to unit spheres (see also Appendix H)). The expression for $l_{\nu,\mu}$ should not be the Euclidean distance norm, but the geodesic (here, the great-circle) distance. Further, the tangents $\tau_{\nu}$ need to lie in the tangent space to their corresponding image. One may correct the tangents for example by a simple projection, orthogonalizing the corresponding $3$-component subvectors with respect to the spins

[TABLE]

After this, the tangent needs to be re-normalized $\tau_{\nu}\to\tau_{\nu}/|\tau_{\nu}|$. This tangent projection is illustrated for a single spin in Fig. 14. As the spring forces are constructed from tangent vectors, they are by definition in the tangent space. Finally, for the energy gradient force, the same scheme as for the tangents can be applied and we write for each spin

[TABLE]

Appendix G GNEB Parameters of the single spin system

The energy surface of the single-spin system, shown in Fig. 7 in the main text to illustrate the geodesic nudged elastic band method is defined for a single spin as a sum of Gaussians of the form

[TABLE]

with parameters given in Table 1.

Appendix H Details on the curved manifold

The following has been detailed in the supplementary material of Ref. mueller_duplication_2018, , but the key ideas are reproduced here. Both the HTST and MMF methods require the calculation and diagonalization of the Hessian matrix. However, when treating Riemannian manifolds, the second derivatives do not have an intrinsic geometrical meaning and therefore need to be treated with special care. Nakahara2003 In a spin system where the spin length is fixed, the manifold $\mathcal{M}_{\mathrm{phys}}$ of physical states is composed of the direct product of $N$ spheres

[TABLE]

Hence, $\mathcal{M}_{\mathrm{phys}}$ is a submanifold of the embedding euclidean space $\mathcal{E}=\mathbb{R}^{3N}$.

It turns out to be convenient to treat the spins and derivatives with respect to their orientations in a $3N$-dimensional cartesian representation. This also avoids problems of other representations, such as the singularities which arise at the poles of spherical coordinates. The derivatives in the embedding space $\mathcal{E}$ are readily calculated by extending the Hamiltonian $\mathcal{H}$, which is defined on $\mathcal{M}_{\mathrm{phys}}$ to a function $\bar{\mathcal{H}}$ on $\mathcal{E}$. While we denote the gradient taken in the embedding space $\mathcal{E}$ as $\partial\bar{\mathcal{H}}$, the gradient taken on the manifold $\mathcal{M}_{\mathrm{phys}}$ has to lie in the tangent space to the manifold, which we write as a projection $P_{x}\partial\bar{\mathcal{H}}$. The Hessian matrix of second derivatives in the embedding space $\mathcal{E}$ is denoted $\partial^{2}\bar{\mathcal{H}}$.

In this extrinsic view onto the spin manifold, the covariant second derivatives can be extracted from a projector approach, absil_extrinsic_2013 where for any scalar function $f$ on the manifold $\mathcal{M}_{\mathrm{phys}}$, the covariant Hessian is defined as

[TABLE]

$W_{x}$ denotes the Weingarten map, which, for a spherical manifold, for any vector $v$ at a point $x$ is given by

[TABLE]

where $z$ is a tangent vector to the sphere at $x$. To calculate the Hessian, we insert $v=P_{x}^{\perp}\partial\bar{\mathcal{H}}$ and retrieve

[TABLE]

where $x^{T}\partial\bar{\mathcal{H}}$ is the scalar product of the spin with the gradient.

To illustrate the implementation in Spirit, spirit we switch notation to matrix representation and drop the subscript $x$. For spin indices $i$ and $j$, the gradient $\partial\bar{\mathcal{H}}$ can be written as a $3$-dimensional object $\nabla_{i}\bar{\mathcal{H}}$ and the second derivative $\partial^{2}\bar{\mathcal{H}}$ as a matrix $\bar{H}$. In Euclidean representation, the Hessian of Eq. (55) becomes as a ${3N\times 3N}$ matrix

[TABLE]

consisting of $N^{2}$ blocks, each corresponding to a different spin-spin subspace. It is obtained by acting with Eq. (55) on the euclidean basis vectors of the embedding space $\mathcal{E}$. These subspace matrices of size ${3\times 3}$ are given by

[TABLE]

where $I$ denotes the ${3\times 3}$ unit matrix.

The matrix $H|^{3N}$ of course describes $3N$ degrees of freedom, while there can only be $2N$ physical eigenmodes of the spins, spanning the tangent space to the spin configuration. In order to remove the unphysical degrees of freedom in the embedding space $\mathcal{E}$, is is sufficient to transform the matrix into a tangent space basis, which we can write as $H_{ij}=T_{i}^{T}H_{ij}|^{3N}T_{j}$, where $T_{i}$ is the basis transformation matrix of spin $i$ fulfilling $T^{T}P=T^{T}$ and $T^{T}T=I|^{2N}$. The true Hessian $H=(H_{ij})$ of Eq. (55) in the ${2N\times 2N}$ matrix representation, containing only the physical degrees of freedom, is therefore defined as

[TABLE]

Note that this reduction of dimensionality also improves the numerical efficiency of the diagonalization. As the eigenmodes $\lambda|^{2N}$ are represented in the tangent basis, the $3N$ representation needs to be calculated by $\lambda|^{3N}=T\lambda|^{2N}$.

While the ${3\times 2}$ basis matrix $T_{i}$ can be calculated quite arbitrarily by choice of two orthonormal vectors, tangent to the spin $\bm{n}_{i}$, we found it convenient to use the unit vectors of spherical coordinates $\theta$ and $\varphi$

[TABLE]

where $r_{xy}=\sin\theta=\sqrt{1-z^{2}}$. Note that the poles need to be excluded, but since the basis does not need to be continuous over the manifold, one may e.g. orthogonalize $\bm{e}_{x}$ and $\bm{e}_{y}$ with respect to the spin vector to obtain suitable tangent vectors.

Finally, the Hessian matrix in the embedding space $\mathcal{E}=\mathbb{R}^{3N}$ is needed, denoted $\bar{H}_{ij}|^{3N}$. As the atomistic Hamiltonian can generally be written in matrix form

[TABLE]

where $A_{ij}$ are matrices of size ${3\times 3}$ describing the linear contributions, such as the Zeeman term, and $B_{ij}$ are matrices describing the quadratic contributions, such as anisotropy, Exchange, DMI and dipolar interactions. The Hessian matrix is then naturally given by

[TABLE]

Appendix I Minimum modes in the interacting spin system

The following Fig. 15 illustrates how the minimum eigenmode unit vector $\hat{\lambda}$ is oriented and in which direction, therefore, the gradient force is inverted. Recall Eq. (22), which can be written

[TABLE]