Duplication, collapse and escape of magnetic skyrmions revealed using a systematic saddle point search method
Gideon P. Müller Pavel F. Bessarab Sergei M. Vlasov Fabian R. Lux Nikolai S. Kiselev Stefan Blügel Valery M. Uzdin Hannes Jónsson
















	
 	
        

    
	

	
	

	
	

	

 	

	    
	

 	[]hj@hi.is
	
	







Various transitions that a magnetic skyrmion can undergo are found in calculations using a method for climbing up the energy surface and converging onto first order saddle points. In addition to collapse and escape through a boundary, the method identifies a transition where the skyrmion divides and forms two skyrmions. The activation energy for this duplication process can be similar to that of collapse and escape. A tilting of the external magnetic field for a certain time interval is found to induce the duplication process in a dynamical simulation. Such a process could turn out to be an important avenue for the creation of skyrmions in future magnetic devices.












Duplication, collapse and escape of magnetic skyrmions revealed using a systematic saddle point search method
    Hannes Jónsson
    December 30, 2023
=============================================================================================================







Localized, non-collinear magnetic states are receiving a great deal of attention, where skyrmions have come under special focus.
Along with interesting transport properties, skyrmions exhibit particle-like behaviour and carry a topological charge enhancing their stability with respect to uniform ferromagnetic background.
In addition to the interest in their intriguing properties, they have been suggested as a basis for technological applications e.g. data storage or even data processing devices <cit.>.
A racetrack design of a memory device has been outlined where a spin polarized current drives a chain of skyrmions past a reading device <cit.>.
The effect of temperature and external magnetic field on the stability of the skyrmions need to be studied, as well as
ways to generate and manipulate them. The effect of defects is also an important consideration <cit.>. 
Two mechanisms for the annihilation of skyrmions have been characterized by theoretical calculations of atomic scale systems: Collapse of a skyrmion to form ferromagnetic state <cit.> and escape of a skyrmion through the boundary of the magnetic domain <cit.>. 

The effect of a non-magnetic impurity has also been calculated <cit.>. 

By using harmonic transition state theory for magnetic systems <cit.>, the lifetime of skyrmions has been estimated <cit.>.
Parameter values obtained from density functional theory <cit.> are found to give
results that are consistent with experimental observations <cit.>.  
The challenge is to design materials where magnetic skyrmions are small enough while being sufficiently stable at ambient temperature, and to develop methods for manipulating them. 

Theoretical calculations can help accelerate this development 

by identifying the various possible transformations that a skyrmion can undergo at a finite temperature on a laboratory time scale. 
This can be achieved by the use of rate theory where the major challenge is to find the relevant transition mechanisms. If the final state of a transition is specified, in addition to the initial state, the geodesic nudged elastic band (GNEB) method <cit.> 
can be used to find the minimum energy path of the transition and, thereby, the activation energy which is the highest rise in energy along the path.

However, the final states of possible transitions are not always known. 

Another category of methods 
for identifying possible transition mechanisms where the final states are not specified, only the initial state, 
has turned out to be highly valuable
in a different context, namely in studies of atomic rearrangements such as chemical reactions and diffusion events 
<cit.>.  
Unexpected transition mechanisms have in many cases turned out to be preferred over mechanisms that seem a priori most likely
<cit.>.

Here, we describe a method that can be applied to identify transition mechanisms in magnetic systems without specifying final states. 
It represents an adaptation of a method that has been used extensively in studies of atomic rearrangements. 
A complication arises from the fact that magnetic systems are characterized by the 
orientation of the magnetic moments while the length of the magnetic moments is either fixed or obtained from 
self-consistency calculations <cit.>.  
The configuration space is therefore curved.  
The method is first described briefly, with more detailed information in Supplemental Material <cit.>, and then an application to a magnetic skyrmion is presented
where, in addition to collapse and escape, the method gives a mechanism and activation energy 
that has not been reported before: Duplication of a skyrmion.
Finally, a dynamical simulation is described where a time dependent external field is used to induce such an event. 





Within harmonic transition state theory <cit.>, 
the mechanism and rate of a thermally induced transition is characterized by the first order saddle point <cit.>
representing the bottleneck for the transition. 
Given an initial state corresponding to a local energy minimum, the various possible transitions the system may undergo can be identified by climbing up the energy surface and converging on the various first order saddle points on the energy ridge surrounding the minimum. 

