Topological defect-mediated skyrmion annihilation in three dimensions

The creation and annihilation of magnetic skyrmions are mediated by three-dimensional topological defects known as Bloch points. Investigation of such dynamical processes is important both for understanding the emergence of exotic topological spin textures, and for future engineering of skyrmions in technological applications. However, while the annihilation of skyrmions has been extensively investigated in two dimensions, in three dimensions the phase transitions are considerably more complex. We report field-dependent experimental measurements of metastable skyrmion lifetimes in an archetypal chiral magnet, revealing two distinct regimes. Comparison to supporting three-dimensional geodesic nudged elastic band simulations indicates that these correspond to skyrmion annihilation into either the helical and conical states, each exhibiting a different transition mechanism. The results highlight that the lowest energy magnetic configuration of the system plays a crucial role when considering the emergence and stability of topological spin structures via defect-mediated dynamics. Skyrmions are topologically non-trivial, vortex-like magnetic structures the dynamics of which have been mostly studied in 2D systems, but they are also able to exist as 3D tube-like structures. Here, the authors report a combination of experimental and computational results investigating the annihilation dynamics of 3D skyrmion structures in order to better understand how to stabilise topological structures in other bulk magnetic systems.

T opological defects have been shown to mediate phase transitions by breaking local symmetries. For example, Kosterlitz and Thouless suggested that the proliferation and propagation of dislocations and disclinations, respectively breaking translational and rotational symmetry, are responsible for the melting of a 2D crystalline solid into a liquid 1 . Such discontinuities are excited states, and can be created either in opposing pairs within the system, or singly at the system's boundaries 2 . Topological defects have been observed in a range of physical systems, such as in the nematic phase of liquid crystals 3 , as fluxons in type-II superconductors 4 , or as magnetic vortices in the two-dimensional X-Y spin model 1 .
Topological phases have been the focus of recent research due to their unique behaviours and properties. In particular, magnetic skyrmions, stabilised by the Dzyaloshinskii-Moriya interaction (DMI), have been investigated due to their potential applications in future spintronic devices [5][6][7] . Similar to the melting of solids by dislocations and disclinations, magnetic skyrmions are annihilated by the proliferation and subsequent propagation of magnetic Bloch points 8 . Such three-dimensional topological defects have been called hedgehog defects, or emergent magnetic monopoles 9 , due to their singular nature. Finding methods of controlling this behaviour will be crucial for realising reading and writing protocols in proposed skyrmion devices 10-12 . Due to the focus on thin-film applications, previous investigations into the stability and annihilation mechanisms of skyrmions have largely focused on isolated skyrmions in twodimensional systems [13][14][15] . In these confined systems, the energy barrier required to annihilate a skyrmion at the boundary of the system is typically lower than that required for its destruction by direct collapse 16,17 . Further calculations have demonstrated that entropic considerations, which concern the number of available paths across the annihilation energy barrier, play a more prominent role in determining the lifetime of skyrmions in twodimensional systems than the energy barrier itself [18][19][20][21] . This phenomenon has been demonstrated experimentally 18,22 .
However, while they are commonly portrayed as twodimensional objects forming a hexagonal skyrmion lattice (SkL), in three dimensions magnetic skyrmions exist as extended tube-like objects 23,24 . In a bulk skyrmion system, such a skyrmion tube (SkT) can exist within either a helical or conical state, as shown by the visualisations in Fig. 1a, b. These two chiral magnetic states consist of spin spirals aligned along either the easy axes of the system (helical), or along the applied field direction (conical) 7,25 . In each configuration, the SkT is thought to be annihilated by a different Bloch-point-mediated mechanism, as depicted in Fig. 1c. On the left, the SkT connects with the edge of the nearby helical domain, forming a pair of anti-hedgehog-like Bloch points (H+ and H−), whose subsequent motion joins the SkT to the neighbouring helical domain 8 . On the other hand, on the right side, the SkT breaks in two, forming a pair of spiral-like Bloch points (C+ and C−), which unwind the SkT into the conical state 26 . Cross-sections of the system are shown in Fig. 1d, including insets showing the spin configuration surrounding each Bloch point. The chirality of these Bloch points is defined by the local topology of the spin structures (see Supplementary Fig. S1, Supplementary Note 1).
While there have been previous computational studies into the topological phase transitions governing three-dimensional skyrmion annihilation [27][28][29][30] , in this work, we undertake a combined experimental and computational investigation into the annihilation dynamics in bulk skyrmion systems. Beyond being vital for technological implementations, such Bloch-point-mediated mechanisms likely hold the key to understanding the stabilisation of more exotic topological phases, for example: the emergence of the low-temperature skyrmion phase from the tilted conical state in Cu 2 OSeO 3 31-33 , the transformation of hexagonal skyrmion lattices to meron-antimeron lattices in Co 8 Zn 9 Mn 3 34 , or future experimental realisation of magnetic hopfion states 35 .

