Spirit: Multifunctional Framework for Atomistic Spin Simulations
The 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.
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 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, .
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 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
In Spirit, we implemented an extended Heisenberg Hamiltonian aharoni_introduction_2000 (); rado_magnetism_1963 () of classical spins of unit length located at lattice sites giving rise to the magnetic moment . The general form
includes (i) the Zeeman term describing the interaction of the spins with the external magnetic field , (ii) the single-ion magnetic anisotropy, where are the axes of the uniaxial anisotropies of the basis cell with the anisotropy strength , (iii) the symmetric exchange interaction given by and the antisymmetric exchange, also called Dzyaloshinskii-Moriya interaction, given by vectors , where denotes the unique pairs of interacting spins and , (iv) the dipolar interactions, where 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
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 , where 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 . 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 which is composed of the upper edges of the cube:
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 with increasing resolution of the grid, where a cube with lattice sites already gives an agreement within and the shown example with sites a discretization error of less than with respect to the continuum solution.
Topological charge. In the past years, we witnessed the characterization of smooth magnetization fields in terms of topological classes. In this case, the winding number defined as
It should be noted that a secondary topological charge can be considered as the defining index to distinguish between skyrmion () and antiskyrmion () independently of the background state. It is defined as
and is also referred to as the winding number in comparison to the conventional topological charge . Here, 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 . For sufficiently regular fields decaying to the background state , the index is defined relative to this background state, .hoffmann_antiskyrmions_2017 () In a typical situation where the horizontal magnetization field vanishes at a single point (skyrmion center) the relative skyrmion number agrees with the index of the horizontal magnetization field at the skyrmion center. It is customary to fix the background state , which leads to the characterization for skyrmions and for antiskyrmions.
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 -axis and choosing a new spin direction by uniformly distributed random variables and , where is the opening angle of the cone. The trial step is accepted with a probability
where 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 ()
where is the average magnetization of the sample, while the specific heat relates to fluctuations of the energy
where both should diverge at the critical temperature for a phase transition, e.g. to the paramagnetic phase. The order Binder cumulant binder_finite_1981 (), which is often used to avoid finite size scaling effects, is defined as
Fig. 4 shows these quantities as results of a Monte Carlo calculation of a cube of lattice sites for 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 () 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 .
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
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. is the electron gyromagnetic ratio, is the damping parameter, is the effective field, is a non-adiabatic parameter, is the electron current vector, and is the spatial gradient acting on the spin orientation. The effective field always contains a component related to the energy gradient , but in this notation for the LLG equation, the effective field may contain also a stochastic thermal field, i.e. , given by
where the magnitude is given by the fluctuation-dissipation theorem and is white noise, such that and . To achieve these properties in an implementation, the vectors 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 . 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.
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 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
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 . The chain is oriented along the -axis and the first and the last spin are fixed in and direction, respectively. As a subset of the general Hamiltonian (1), the Hamiltonian for this example can be written as follows:
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 over applied current in normalised units.
The approximate prediction thiaville_domain_2004 () 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 , which is in close agreement with the reported value of . Note, for and currents larger than and for and currents larger than , the wall starts rotating around the -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
where is the image index along the chain, is a spring force, and is an energy gradient force. The forces in this section are -dimensional vectors. A simple definition of the spring force, which gives an equidistant distribution of images in phase space, is given by
where is a measure of distance between images and and is the (normalized) path tangent at image . The 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
where . 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:
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 and , meaning that the incommensurate spin spiral has a period of . 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 , while Ref. rybakov_new_2015, – not using climbing images – reports a ratio of .
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 , 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 and the energy barrier of the transition :
where the M and S superscripts indicate the minimum and first order saddle point of the transition. The are eigenvalues of the Hessian matrix (see Appendix H), are the phase space volumes of zero modes (if present, otherwise ), are the number of zero modes – modes with zero eigenvalue – and 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 are in fact velocities: the first row of the dynamical matrix transformed into the eigenbasis of the Hessian according to
where is a basis matrix of the tangent space and denote the matrix of the Hessians eigenvectors (in -representation, i.e. in the basis ). 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 of an isolated skyrmion in a two-dimensional system. As parameters we chose and 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 T and T and temperature between K and 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
where is the negative gradient of the energy and 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 -dimensional for a system with spins.
The calculation of second derivatives requires further attention, as the requirement of constant length effectively constrains the spins to a sub-manifold of an embedding space . 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 Hessian matrix can be represented as
where and are spin indices, is the smooth continuation of the Hamiltonian to the embedding space, , is the unit matrix and is a matrix that transforms into a tangent space basis of spin . As the Hessian matrix (23) is represented in the -dimensional tangent basis, the evaluation of an eigenmode in the -representation of the embedding space requires a transformation back, i.e. . 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 ,
where the anisotropy 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.
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 , where runs over all lattice sites, we follow the definition given by Berg and Lüscher,Berg81 () and arrive at the following expression:
where runs over all elementary triangles of the hexagonal lattice, and is the solid angle, i.e. the area of the spherical triangle with vertices , , , see Fig. 12. The sign of is determined as .
The sites , , of each elementary triangle are numbered in a counter-clockwise sense relative to the surface normal vector pointing in positive direction of the -axis. The latter means that the numbering should satisfy the condition , where is a connection vector directed from lattice site to .
The parameter can be thought of as local topological charge, which takes values in the range of . According to Berg and Lüscher,Berg81 () there is a set of exceptional spin configurations for which is not defined but still measurable as 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 . In this case the orientation of becomes ambiguous and the position of these elementary triangles are considered as exceptional configurations or topological defects of a two-dimensional magnetic structure.
These topological defects satisfy the following condition:
The elementary triangle , 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
where is the set of all spins and we keep the explicit time-dependence of , 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 , for an ODE of the form
the predicted value is first calculated as
and then the approximation for the next step as
When applied to the LLG equation, where , this integration scheme obviously does not intrinsically preserve the spin length, requiring the re-normalization of the vectors 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 , which is calculated directly from the field . Applied to Heun’s method, the prediction step (30) reads
To perform the correction step (31), one needs the correction field , which is calculated from the average of the initial and predicted fields:
From this, in turn, the rotation matrix for the correction step is obtained and the final step of the scheme reads
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
from which the corrector step can be analogously calculated as
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 , (see Eq. (8) of the main text) and an iteration step is defined as
For the LLG equation (28) and a time step this leads us to
The semi-implicit scheme B (SIB) mentink_stable_2010 () uses a predictor to reduce the implicitness of the equation above by removing the dependence of on . To preserve the spin length the predictor is obtained with the IMP structure.
Eq. (39) can be rewritten as:
with the matrix
The right hand side of Eq. (40) can be easily calculated as:
To solve Eq. (40) we use Cramer’s rule. The components with of are calculated with
where is the same matrix as but column is replaced with the vector , for example
We now use the predictor in the IMP step (38) to calculate :
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 , , , such as Newtons equation of motion. One formulation of this method is to increment both the position and the velocity at each time step
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
Note that the dot product and norm in this equation denote those of -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
As the method does not conserve the length of the spins, they should be renormalized after each iteration
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 K. Fig. 13 shows the results.
The results shown in Fig. 13 demonstrate the validity of the implementation, as the expected critical temperature of K is matched with an error of only . Note, however, the higher number of samples (compared to Monte Carlo) required to obtain this result: at each temperature k thermalisation steps were made before taking 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 should not be the Euclidean distance norm, but the geodesic (here, the great-circle) distance. Further, the tangents 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 -component subvectors with respect to the spins
After this, the tangent needs to be re-normalized . 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
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
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 of physical states is composed of the direct product of spheres
Hence, is a submanifold of the embedding euclidean space .
It turns out to be convenient to treat the spins and derivatives with respect to their orientations in a -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 are readily calculated by extending the Hamiltonian , which is defined on to a function on . While we denote the gradient taken in the embedding space as , the gradient taken on the manifold has to lie in the tangent space to the manifold, which we write as a projection . The Hessian matrix of second derivatives in the embedding space is denoted .
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 on the manifold , the covariant Hessian is defined as
denotes the Weingarten map, which, for a spherical manifold, for any vector at a point is given by
where is a tangent vector to the sphere at . To calculate the Hessian, we insert and retrieve
where 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 . For spin indices and , the gradient can be written as a -dimensional object and the second derivative as a matrix . In Euclidean representation, the Hessian of Eq. (55) becomes as a matrix
consisting of 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 . These subspace matrices of size are given by
where denotes the unit matrix.
The matrix of course describes degrees of freedom, while there can only be 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 , is is sufficient to transform the matrix into a tangent space basis, which we can write as , where is the basis transformation matrix of spin fulfilling and . The true Hessian of Eq. (55) in the matrix representation, containing only the physical degrees of freedom, is therefore defined as
Note that this reduction of dimensionality also improves the numerical efficiency of the diagonalization. As the eigenmodes are represented in the tangent basis, the representation needs to be calculated by .
While the basis matrix can be calculated quite arbitrarily by choice of two orthonormal vectors, tangent to the spin , we found it convenient to use the unit vectors of spherical coordinates and
where . 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 and with respect to the spin vector to obtain suitable tangent vectors.
Finally, the Hessian matrix in the embedding space is needed, denoted . As the atomistic Hamiltonian can generally be written in matrix form
where are matrices of size describing the linear contributions, such as the Zeeman term, and are matrices describing the quadratic contributions, such as anisotropy, Exchange, DMI and dipolar interactions. The Hessian matrix is then naturally given by
Appendix I Minimum modes in the interacting spin system
- (1) M. Hoffmann, B. Zimmermann, G. P. Müller, D. Schürhoff, N. S. Kiselev, C. Melcher, and S. Blügel, "Antiskyrmions Stabilized at Interfaces by Anisotropic Dzyaloshinskii-Moriya Interactions," Nat. Commun. 8, 308 (2017).
- (2) I. Žutić, J. Fabian, and S. Das Sarma. "Spintronics: Fundamentals and applications," Rev. Mod. Phys. 76 (2004).
- (3) S. D. Bader. "Colloquium: Opportunities in nanomagnetism," Rev. Mod. Phys. 78 (2006).
- (4) W. F. Brown, "Micromagnetics," Interscience Publishers (1963).
- (5) Spirit – spin simulation framework (see https://spirit-code.github.io).
- (6) M.J. Donahue and D.G. Porter, "OOMMF User’s Guide, Version 1.0," Interagency Report NISTIR 6376, National Institute of Standards and Technology, Gaithersburg, MD (1999).
- (7) A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, "The design and verification of MuMax3," AIP Advances 4 107133 (2014).
- (8) B. Skubic, J. Hellsvik, L. Nordström, and O. Eriksson, "A Method for Atomistic Spin Dynamics Simulations: Implementation and Examples," J. Phys.: Condens. Matter 20 315203 (2008).
- (9) R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. Ellis, and R. W. Chantrell, "Atomistic Spin Model Simulations of Magnetic Nanomaterials," J. Phys. Cond. Mat. 26 10, 0953-8984 (2014).
- (10) F. N. Rybakov and N. S. Kiselev, "Chiral Magnetic Skyrmions with Arbitrary Topological Charge ("skyrmionic sacks")," arXiv:1806.00782 (2018).
- (11) U. Nowak, "Thermally Activated Reversal in Magnetic Nanostructures," Annual Reviews of Computational Physics IX, 105 (2001).
- (12) P. F. Bessarab, V. M. Uzdin, and H. Jónsson, "Method for Finding Mechanism and Activation Energy of Magnetic Transitions, Applied to Skyrmion and Antivortex Annihilation," Comp. Phys. Comm. 196 335 (2015).
- (13) G. P. Müller, P. F. Bessarab, S. M. Vlasov, F. R. Lux, N. S. Kiselev, S. Blügel, V. M. Uzdin, and H. Jónsson, "Duplication, Collapse, and Escape of Magnetic Skyrmions Revealed Using a Systematic Saddle Point Search Method," Phys. Rev. Lett. 121 19 197202 (2018).
- (14) P. F. Bessarab, V. M. Uzdin, and H. Jónsson, "Harmonic Transition-State Theory of Thermal Spin Transitions," Phys. Rev. B 85 184409 (2012).
- (15) H.-B. Braun, "Topological effects in nanomagnetism: from superparamagnetism to chiral quantum soliton," Advances in Physics 61 1 (2012).
- (16) VFRendering – A vector field rendering library (see https://github.com/FlorianRhiem/VFRendering).
- (17) Y. Liu, R. Lake, and J. Zang, "Binding a Hopfion in Chiral Magnet Nanodisk," arxiv:1806.01682 (2018).
- (18) F. Zheng, F. N. Rybakov, A. B. Borisov, D. Song, S. Wang, Z.-A. Li, H. Du, N. S. Kiselev, J. Caron, A. Kovács, M. Tian, Y. Zhang, S. Blügel, and R. E. Dunin-Borkowski, "Experimental Observation of Chiral Magnetic Bobbers in B20-Type FeGe," Nat. Nanotechnol. 13, 451-455 (2018).
- (19) H. Du, X. Zhao, F. N. Rybakov, A. B. Borisov, S. Wang, J. Tang, C. Jin, C. Wang, W. Wei, N. S. Kiselev, Y. Zhang, R. Che, S. Blügel, and M. Tian, "Interaction of Individual Skyrmions in a Nanostructured Cubic Chiral Magnet," Phys. Rev. Lett. 120, 197203 (2018).
- (20) J. Hagemeister, A. Siemens, L. Rózsa, E. Y. Vedmedenko, and R. Wiesendanger, "Controlled Creation and Stability of Skyrmions on a Discrete Lattice," Phys. Rev. B 97, 174436 (2018).
- (21) M. Redies, F. R. Lux, P. M. Buhl, G. P. Müller, N. S. Kiselev, S. Blügel, and Y. Mokrousov, "Distinct Magnetotransport and Orbital Fingerprints of Chiral Bobbers," arXiv:1811.01584 (2018).
- (22) Web interface for Spirit (see https://juspin.de).
- (23) A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, "The Atomic Simulation Environment—A Python library for working with atoms," J. Phys.: Condens. Matter 29 273002 (2017). (see https://wiki.fysik.dtu.dk/ase)
- (24) G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzaria, and B. Kozinsky, "AiiDA: automated interactive infrastructure and database for computational science," Comp. Mat. Sci. 111 218 (2016). (see also http://www.aiida.net/)
- (25) A. Aharoni, "Introduction to the Theory of Ferromagnetism," Oxford University Press (2000).
- (26) G. T. Rado, Magnetism: a treatise on modern theory and materials. 3. Spin arrangements and crystal structure, domains, and micromagnetics. Academic Press (1963).
- (27) M. Hoffmann and S. Blügel, "Systematic derivation of realistic spin-models for beyond-Heisenberg solids from microscopic model," arXiv:1803.01315 (2018).
- (28) A. Szilva, M. Costa, A. Bergman, L. Szunyogh, L. Nordström, and O. Eriksson, "Interatomic Exchange Interactions for Finite-Temperature Magnetism and Nonequilibrium Spin Dynamics," Phys. Rev. Lett. 111 (2013).
- (29) A. Krönlein, M. Schmitt, M. Hoffmann, J. Kemmer, N. Seubert, M. Vogt, J. Küspert, M. Böhme, B. Alonazi, J. Kügel, H. A. Albrithen, M. Bode, G. Bihlmayer, and S. Blügel, "Magnetic Ground State Stabilized by Three-Site Interactions: Fe/Rh(111)," Phys. Rev. Lett. 120, 207202 (2018).
- (30) S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel, "Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions," Nat. Phys. 7 713 (2011).
- (31) N. Hayashi, K. Saito, and Y. Nakatani, "Calculation of Demagnetizing Field Distribution Based on Fast Fourier Transform of Convolution," Japanese Journal of Applied Physics 12 (1996).
- (32) M. Frigo and S. G. Johnson, "The Design and Implementation of FFTW3," Proceedings of the IEEE 93 (2), 216 (2005). Invited paper, Special Issue on Program Generation, Optimization, and Platform Adaptation
- (33) cuFFT https://developer.nvidia.com/cufft
- (34) A. Hubert and W. Rave, "Systematic Analysis of Micromagnetic Switching Processes," Phys. Status Solidi B 211 2 815-829 (1999).
- (35) K. Binder and D. W. Heermann, "Monte Carlo Simulation in Statistical Physics," Springer-Verlag, Berlin, (1997).
- (36) D.Hinzke and U.Nowak, "Monte Carlo simulation of magnetization switching in a Heisenberg model for small ferromagnetic particles," Comp. Phys. Comm. 121 334 (1999).
- (37) D. P. Landau and K. Binder, "A guide to Monte Carlo simulations in Statistical Physics," Cambridge University Press, New York, NY. (2005).
- (38) K. Binder, "Finite Size Scaling Analysis of Ising Model Block Distribution Functions," Z. Phys. B 43, 119 (1981).
- (39) G. A. Baker, Jr., H. E. Gilbert, J. Eve, and G. S. Rushbrooke. High-temperature expansions for the spin 1/2 Heisenberg model," Phys. Rev. 164, 800 (1967).
- (40) Y. Rocio, "Modeling of Macroscopic Anisotropies Due to Surface Effects in Magnetic Thin Films and Nanoparticles," PhD thesis (2011).
- (41) R. H. Swendsen and J.-S. Wang, "Replica Monte Carlo Simulation of Spin-Glasses," Phys. Rev. Lett. 57 2607 (1986).
- (42) K. Hukushima and K. Nemoto, "Exchange Monte Carlo Method and Application to Spin Glass Simulations," J. Phys. Soc. Jap. 65 1604 (1996).
- (43) M. Böttcher, S. Heinze, S. Egorov, J. Sinova, and B. Dupé, "- Phase Diagram of Pd/Fe/Ir(111) Computed with Parallel Tempering Monte Carlo," arXiv:1707.01708 (2017).
- (44) L. D. Landau and E. M. Lifshitz, "On the Theory of the Dispersion of Magnetic Permeability in Ferromagnetic Bodies," Physik Z. Sowjetunion 8 153 (1935).
- (45) T. L. Gilbert, "A phenomenological theory of damping in ferromagnetic materials," IEEE Transactions on Magnetics, 40 (6), (2004).
- (46) W. F. Brown, "Thermal Fluctuations of a Single-Domain Particle," Physical Review 130 5, 1677–86 (1963).
- (47) C. Schieback, M. Kläui, U. Nowak, U. Rüdiger, and P. Nielaba, "Numerical Investigation of Spin-Torque Using the Heisenberg Model," The European Physical Journal B 59 4, 429–33 (2007).
- (48) B. Berg, M. Lüscher, "Definition and statistical distribution of a topological number in the Lattice -Model*," Nuclear Physics B 190[FS3] 412 (1981).
- (49) J. H. Mentink, M. V. Tretyakov, A. Fasolino, M. I. Katsnelson, and T. Rasing, "Stable and Fast Semi-Implicit Integration of the Stochastic Landau–Lifshitz Equation," J. Phys. Cond. Mat. 22 17, 176001 (2010).
- (50) D. S. G. Bauer, P. Mavropoulos, S. Lounis, and S. Blügel, "Thermally activated magnetization reversal in monatomic magnetic chains on surfaces studied by classical atomistic spin-dynamics simulations," J. Phys.: Condens. Matter 23 394204 (2011).
- (51) L. Rózsa, L. Udvardi, and L. Szunyogh, "Langevin spin dynamics based on ab initio calculations: numerical schemes and applications," J. Phys.: Condens. Matter 26 216003 (2014).
- (52) P. Depondt and F. G. Mertens, "Spin Dynamics Simulations of Two-Dimensional Clusters with Heisenberg and Dipole-Dipole Interactions," J. Phys.: Condens. Matter 21 33 (2009).
- (53) A. Thiaville, Y. Nakatani, J. Miltat, and N. Vernier, "Domain Wall Motion by Spin-Polarized Current: A Micromagnetic Study," J. Appl. Phys. 95 11 (2004).
- (54) N. L. Schryer and L. R. Walker, "The motion of domain walls in uniform dc magnetic fields," J. Appl. Phys. 45 (12):5406–5421 (1974)
- (55) G. Henkelman and H. Jónsson, "Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points," J. Chem. Phys. 113, 9978 (2000).
- (56) G. Henkelman, G. P. Uberuaga, and H. Jónsson, "A climbing image nudged elastic band method for finding saddle points and minimum energy paths," J. Chem. Phys. 113, 9901 (2000).
- (57) F. N. Rybakov, A. B. Borisov, S. Blügel, and N. S. Kiselev, "New type of particlelike state in chiral magnets," Phys. Rev. Lett. 115 117201 (2015).
- (58) E. Wigner, "The Transition State Method," Transactions of the Faraday Society 34, 29 (1938).
- (59) D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, "Current Status of Transition-State Theory," J. Phys. Chem. 100 12771 (1996).
- (60) J. S. Langer, "Statistical theory of the decay of metastable states," Annals of Physics 54 258 (1969).
- (61) P. F. Bessarab, G. P. Müller, I. S. Lobanov, F. N. Rybakov, N. S. Kiselev, Nikolai, H. Jónsson, V. M. Uzdin, S. Blügel, L. Bergqvist, and A. Delin, "Lifetime of Racetrack Skyrmions," Sci. Rep. 8 3433 (2018).
- (62) L. Desplat, D. Suess, J.-V. Kim, and R. L. Stamps, "Thermal stability of metastable magnetic skyrmions: Entropic narrowing and significance of internal eigenmodes," Phys. Rev. B 98, 134407 (2018).
- (63) S. von Malottki, P. F. Bessarab, S. Haldar, A. Delin, and S. Heinze, "Skyrmion lifetimes in ultrathin films," arXiv:1811.12067 (2018).
- (64) P.-A. Absil, R. Mahony, and J. Trumpf, "An Extrinsic Look at the Riemannian Hessian," In Geometric Science of Information 361-68, Lecture Notes in Computer Science, Springer, Berlin, Heidelberg (2013).
- (65) M. Nakahra, "Geometry, topology and physics," CRC Press (2003).