A version of this approach, referred to as minimum mode following 
<cit.>, is illustrated for a single-spin test system in Fig. 1. 
The method is based on the evaluation of 


eigenvalues and corresponding eigenvectors of the Hessian, .
First, the region near the minimum, the convex region where all eigenvalues of the Hessian are positive, is escaped by following, for example, the gradient of the energy. Alternatively, a random vector or one of the eigenvectors of  can be followed to escape the convex region. Once an eigenvalue turns negative, an effective force 




    F^eff = F - 2 (λ̂·F) λ̂ ,


is followed, where F=-∇ℋ is the negative gradient of the energy

and λ̂ is the normalized eigenvector corresponding to the negative eigenvalue. Note that these vectors are 3N-dimensional for a system with N spins.
As the system is displaced in the direction of the effective force,



it moves to higher energy along λ̂ but to lower energy along the orthogonal degrees of freedom. 
Eventually, this brings the system to a first order saddle point on the energy surface. 
The final state of the transition can be obtained by a slight displacement further along λ̂, 
followed by energy minimization.





This approach has not, to the best of our knowledge, been applied previously to a magnetic system. 
Here, the configuration space ℳ_phys 
is given by the direct product of N spheres.
In order to apply the mode following method, knowledge of second order derivatives is required
but, as is well-known from the theory of Riemannian manifolds,  
they need
to be treated with special care.

The Hessian
can be calculated by application of covariant derivatives, but

their evaluation in spherical coordinates is usually cumbersome and suffers from singularities at the poles.

A more convenient approach 
is offered by viewing the configuration space as being embedded in a surrounding euclidean space ℰ⊃ℳ_phys. In this larger space, second order derivatives are easily performed. The Hessian in the physical subspace can then be reconstructed by a projection operator approach <cit.>. For any scalar function f on the manifold ℳ_phys, this true Hessian is defined as

    Hess f(x)[z] = P_x∂^2 f̅(x) z + W_x (z, P_x^⊥∂f̅),


where z is a vector tangent to the manifold, P_x and P_x^⊥ are the projectors onto the tangent space and onto the normal space, respectively, to the surface of the manifold at a point x on the sphere, f̅ is the smooth extension of f to the Euclidean space, and W_x is the Weingarten map at x. 
The Weingarten map, sometimes also referred to as the shape operator, describes the curvature of 2D surfaces in terms of an embedding space.



For a Hamiltonian ℋ of a spin system, the matrix representation of the Hessian is obtained by its action on the basis vectors

(see <cit.>), i.e.,

    _ij = T_i^T _ij T_j - T_i^T I (n⃗^j·∇^jℋ̅) T_j,


where the indices i and j denote spins, =∂^2 ℋ̅(x), I is the 3×3 unit matrix and T_i is the 3×2 matrix that transforms into the tangent space of spin i.



As the Hessian matrix given by Eq. (<ref>) is represented in 2N, 

the evaluation of an eigenmode in the 3N-representation
requires a


transformation, i.e. λ|^3N = Tλ|^2N (see <cit.>).




The formulation of the constrained Hessian given by Eq. (<ref>) and the saddle point search method described above have been implemented in the Spirit <cit.> software and used here to analyse transitions from a magnetic skyrmion state.  


The energy of the system is described by an extended Heisenberg model


    ℋ = -μ_S∑_i=1^N
    	
    
    	H⃗·n⃗_i - J ∑_⟨ij|⟩ n⃗_i·n⃗_j - D ∑_⟨ij|⟩d̂⃗̂_ij· (n⃗_i×n⃗_j),


where H⃗ is a uniform external magnetic field, n⃗_i is the magnetic moment of 
spin i, J is the exchange coupling between nearest neighbor spins, 

the parameter D 
and unit vector d̂⃗̂_ij give the Dzyaloshinskii-Moriya vector in the plane of the lattice parallel to the vector connecting the two nearest neighbors i and j.
The sums include only distinct nearest neighbor pairs. 
The system consists of 40×40 spins on a square lattice with free boundary conditions and
the parameter values 
and field strength
are chosen to be the same as in a recent
theoretical study <cit.> (values are given in <cit.>).
The Bloch skyrmion, shown in Fig. <ref>, is metastable with respect to the ferromagnetic phase.






Fig. <ref> shows an illustration of the eigenvectors corresponding to the three lowest eigenvalues. They correspond to translation, breathing and elliptical distortion

 <cit.>. 