Results
Metastable skyrmion phase diagram. In bulk chiral crystals, skyrmions typically exist in a small range of temperature and applied magnetic field close to the Curie temperature, T C , of the material 25 . Previous studies have shown that a metastable skyrmion state can be realised outside this equilibrium region by cooling the sample under an applied magnetic field [36][37][38] . Being metastable, this state has a finite, temperaturedependent lifetime 39 , which can be studied with time-resolved measurements 40 . It has been noted that disorder and defects in the underlying crystal lattice increase the lifetime of the metastable skyrmion state 36 . Therefore, in order to stabilise metastable skyrmions over the widest possible range of temperature and applied magnetic fields, we chose to investigate a Zn-substituted Cu 2 OSeO 3 sample, which exhibits such an enhanced lifetime 41 . Figure 2a shows the magnetic phase diagram of the chosen (Cu 1−x Zn x ) 2 OSeO 3 (x = 0.02) single crystal, with a T C = 57.5 K (see Supplementary Figs. S2-3, Supplementary Note 2). Being multiferroic, the phase diagram was determined by measurements of the electric polarisation, P, along the [001] crystal axis, as a function of the magnetic field applied along the [110] axis 42 . We chose this sample configuration due to the relative ease of modelling the magnetic phase transitions for these domain orientations (see Supplementary Note 3). At low fields, the ground state consists of helical domains, with their spin modulations oriented along the 〈001〉 axes due to the cubic anisotropy 43 . With increasing applied field, the conical state is realised, with a propagation vector aligned with the magnetic field direction. The phase diagram illustrates the stabilisation of the equilibrium SkL state over a short range of temperature and applied field around~56 K, and the comparably larger extent of the metastable SkL state, realised by field cooling at 20 mT.
Corresponding AC susceptibility measurements are displayed in Fig. 2b, c, performed after zero field cooling (ZFC), high field cooling (HFC) and field cooling (FC) the sample to the target temperature (see Methods). In Fig. 2b, the decrease in the real component of the AC susceptibility, χ 0 , around 20 mT is characteristic of the formation of the SkL state 44 , as highlighted in yellow. Upon FC at 20 mT to 48 K, a similar characteristic decrease in χ 0 is indicative of the presence of metastable skyrmions, as shown in Fig. 2c.
Examination of the magnetic phase diagram reveals that the metastable SkL state overlays both the conical and the helical ground states, as indicated by the yellow hatched and dotted sections on the phase diagram. Therefore, one might expect that the metastable SkL state will decay into either the conical (SkL → C) or helical (SkL → H) states depending on the field applied to the sample. This affords us the opportunity to study the dynamics of both skyrmion annihilation mechanisms proposed in Fig. 1.
Metastable skyrmion lifetime. Measurements of the AC susceptibility, χ 0 , as a function of time allow the annihilation of the metastable skyrmion state to be measured. Specifically, we can expect that the difference in the value of χ 0 between the FC and ZFC or HFC measurements is proportional to the population of metastable skyrmions present in the sample, as has been observed in previous works 41,45 . Thus, when the metastable skyrmions annihilate, the value of χ 0 will relax over time, as indicated by the green vertical arrow in Fig. 2c. We carried out such time-resolved skyrmion lifetime measurements as a function of the applied magnetic field. A selection of the resulting data is displayed in Fig. 3a (full dataset shown in Supplementary Fig. S4 and Supplementary Note 4). Here, the AC susceptibility has been normalised, where χ 0 is the initial value of χ 0 at t = 0, and χ f is the value the system is tending to as t → ∞.
The lifetime of the metastable skyrmion state at each field, τ, was extracted by fitting the data to a stretched exponential decay function 46 , A stretching parameter of β ≠ 1 indicates that the system displays a range of lifetimes. Such a distribution of lifetimes may be due to inhomogeneous magnetic fields within irregularly shaped samples 47,48 . Nevertheless, for our measurements, the inclusion of β is necessary to fit the data, and varies between 0.3 and 0.6 for different applied fields.
The resulting metastable skyrmion lifetimes are plotted as a function of magnetic field in Fig. 3b. The lifetime is at a maximum at 24 mT, and decreases at magnetic fields above and below this point. We note that the larger errors bars at both low and high values are due to the difficulty of fitting either very short or long lifetimes-we are effectively exploring the full range of measurable lifetimes for our chosen measurement time. The two distinct exponential trends point towards observation of the two anticipated regimes of the SkL → H and SkL → C annihilation mechanisms, and we tentatively label the corresponding regions of the figure accordingly. In reality, since the decay is a statistical process, one might expect a crossover region in the middle of the field range, with the simultaneous occurrence of both decay pathways.
A modified Arrhenius' law can be employed to describe the temperature dependence of the skyrmion lifetimes, allowing the energy barrier E B governing the decay mechanism to be determined, Here, the temperature dependence of the energy barrier is assumed to be approximately linear, with a as the proportional constant, following the relationship E B /k B = a(T − T s ) 40 . The lifetime can be expected to reach a minimum just below the lowest temperature extent of the equilibrium skyrmion phase, T s , which we determined to be approximately 4 K below T C (see Supplementary Fig. S2). Considering equation (2), the two distinct exponential trends in Fig. 3b suggest both a negative or positive linear dependence of a with the applied magnetic field. However, one must also consider the contribution of the τ 0 term. This prefactor is typically thought of as a characteristic attempt frequency with which the system attempts to overcome the energy barrier. The term also includes the aforementioned entropic correction which has been shown to play a crucial role in the stability of magnetic skyrmions in thin lamellae 22 .
In order to separate the contributions of a and τ 0 , we performed further, temperature-dependent lifetime measurements at five applied magnetic fields between 16 and 34 mT (full data shown in Supplementary Fig. S4). The resulting extracted lifetimes are plotted as a function of temperature in Fig. 3c. By fitting these lifetimes with the modified Arrhenius law in equation (2), the values of a (the gradient) and τ 0 (the value of τ at T s ) were determined, and are plotted in Fig. 3d.
Immediately, one can see that a appears to decrease either side of 24 mT, as was indicated by the initial lifetime measurements in Fig. 3b. This supports the interpretation that the energy barrier appears to vary linearly, with either a positive or negative gradient, in two regimes which correspond to the SkL → H and SkL → C transitions respectively. However, the contribution from the τ 0 term appears to decrease exponentially across the measured field range. This relationship is significantly different to the one determined by Wild et al. in their thin lamella sample 22 . Following the phenomenological Meyer-Neldel compensation rule 21 , they demonstrated that τ 0 varied exponentially with the measured energy barrier-concluding that the entropy therefore varies linearly with the energy barrier. However, when plotting our measured a and τ 0 values in the inset of Fig. 3d, it is clear that such a dependence is not recovered for our bulk sample, and entropic considerations therefore appear to play a less prominent role.
The discrepancy can be explained by considering that in bulk measurements, the annihilation of skyrmions consists of two processes: the nucleation of Bloch-point pairs, and the subsequent motion of these Bloch points along the length of the skyrmion tube. In bulk systems, τ 0 appears to be correlated with increased disorder and defects within the underlying crystal lattice 45 . Therefore, the decrease of τ 0 with increasing magnetic field in our results indicates that the thermally-activated depinning of the Bloch-point from structural defects becomes more energetically favourable at higher applied magnetic fields. Any entropic contributions may well be obscured by these pinning effects.
Energy barrier simulations. To support the interpretation of the experimental measurements, we modelled a chiral magnet by means of a classical spin model with periodic boundaries in the xy plane and open boundaries in the z direction (see Methods). To begin, suitable magnetic configurations for the helical and conical states were prepared and their energy was numerically minimised over the field range. This was repeated for single SkTs embedded within helical and conical states, resulting in the configurations depicted in Fig. 1a To investigate the transitions of isolated SkTs into helical and conical states, the geodesic nudged elastic band method (GNEBM) was employed 13,49 . The GNEBM finds minimum energy paths (MEPs) through configuration space, between two equilibrium states by minimising the energy of a chain of copies of the system known as images. This allows the determination of the energy barrier for the transition. (see Methods, Supplementary Notes 10 and 11). Comparing the energies of magnetic states is not sufficient to describe their stability against one another-it is only by considering the microscopic decay pathway, and its associated energy barrier, that one can gain insight into skyrmion stability.
We first applied this algorithm to determine a MEP for a SkT annihilating within a helical state (SkL → H) for applied magnetic fields between b z = 0.04 and 0.16 (in units of exchange, see "Methods"), as shown in Fig. 4a-i (full results shown in Supplementary Fig. S7, Supplementary Movie 1). Since previous measurements suggest that the skyrmion state prefers to annihilate into helical domains with magnetic modulation in the same plane as the skyrmion lattice state 8,50 , we chose to utilise helical domains aligned along the x-axis accordingly. We found it was necessary to apply a cubic anisotropy of k c = 0.05 in order to stabilise the helical domains against the conical state over the investigated field range.
As shown by the selected visualisations of the simulated images in Fig. 4c, the SkT decays by forming a pair of Bloch points, whose subsequent motion connects the SkT and the helical domain. The energy of the system along the determined MEP, parameterised by the reaction coordinate X through configuration space, is plotted in Fig. 4a, where each data point is associated with an image of the system. The trajectories of the Bloch points along the z-axis are plotted in Fig. 4b, with the vertical dashed lines indicating images where Bloch points are created or destroyed. By comparison of these panels, one can see that the two energy barriers along the MEP are each associated with the formation of a pair of Bloch points.
We repeated the GNEBM simulations at a range of applied magnetic fields. At some applied fields, both pairs of Bloch points form and annihilate simultaneously, as shown by the simulation performed at b z = 0.10 in Fig. 4d-f. Furthermore, at b z = 0.14 in Fig. 4g-i, we found a simultaneous formation of four pairs of Bloch points. At b z = 0.14 and 0.16, the final state was found to be higher in energy compared to the SkT in helical state, suggesting that, for the chosen parameters, the SkT within the helical state was the lowest energy configuration over a small range of field. Due to the complexity of the transition, there are likely many similar paths through configuration space, all with energy close to the MEP, and therefore we speculate that the formation of different numbers of Bloch points, and whether they form on the left or right of the skyrmion tube, is likely due to the original path interpolation between the specified initial and final states.
There has been a previous work utilising Ginzburg-Landau analysis to investigate energy barriers for a SkT decaying within a helical state. However, the decay mechanism the authors found was more akin to the breaking of the SkT within a conical state 27 , rather than the joining of the SkT to the helical domains as seen in this work, and other experimental works 8,22 .
We performed similar GNEBM simulations for the SkT in a conical state (SkL → C) at applied fields between b z = 0.36 and 0.50, as shown in Fig. 5a-i (full results shown in Supplementary  Fig. S8, Supplementary Movie 2). With the same cubic anisotropy utilised for the helical simulations, we found that the conical state was only the lowest energy state over a short range of applied field (see Supplementary Fig. S9, Supplementary Movie 3). Therefore, in order to simulate the SkL → C transition over a wide field range, we performed this set of simulations with no cubic anisotropy present. Figure 5a-c shows the results of the SkL → C simulation under an applied field of b z = 0.36. In comparison to the SkL → H  simulations, the energy landscape is considerably more complex, with many more metastable states possible along the MEP. In Fig. 5c, images of the system along the MEP show the breaking of the SkT into three sections by the formation of two pairs of Bloch points. The sections consist of two chiral bobbers at the open boundary surfaces of the system 51 , and a short SkT capped by two Bloch points, known as a toron 52,53 . Following along the MEP, the formation of each Bloch-point pair is once again associated with an energy barrier, as shown in Fig. 5a. In addition, at this lower field there is an energy barrier associated with the collapse of the toron (image 8), and the two chiral bobbers (images 10 and 11), suggesting they are metastable states. We confirmed that these states are indeed local energy minima by energy minimisation. The system behaves similarly at b z = 0.42, as shown in Fig. 5d-f. At b z = 0.48 and above, the applied field is sufficient to transform the conical state into the uniformly magnetised state. From this point, only one pair of Bloch points is formed during the SkT annihilation, as shown in Fig. 5g-i, which resembles previous computational works 28,29 . Furthermore, there is no longer an energy barrier associated with the destruction of the torons and chiral bobbers, suggesting they are not metastable configurations at higher applied fields.
To allow comparison of these simulations with the experimental results, we calculated the energy barrier required for the formation of the Bloch-point pairs in each simulation, as indicated by the arrow in Fig. 4d. In cases where two, or four, pairs of Bloch points were formed simultaneously, we divided the energy by the number of pairs. Energy barriers associated with the formation of Bloch-point pairs for the SkL → H and SkL → C transitions are shown as a function of applied field in Fig. 6a, b respectively. In addition, energy barriers for the formation of Bloch-point pairs for SkT tubes annihilating within the uniformly magnetised state (SkL → UM) with cubic anisotropy of k c = 0.05 are shown in Fig. 6c (see Supplementary Figs. S7-S10). Finally, energy barriers associated with the annihilation of toron and chiral bobber states within the conical and uniformly magnetised simulations are shown in Fig. 6d.

Discussion
Considering Fig. 6a, b, one can see that the energy barrier associated with Bloch-point formation increases as a function of the applied field for the SkT within the helical state, whereas it decreases as a function of field for the SkT within the conical state. This qualitatively agrees with the conclusions drawn from the experimental data in Fig. 3, where we anticipated the energy barriers varying linearly with the applied field, supporting our labelling of the two linear regimes in the lifetime data to the two distinct SkT annihilation mechanisms.
Our results of SkT annihilation to the conical state can be compared to previous results of isolated skyrmions in twodimensional systems, which typically show two transition paths: either the skyrmion is annihilated at the boundary of the system, or it is directly annihilated within the system's interior. For direct annihilation, this may occur either by shrinking and flipping the central spin 16 , or by the newly-coined chimera collapse, where the spin flip occurs at the edge of the skyrmion under an in-plane field 54,55 . The shrinking annihilation mechanism resembles our results in Fig. 5, where the SkT narrows prior to the nucleation of the Bloch-point pair, at the point where the spin in the centre of the skyrmion flips to align with the applied field. We also investigated the possibility of surface Bloch-point nucleation, as has been done in previous works 29,56 (see Supplementary Fig.  S11, Supplementary Note 12, Supplementary Movie 4), but when applying the GNEBM, the algorithm appeared to prefer forming Bloch-point pairs within the interior of the system over the simulated field range.
There are a number of limitations in our model which prevent us from reaching a direct quantitative comparison to the experiment. Firstly, in the SkL → H model, we only consider skyrmions decaying directly into helices aligned to the [001] crystal axis. This is because it is simple to conceive of an initial state for a skyrmion embedded in a helical domain with a winding direction in the plane of the skyrmion lattice 8 , and this also seemed to be the dominant decay pathway in a previous neutron scattering experiment 50 . This was also our motivation for choosing to measure the sample with the field applied along the [110] axis, rather than the [111] axis, where the helices would all be oriented at some angle to the applied field. However, it is possible that in the measured bulk sample there are additional decay pathways for skyrmions decaying into the helical domains oriented along the [100] and [010] crystal axes, which we have not considered.
Second, due to the considerable complexity of the model, we only consider the annihilation of an isolated skyrmion in our system, whereas in the experimental bulk sample we are measuring the decay of a skyrmion lattice. One could conceive that the lifetime of an isolated skyrmion may differ from one surrounded by other skyrmions 57 . Third, the decay dynamics in the experimental sample are, by design, influenced by the pinning introduced by defects and impurities in the Zn-substituted sample. Thus, it would make an interesting follow-up study to simulate the effect of defects on the determined energy barriers, as has been done for 2D systems 17,58,59 .
To summarise, our experimental measurements of fielddependent metastable skyrmion lifetimes indicated two distinct regimes, where skyrmion annihilation is dominated by decay into either the helical or conical states. This conclusion is supported Fig. 6 Field dependence of the annihilation energy barriers. a-c The field dependence of the simulated energy barrier which must be overcome to form one pair of Bloch points is plotted as a function of the applied field for the skyrmion lattice to helical (SkL → H) (blue), skyrmion lattice to conical (SkL → C) (red) and skyrmion lattice to uniformly magnetised (SkL → UM) (yellow) simulations. d The simulated energy barrier for the annihilation of a chiral bobber (ChB, squares) or a toron (triangles) is plotted as a function of the applied magnetic field for the SkL → C (red) and SkL → UM (yellow) decay mechanisms. Each dataset has been fitted with a linear trend. by GNEBM simulations, where the determined energy barriers associated with Bloch-point-mediated skyrmion annihilation to the helical and conical states revealed the same field dependence. In contrast to previous studies of two-dimensional systems, in our bulk sample the effects of defect pinning appear to dominate over entropic contributions. The results highlight that, for technological applications, consideration of the competing magnetic ground state and the dimensionality of the system will be crucial when engineering skyrmion stability for reading and writing processes. Further utilisation of the GNEBM, and similar computational methods, will no doubt prove to be a powerful tool in the understanding of topological defect-mediated phase transitions, particularly when considering the stabilisation of future complex three-dimensional spin textures.