In order to escape from the convex region, we have chosen here to follow each one of these three modes until the corresponding eigenvalue becomes negative. After that, the system is displaced in the direction of the effective force given by 
Eq. (<ref>), where λ̂ is the eigenvector of the initially selected mode rather than the one with lowest eigenvalue.












By following the translational mode, the mode that has lowest eigenvalue initially, the skyrmion moves towards the edge of the system and eventually, as it is pushed along the effective force, converges on a saddle point corresponding to an escape through the edge. Relaxation of the system after a slight displacement along the unstable mode at the saddle point brings the system to the ferromagnetic state.
Similarly, by following the breathing mode, the second lowest eigenvalue mode, the skyrmion shrinks and eventually, 

converges to a saddle point corresponding to collapse. Again the final state of the transition is the ferromagnetic state.  
However, by following the elliptical mode, the skyrmion becomes stretched and 

converges on a saddle point corresponding to a division of the skyrmion and formation of two skyrmions. While the initial state contains one skyrmion, the final state contains two. 

The energy along the minimum energy paths for the three transitions is shown in Fig. <ref>. These were calculated using the geodesic nudged elastic band method <cit.> using an initial path formed by linear interpolation between the intial state and the saddle point, as well as between the saddle point and the final state. The activation energy for the collapse and escape is quite similar, but the activation energy for the duplication is higher 

in this case.
The relative height of the energy barriers for the three processes depends on the parameter values. 
A calculation using a 10% smaller field and 60% larger value of D/J gives a lower activation energy for duplication than the other two transitions <cit.>.

The shape of the three curves is quite different. The total displacement along the paths
from the initial state to the saddle point is shortest for the duplication, while the displacement of the skyrmion to the boundary and to the saddle point for escape involves the largest displacement. The energy profile for the collapse is similar to what has been presented before, a gradual shrinkage of the skyrmion to the saddle point followed by abrupt drop of the energy to that of the ferromagnetic state  <cit.>.

The duplication of the skyrmion leads to an increase in the energy of the system because the skyrmion is only metastable with respect to the ferromagnetic state for this set of parameters, while the other two transitions lead to a decrease in the energy.
The energy along the minimum energy path past the saddle point for duplication contains important information about skyrmion interaction. It shows how the repulsive interaction between skyrmions varies with the distance between them. 

The minimum energy path going from right to left along the mininum energy path in Fig. 2 shows how two skyrmions can merge to form one skyrmion.




The duplication transition identified here has not been described previously, but 
may turn out to be an important mechanism for generating skyrmions.

The first order saddle point on the energy surface corresponding to this mechanism can be found for a wide range of parameter
values, as will be described in detail in a later publication.
For example, we have found the duplication saddle point for a Pd/Fe bilayer on an Ir(111) substrate, a system that has been studied
extensively, using parameter values that give close agreement with experimental results <cit.>.
There, the energy barrier for duplication turns out to be slightly lower than that of collapse, 78 meV vs. 80 meV.


The question now arises whether it is possible to induce the duplication event dynamically.
We address this by performing a dynamical simulation where a short magnetic pulse is applied in addition to the stationary magnetic
field. 
The pulse is represented by a uniform magnetic field which has both in-plane and out-of-plane components. 
For the duration of the pulse, 200 ps, the total magnetic field is tilted (see <cit.>).  
This is known to be an efficient way of exciting nonlinear skyrmion dynamics <cit.>.
The saddle point searches show that an excitation of the elliptical mode can induce duplication. 
Fig. <ref> shows shapshots from the simulation.

During the pulse, an elongated structure is formed.  After the pulse has been turned off, it splits up and eventually two skyrmions remain in the system.

In addition to the elliptic elongation, the duplication also requires the deformed structure to be bent in order for it to divide up. The bending corresponds to another eigenmode of the elongated skyrmion.


More efficient and/or reliable ways of achieving the duplication process could be devised, for example by pinning the skyrmion, introducing defects or by other more elaborate external stimuli. 
The simulation described here is merely meant to show that it is possible to induce skyrmion duplication in a rather simple way.