Methods
Sample growth and preparation. A single crystal of (Cu 1−x Zn x ) 2 OSeO 3 was grown from polycrystalline powders using the chemical vapour transport method 60 . The composition of the resulting crystal was determined by inductively coupled plasma mass spectroscopy, and energy-dispersive X-ray spectroscopy, yielding x = 0.02 (2% Zn substitution), as reported in previous work 60 . For electric polarisation measurements, the sample was aligned and cut into a plate shape, with two large parallel 111 h i faces.
AC susceptibility magnetometry. AC susceptibility measurements were performed with a MPMS3 Quantum Design magnetometer, fitted with the AC option. The sample was attached to a quartz rod using GE varnish and mounted in the instrument with the [100] crystal direction parallel to the applied magnetic field. AC susceptibility measurements were performed with an oscillating magnetic field amplitude of 0.1 mT and a frequency of 10 Hz. All cooling procedures were performed at a rate of 40 K/min. For the ZFC, FC and HFC measurements, the sample was cooled from 60 K to the target temperature under an applied field of 0, 20 or 200 mT respectively, and measurements subsequently performed as a function of increasing or decreasing magnetic field. For the lifetime measurements, the sample was initialised following the FC process, and the AC susceptibility was then measured as a function of time under stable conditions for three or more hours.
Electric polarisation measurements. For the evaluation of electric polarisation P, silver paste was painted on a pair of [111] surfaces of the Zn-substituted Cu 2 OSeO 3 crystal. With the field applied along the [110] crystal axis, the pyroelectric current was measured during a constant rate of magnetic field sweep using an electrometer (Keithley 6517B), and integrated over time to deduce P 42 .
Computational modelling. The system was modelled using a cubic lattice of N = 30 × 30 × 30 classical spins with periodic boundaries in the xy plane. The model parameters were chosen in order to fit three helical domains within the specified system size, as has been done in previous works 52 , and therefore does not directly relate to the experimental values. We are aiming for a qualitative description of the universal skyrmion annihilation mechanisms, rather than a direct quantitative comparison. The system follows a Heisenberg-like Hamiltonian where energy constants of the different interactions are made dimensionless by normalising by the exchange constant, Here, s i are unit vectors denoting spin orientations, d = D/J > 0 is the DMI constant (for a cubic helimagnet), k c = K c /J > 0 is the cubic anisotropy constant (easy 〈001〉 axes) and b z = (μ s /J)B z is the z-component of the applied field, specified in units of magnetic moment per exchange constant μ s /J. The angled brackets in the first two sums denote counting pairs of spins only once. The magnitude of d ≈ 0.727 was set such that the periodicity of a spin spiral at zero field and in the absence of any anisotropy is 10 lattice sites 61,62 . The value of the saturation field in the absence of anisotropy is defined as B D = D 2 /(μ s J) = d 2 ≈ 0.53 62 . To stabilise helical configurations, a cubic anisotropy constant of k c = 0.05 was chosen, while for the conical simulations, k c = 0.00 was specified. Equilibrium states were obtained by initialising the system with a suitable configuration and minimising the energy using an over-damped Landau-Lifshitz-Gilbert (LLG) equation (see Supplementary Figs. S5, S6).
The GNEBM finds minimum energy paths (MEPs) through configuration space between two magnetic states 13 . Using a suitable initial state, a series of images along the transition path are calculated. During the algorithm, the images follow gradients of energy and are kept distanced in configuration space until reaching a MEP between the initial and final states. Saddle points along the MEP determine the energy barriers of the transitions. In this study, we use a rotation formula for the spin orientations to initiate the algorithm 13 .
All simulations and GNEBM calculations were obtained using Fidimag 16,49 . Further details about the simulations and the applied convergence criteria to stop the energy minimisation for both the LLG equation and GNEBM are discussed in the Supplementary Information.

Data availability
All experimental data analysed and presented in the main text and supplementary materials can be found in an online repository 63 https://doi.org/10.5281/zenodo.4384569. Further material is available from the corresponding author upon reasonable request.

Code availability
All code utilised in the simulations presented in the main text, supplementary text, and Supplementary Movies 1-4 can be found in an online repository 63 https://doi.org/ 10.5281/zenodo.4384569. Further material is available from the corresponding author upon reasonable request.