We have presented here a method that can be used to search for transition mechanisms and determine the activation energy for transitions in magnetic systems. It represents an adaptation of mode following methods that have been used for several years in studies of atomic rearrangements. The fact that magnetic transitions involve rotations of the magnetic vectors leads to constraints that need to be incorporated into the Hessian matrix in order to determine the low lying eigenvalues and corresponding eigenvectors. 
Saddle point searches where only the initial state of the transition is specified, such as the ones carried out here, are more challenging than a minimum energy path calculation where both the initial and final states are specified, 
but they have the advantage of being able to reveal unexpected transitions and unknown final states.

In the application calculations presented here, we have chosen to drive the system out of the convex region by following a particular
mode.
We note that, eigenvalues may cross and modes change direction 

as the system is pushed along a specific mode, so care should be taken to remain on the same mode throughout the saddle point search
(see <cit.>). 





Applications of the saddle point search method presented here to three-dimensional systems can be expected to yield an even larger variety of mechanisms than for two-dimensional systems since more possibilities open up with the additional dimension. For example, cylindrical skyrmions can contract and form bobbers <cit.> and possibly other previously unknown magnetic structures.

A saddle point search method that only requires the initial state as input 

(unlike the GNEB method where the final state also needs to be specified),
can be used to sample an energy surface in a systematic way and search for possible states of the system, each state corresponding to a local energy minimum.
The search takes the system from one local minimum to another via first order saddle points. This can be used as the basis of a simulation of the long time scale evolution of a system that undergoes thermally activated transitions that may be enhanced by an external field.
For each state visited, multiple saddle point searches are carried out to identify the most relevant, low energy saddle points. The transition rate corresponding to each saddle point can be estimated using harmonic transition state theory, and a random number used to choose which transition will occur next based on the normalized rates. This is referred to as adaptive kinetic Monte Carlo algorithm <cit.>. The simulated time can be estimated from the sum of rates.
This approach has been used as the basis for simulations of long time scale evolution of complex atomic systems where atomic rearrangement events can involve non-intuitive displacements of multiple atoms <cit.>. The adaptation presented here of the saddle point search method to spins opens the possibility of carrying out such long time scale simulations of magnetic systems.





	This work was funded by Icelandic Research Fund (grants 185405-051 & 184949-051),
Academy of Finland (grant 278260), Russian Foundation of Basic Research (grant RFBR 18-02-00267 A),
EU-H2020 MAGicSky (grant 665095), DARPA TEE 
(#HR0011831554) from DOI, and DFG CRC 1238 (project C01). 






99



Kiselev11
N. S. Kiselev, A. N. Bogdanov, R. Schäfer, and U. K. Rössler, 
Chiral skyrmions in thin magnetic films: new objects for magnetic storage technologies? 
J. Phys. D: Appl. Phys. 44, 392001 (2011).


Fert17
A. Fert, N. Reyren and V. Cros,
Magnetic skyrmions: advances in physics and potential applications,
Nature Reviews Materials 2, 17031 (2017).


Fert13
A. Fert, V. Cros, and J. Sampaio,
Skyrmions on the track,
Nature Nanotechnology 8, 152 (2013).


Muller17
J. Müller, 
Magnetic skyrmions on a two-lane racetrack. 
New J. Phys. 19, 025002 (2017).


Muller15
J. Müller, and A. Rosch, 
Capturing of a magnetic skyrmion with a hole,
Phys. Rev. B 91, 054410 (2015).


Hanneken16
C. Hanneken, A. Kubetzka, K. von Bergmann, R. Wiesendanger, 
Pinning and movement of individual nanoscale magnetic skyrmions via defects,
New J. Phys. 18, 055009 (2016).


bessarab_method_2015
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. Commun. 196, 335 (2015).



Lobanov16
I. Lobanov, H. Jónsson, and V. M. Uzdin,
Mechanism and activation energy of magnetic skyrmion annihilation obtained from minimum energy path calculations,
Phys. Rev. B 94, 174418 (2016).



Stosic17
D. Stosic, J. Mulkers, B. Van Waeyenberge, T. Ludermir, and M.V. Milosević, 
Paths to collapse for isolated skyrmions in few-monolayer ferromagnetic films. 
Phys. Rev. B 95, 214418 (2017).


Cortes-Ortuno17
D. Cortés-Ortuno et al.,
Thermal stability and topological protection of skyrmions in nanotracks,
Scientific Reports 7, 4060 (2017).


Uzdin18
V.M. Uzdin, M.N. Potkina, I.S. Lobanov, P.F. Bessarab and H. Jónsson,
Energy surface and lifetime of magnetic skyrmions,
Journal of Magnetism and Magnetic Materials 459,  236 (2018).


bessarab_2018
P. F. Bessarab, G. P. Müller, I. S. Lobanov, F. N. Rybakov, N. S. Kiselev, H. Jónsson, V. M. Uzdin, S. Blügel, L. Bergqvist, and A. Delin, 
Lifetime of Racetrack Skyrmions,
Sci. Rep. 8, 3433 (2018).



Malottki_2018
S. von Malottki, B. Dupé, P. F. Bessarab, A. Delin, and S. Heinze, 
Enhanced skyrmion stability due to exchange frustration. 
Sci. Rep. 7, 12299 (2017).


bessarab_harmonic_2012
P. F. Bessarab, V. M. Uzdin, and H. Jónsson, 
Harmonic Transition-State Theory of Thermal Spin Transitions,
Phys. Rev. B 85, 184409 (2012).



bessarab_potential_2013
P. F. Bessarab, V. M. Uzdin, and H. Jónsson, 
Potential Energy Surfaces and Rates of Spin Transitions.
Zeitschrift Für Physikalische Chemie 227, 1543 (2013).



Romming_2013
N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger,
Writing and deleting single magnetic skyrmions,
Science 341, 636 (2013).


Kubetzka_2017
A. Kubetzka, C. Hanneken, R. Wiesendanger, and K. von Bergmann, 
Impact of the skyrmion spin texture on magnetoresistance,
Phys. Rev. B 95, 104433 (2017).


Bessarab_2017
P. F. Bessarab, 
Comment on ”Path to collapse for an isolated Néel skyrmion”, 
Phys. Rev. B 95, 136401 (2017).


Peters17
B. Peters,
Reaction Rate Theory and Rare Events
(Elsevier Science & Technology, 2017).


Jonsson11
H. Jónsson,
Simulation of Surface Processes, 
Proceedings of the National Academy of Sciences 108, 944 (2011).


Johannesson01
G.H. Jóhannesson and H. Jónsson,
Optimization of Hyperplanar Transition States.
J. Chem. Phys.  115, 9644 (2001).


Bessarab_2014
P. F. Bessarab, V. M. Uzdin, and H. Jónsson,
Calculations of magnetic states and minimum energy paths of transitions using a non-collinear extension of the Alexander-Anderson model and a magnetic force theorem,
Phys. Rev. B 89, 214424 (2014).



SupportingInfo
See Supplemental Material

for details on the methodology, description of the 2-dimensional test problem, 
parameter values used in the calculation described in the main text, 
as well as another set of values and the corresponding minimum energy paths. 


FirstOrderSP
At a first order saddle point the gradient of the energy vanishes and the Hessian, 
the matrix of second derivatives of the energy, has one and only one negative eigenvalue.
 

Pedersen09
A. Pedersen and H. Jónsson,
Simulations of Hydrogen Diffusion at Grain Boundaries in Aluminum,
Acta Materialia 57, 4036 (2009).


Henkelman99
	G. Henkelman, and H. Jónsson.
A Dimer Method for Finding Saddle Points on High Dimensional Potential Surfaces Using Only First Derivatives,
J. Chem. Phys. 111, 7010 (1999).



Gutierrez17
	M. P. Gutiérrez, C. Argáez, and H. Jónsson. 
Improved Minimum Mode Following Method for Finding First Order Saddle Points.
    J. Chem. Theo. Comput. 13, 125 (2017).



    
absil_extrinsic_2013
P.-A. Absil, R. Mahony and J. Trumpf. 
An Extrinsic Look at the Riemannian Hessian. 
In Geometric Science of Information 361-368. 
Lecture Notes in Computer Science.
Springer, Berlin, Heidelberg (2013).


    








spirit
Spirit - a spin simulation framework (see <https://spirit-code.github.io> )

Rybakov_15
Rybakov, F. N., Borisov, A. B., Blügel, S., Kiselev, N. S. 
New type of particlelike state in chiral magnets.
Phys. Rev. Lett. 115 117201 (2015).



lin_internal_2014
S-Z. Lin, C. D. Batista, and A. Saxena, 
Internal Modes of a Skyrmion in the Ferromagnetic State of Chiral Magnets
Phys. Rev. B 89, 24415 (2014).
















Heo16
C. Heo, N. S. Kiselev, A. K. Nandy, S. Blügel, and T. Rasing,  
Switching of chiral magnetic skyrmions by picosecond magnetic field pulses via transient topological states, 
Sci. Rep. 6, 27146 (2016).


Henkelman_2001
G. Henkelman and H. Jónsson, 
Long time scale kinetic Monte Carlo simulations without lattice approximation and predefined event table, 
J. Chem. Phys. 115, 9657 (2001).


Chill14
S. T. Chill, M. Welborn, R. Terrell, L. Zhang, J-C. Berthet, A. Pedersen, H. Jónsson and G. Henkelman,
EON: Software for long time simulations of atomic scale systems,
Materials Science and Engineering 22, 055002 (2014).


Pedersen09
A. Pedersen and H. Jónsson, 
Simulations of Hydrogen Diffusion at Grain Boundaries in Aluminum, 
Acta Materialia 57, 4036 (2009).





















Supplemental Material for “Duplication, collapse and escape of 
magnetic skyrmions revealed using a systematic saddle point search method”









0.3 true cm

Gideon P. Müller,^1,2,3 Pavel F. Bessarab,^1,4 Sergei M. Vlasov,^1,4 Fabian Lux^2,3
Nikolai S. Kiselev,^2 Stefan Blügel,^2,3 Valery M. Uzdin,^4,5 and Hannes Jónsson^1,6




















0.1 true cm
^1Science Institute and Faculty of Physical Sciences, University of Iceland, VR-III, 107 Reykjavík, Iceland
^2Peter Grünberg Institut and Institute for Advanced Simulation, Forschungszentrum Jülich and JARA, 52425 Jülich, Germany
^3RWTH Aachen University, D-52056 Aachen, Germany
^4ITMO University, 197101, St. Petersburg, Russia
^5Dpt. of Physics, St. Petersburg State University, St. Petersburg 198504, Russia
^6Dpt. of Applied Physics, Aalto University, FIN-00076 Espoo, Finland






 
0.5 true cm


section


§ EVALUATION OF THE HESSIAN FOR SPIN SYSTEMS


Eigenvectors of the Hessian matrix are needed in the saddle point search method.
However, 

the second derivatives do not have a direct

geometrical meaning and

need to be treated within the theory of Riemannian manifolds
 <cit.>.

Below, a representation is derived for the Hessian of spin systems where the spin length is fixed,

i.e where the manifold ℳ_phys of physical states is composed of the direct product of N spheres

    ℳ_phys = ⊗_i=1^N S^2 ⊂ℝ^3N.

Therefore, ℳ_phys is a submanifold of the embedding euclidean space ℰ = ℝ^3N.


Due to the singularities that can arise when using spherical coordinates, we choose to use a 3N cartesian representation of the spin coordinates. 


As the Hamiltonian ℋ can be extended to a function ℋ̅ which is defined on ℰ, the derivatives can be calculated easily.
We denote the gradient taken in the embedding space ℰ as ∂ℋ̅ and the gradient taken on the manifold  ℳ_phys as P_x∂ℋ̅, i.e., the projection of the gradient onto the tangent space.
The matrix of second derivatives in the embedding space is denoted ∂^2ℋ̅.

In order to calculate the second derivatives on ℳ_phys while keeping the description in Cartesian coordinates, the projector approach described in Ref. <cit.> is used.
For any scalar function f on the manifold ℳ_phys, this true, covariant Hessian is defined as

    Hess f(x)[z]  = P_x∂^2 f̅(x) z + W_x (z, P_x^⊥∂f̅).



Here, W_x denotes the Weingarten map which for any vector v at a point x of a spherical manifold is given by

    W_x(z,v) = -z x^T v,

where z is a tangent vector to the sphere at x.

In order to calculate the Hessian, v=P_x^⊥∂ℋ̅ is inserted to give

    W_x (z, P_x^⊥∂ℋ̅)
          = -z x^T P_x^⊥∂ℋ̅
          = -z x^T x x^T ∂ℋ̅
          = -z x^T ∂ℋ̅,

where x^T ∂ℋ̅ is the scalar product of the spin with the gradient. 


In Spirit <cit.>, the Hessian is implemented in matrix representation, so we switch notation and drop the subscript x.
For spin indices i and j, we write the gradient ∂ℋ̅ as ∇_iℋ̅ (now to be understood as a three-dimensional object) and the second derivative ∂^2ℋ̅ as H̅.

We reformulate Eq. (<ref>) within the Euclidean representation as a 3N×3N matrix 

    H|^3N = (H_ij|^3N ) = 
    [ H_11|^3N H_12|^3N        ⋯; H_21|^3N H_22|^3N        ⋯;        ⋮        ⋮        ⋱ ]

which consists of N^2 blocks that correspond to the different spin-spin subspaces. This matrix representation  of the Hessian is obtained by acting with Eq. (<ref>) on the euclidean basis vectors of the extended space. Here, each of the 3× 3 subspace matrices is defined as

    H_ij|^3N = P_i H̅_ij - δ_ij I n⃗_j·∇_jℋ̅,

where I denotes the 3×3 unit matrix and n⃗_j the normalized direction of spin j. 
The resulting matrix H|^3N will, however, have N eigenvectors orthogonal to the tangent space of ℳ_phys, representing unphysical degrees of freedom in the embedding space ℰ.

In order to remove these unphysical degrees of freedom, 

the matrix is transformed using the tangent basis to ℳ_phys and written as
H_ij = T_i^T H_ij|^3N T_j, where T_i is the basis transformation matrix of spin i fulfilling T^TP = T^T and T^TT = I|^2N.

Thereby, the true Hessian of Eq. (<ref>), in the 2N×2N matrix representation, becomes H = (H_ij) with the spin-spin blocks defined as

    H_ij = T_i^T H̅_ij T_j - T_i^T I (n⃗_j·∇_jℋ̅) T_j,


and which now contains only the physical degrees of freedom.
From this matrix, the eigenmodes λ|^2N are calculated. The 3N representation is obtained as λ|^3N = Tλ|^2N.




The transformation matrix T_i can be a 3×2 matrix of two tangent vectors to spin i, that can be obtained from any random vector and another vector found by orthogonalization.
Here, we have chosen the unit vectors of spherical coordinates θ and φ, which gives



    T = {e⃗_θ, e⃗_φ}
          = [ cosθcosφ    -sinφ; cosθsinφ     cosφ;    -sinθ        0 ]
    
          = [ zx/r_xy -y/r_xy; zy/r_xy  x/r_xy;  - r_xy       0 ]

where r_xy = sinθ = √(1-z^2).

Note that the poles need to be excluded, but one may simply choose e⃗_x and e⃗_y and orthogonalize them with respect to the spin vector.






Finally, the Hessian matrix in the embedding space ℰ = R^3N is needed, denoted H_ij|^3N. 





For a Hamiltonian of the form


    2ℋ =
          - μ_S ∑_i=1^NH⃗·n⃗_i
         - ∑_i=1^N K_i (K̂⃗̂·n⃗_i)^2
         - ∑_⟨ij|⟩  J_ij n⃗_i·n⃗_j
          -  ∑_⟨ij|⟩ D⃗_ij· (n⃗_i×n⃗_j),

where H⃗ is the external magnetic field, K are uniaxial anisotropies 
(set to zero in the calculations presented here), J_ij are the exchange constants
and D⃗_ij the Dzyaloshinskii-Moriya vectors of the unique interaction pairs denoted by ⟨ij|$⟩,
the spin-spin matrix blocks ofH̅=∂^2ℋ̅become

 0.2 true cm

diagonal blocks:

    H̅_i=j = -2K_i
        [ K_xK_x K_xK_y K_xK_z; K_yK_x K_yK_y K_yK_z; K_zK_x K_zK_y K_zK_z ]
off diagonal blocks:

    H̅_i≠ j =
        [      -(J_ij+J_ji)  (-D⃗_ij+D⃗_ji)_z -(-D⃗_ij+D⃗_ji)_y; -(-D⃗_ij+D⃗_ji)_z        -J_ij-J_ji  (-D⃗_ij+D⃗_ji)_x;  (-D⃗_ij+D⃗_ji)_y -(-D⃗_ij+D⃗_ji)_x        -J_ij-J_ji ].
 0.5 true cm

section

§ SADDLE POINT SEARCH METHOD




The saddle point search method is an iterative process where the spin configuration is modified according to
a force acting on the spins.
Each iteration involves finding the Hessian and its eigenvalues and eigenvectors.
An important aspect of the adaptation of this method to spin systems is the evaluation of the Hessian for 
the curved configuration space.


The optimization method involves rotating the spins in the direction of the force using the velocity projection method described in <cit.> until the force is zero.
Outside the convex region, where one or more of the eigenvalues of the Hessian is negative, 
the effective force,F⃗^eff, is given by Eq. (4) in the main text.





There are various approaches for escaping from the convex region around the initial state minimum:



  *  Follow a (constant) random direction.

  *  Make an initial (small) random displacement and then follow the gradient.

  *  Pick one of the eigenvectors of the Hessian and follow it. 




In the illustrative test problem shown in Fig. 1 in the main text, the second choice was made (due to
the reduced dimensionality of the configuration space). In the saddle point searches for the skyrmion the third choice was made.

The iterative optimization usingF⃗^effmay lead back into the convex region.

Also, eigenvalues may become degenerate and the lowest eigenvalue mode may change significantly throughout a calculation. When following a given mode, the scalar product of the current and previous mode is evaluated and, if a significant change in direction has occurred, the search is carried out for the mode that has the largest scalar product with the previous mode. This approach works due to the fact that eigenmodes are mutually orthogonal and should change continuously throughout 
configuration space.

In special cases the gradient may become orthogonal to the mode to within a given tolerance.

This case needs to be treated specially.

The following algorithm was used in the present studies:



  *  When the mode is not orthogonal to the gradient
    
  *  if the mode has a negative eigenvalue: follow F⃗^eff
  *  if the mode has a zero or positive eigenvalue: 

                  follow the mode.
    
  *  When the mode is orthogonal to the gradient
    
  *  if the eigenvalue is zero (to within a tolerance): follow the gradient
    
  *  if the eigenvalue is positive: follow F⃗^eff.
     0.5 true cm


section

§ PARAMETER VALUES
 1.  For single spin illustration 0.3 true cm

The energy surface of the single-spin system shown in Fig. 1 in the main text to illustrate the saddle point search method 
is a sum of 7 Gaussians of the form

    H = ∑_i H_i  = ∑_i a_i exp(  -l_i^2(n⃗)/2σ_i^2),

with parameters given in Table 1.


 0.5 true cm




 2.  For Hamiltonian in the skyrmion calculations 0.3 true cm

The values of the parameters in the Heisenberg type Hamiltonian, Eq. (4) in the main text, 
were chosen to be the same as in a previous theoretical study <cit.>,D=0.45 JandH= 0.8 H_D,







whereH_D = D^2/(μ_S J)is the critical field strength. For the parameter values used here
and setting J=1 meV, givesH_D= 3.5 T.


 0.5 true cm



§ MINIMUM ENERGY PATHS


In addition to the minimum energy paths shown in Fig. 2 in the main text, 
a calculation was carried out with another set of parameters,D=0.7 JandH=0.71 H_D.






In this case, the activation enegy is quite similar for the three transitions, the duplication barrier being lowest, as can be seen from Fig. <ref> in this Supplemental Material.
Here, the skyrmion state is lower in energy than the ferromagnetic state. 


 0.5 true cm



§ DYNAMICAL SIMULATION


In the dynamical simulation of the duplication process, a magnetic pulse is applied over a time period of 200 ps.
The pulse has a magnetic field directed along the vector {x,z}={0.8,-0.61} while
the static magnetic field is directed along the normal along the 
vector {x,z}={0,1} so the pulse partly cancels out the static field and leads to a tilt.
The static field has magnitudeH=0.8 H_Dbut during the pulse the total magnetic field drops toH_p=0.63 H_D.D=0.45 Jas in the 
calculations shown in Fig. 2 in the main text. The time evolution was carried out using the semi-implicit algorithm of Mentink
 et al.<cit.> with a time step of 0.01 ps.





99Nakahara2003
M. Nakahra.
Geometry, topology and physics.
CRC Press (2003).
absil_extrinsic_2013Supp
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).

spiritSuppSpirit - spin simulation framework (see <https://spirit-code.github.io>).

bessarab_method_2015Supp
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. Commun.  196, 335 (2015).

Rybakov_15Supp
Rybakov, F. N., Borisov, A. B., Blügel, S., Kiselev, N. S. 
New type of particlelike state in chiral magnets.
Phys. Rev. Lett.115 117201 (2015).

Mentink10
J.H. Mentink, M.V. Tretyakov, A. Fasolino, M.I. Katsnelson and Th. Rasing,
Stable and fast semi-implicit integration of the stochastic Landau-Lifshitz equation,
J. Phys.: Condens. Matter  22 176001 (2010).