Thermal stability and topological protection of skyrmions in nanotracks

Magnetic skyrmions are hailed as a potential technology for data storage and other data processing devices. However, their stability against thermal fluctuations is an open question that must be answered before skyrmion-based devices can be designed. In this work, we study paths in the energy landscape via which the transition between the skyrmion and the uniform state can occur in interfacial Dzyaloshinskii-Moriya finite-sized systems. We find three mechanisms the system can take in the process of skyrmion nucleation or destruction and identify that the transition facilitated by the boundary has a significantly lower energy barrier than the other energy paths. This clearly demonstrates the lack of the skyrmion topological protection in finite-sized magnetic systems. Overall, the energy barriers of the system under investigation are too small for storage applications at room temperature, but research into device materials, geometry and design may be able to address this.

show a continuous line that connects the discrete data points m z from every lattice site across the nanotrack long side (in the x direction). Small circles indicate points where m z = 0 and using the distance between them we obtain the skyrmion size.

NANOTRACKS
In the main text we specified a range of DMI values to analyse the minimum energy transitions of the skyrmions. Specifically, we analysed skyrmions for DMI values from D = 0.586 meV to 0.811 meV in steps of 0.045 meV. We obtained these skyrmions after relaxing them using the Landau-Lifshitz-Gilbert equation. In Supplementary Fig. S2 we show the skyrmion profile for every case, by plotting the out of plane magnetisation component m z at the middle and across the long edge of the tracks. Using the convention of defining the skyrmion size as the width within the points where m z = 0, we can see that the skyrmion size increases with larger DMI magnitudes. In addition, there is a small tilting at the boundaries which makes a skyrmion on the system to have a magnitude for the total topological charge, slightly smaller than one (or not perfectly zero for the ferromagnetic state).

Initial states
Given an initial energy band, the Nudged Elastic Band Method (NEBM) will determine the lowest energy transition path between two equilibrium configurations that is accessible from the initial band. In general, it is difficult to predict an optimal initial guess for the bands, in particular for a system with such a high degree of complexity. In the context of the skyrmion transition, linear interpolations seem a robust starting point and are easily reproducible and applicable to different magnetic systems. A similar technique for initialising the system, where spin directions are interpolated, is using Rodrigues formula S1 .
After exploring the system numerically, we also noticed the major role of the boundaries in the skyrmion annihilation and it was necessary to create an initial path close to this transition, which yields lower energy transition paths. Another way of guessing an initial path would be, for example, to use the dynamics that the system follows when applying spin polarised currents (as Sampaio et al. do in Ref. S3) or magnetic fields (that is not the best choice since the Zeeman energy changes the energy landscape configuration). For the systems of this study, these paths converge to the same results obtained with linear interpolations .

Energy band images
In our study we found three main energy paths for the skyrmion annihilation. For clarity, we provide here more detailed versions of the bands for every transition. Firstly, in Supple- These images can change since they are allowed to move during the iterations. (a) Evolution using the NEBM and starting with linear interpolations for the initial bands. (b) Evolution using the climbing image NEBM, applied to the bands of (a) after relaxation (last NEBM step).

Skyrmion destruction evolution
According to our implementation of the NEBM, we checked how the maximum energy of the band evolves according to the iterations of the algorithm. We compute this magnitude using the image with the largest energy in the band, which can change since images move as the NEBM evolves. In Supplementary Fig. S5 we show this evolution for the case of the skyrmion destruction discussed in the main text. In particular, Supplementary Fig. S5a refers to the skyrmion destruction where the bands for D = 0.721 meV and above relax to the singularity mediated skyrmion destruction, and for magnitudes of D = 0.676 meV and below, the bands relax to the skyrmion collapse. We can observe that for the smaller DMI values, the maximum energy value oscillates before reaching convergence and this is due to the movement of images around the saddle point (where the last spins defining the skyrmion center flip) whose energy fluctuate dramatically when crossing this region.
Moreover, Supplementary Fig. S5b depicts the evolution of the climbing image NEBM technique applied to the relaxed bands of Supplementary Fig. S5a and all the bands relaxed towards the skyrmion collapse. Bands with D = 0.721 meV and above take a long time to relax, specially when the DMI gets stronger. As we specified in the main text, we believe this is due to the rough shape of the energy landscape around the saddle point.

ENERGY BARRIERS FOR WEAK DMI
We have analysed the energy barriers for DMI magnitudes in the range 0.586 to 0.811 meV. From Supplementary Fig. 5 in the manuscript, the tendency is that the skyrmion collapse curve decreases almost linearly in this range and extrapolating the curve it seems that at some small D value it will cross the boundary annihilation curve. We then computed the values of the energy barriers for DMI magnitudes weaker than 0.586 meV. Decreasing in steps of 0.045 meV, we could not stabilise a skyrmion for D smaller than 0.405 meV, probably due to the strong exchange and anisotropy which are kept fixed. For this small DMI magnitude a skyrmion is only 3.09 nm wide.
We show in Supplementary Fig. S6 the energy barriers computed for the whole range of DMI values in the main study and down to 0.405 meV. In this plot, we show both the energy barriers obtained using the image with the largest image in a band, which we call NEBM curve, and the energy barriers obtained after performing a polynomial interpolation on the bands, named Interpolation curve. In general, the interpolation gives optimal results for the boundary annihilation, which can be seen from Supplementary Fig. 4  point of the band in most cases, but sometimes it can slightly overestimate the energy barrier when there is no good resolution around the saddle point, which is given by the reversal of the spins at the skyrmion nucleus. In particular, this happens for D = 0.541 meV, which is shown in Supplementary Fig. S6b. Nevertheless, the quadratic tendency is clear when using the points from the NEBM curve. Interestingly, when comparing the inset of Supplementary   Fig. S6b with Supplementary Fig. S6a at D = 0.450 meV, where the skyrmion is only 4.16 nm wide, the energy barrier of the skyrmion collapse is slightly smaller. Decreasing the DMI even further we could not obtain a clear value for the boundary annihilation, however in this range the skyrmion collapse barrier is smaller than 0.002 meV.

TOPOLOGICAL CHARGE
The topological charge gives us an idea about the contribution of the spin directions to the mapping of the magnetisation field into a unit sphere. This helps us to characterise non trivial magnetic structures, such as the Bloch point like structure we observed in our simulations. From our results, a full mapping of the spins around this topological singularity would not be optimal, since the full structure usually lies between two images around a saddle point. Nevertheless, it is still possible to observe large variations of the topological charge density when the singularity starts to appear, allowing us to distinguish it in the spin field.
In micromagnetics, instead of calculating a topological charge, it is standard to compute the so called skyrmion number, which requires the calculation of the derivatives of the magnetization field in a two dimensional plane. For a discrete spins lattice we use the topological charge definition rather than discretising the derivatives along the 2D-lattice directions (since, for the sake of numerical calculations, we have to anyways discretise the The spherical angle at the µ-th site Ω µ,1,2 , is defined as with respect to the neighbouring spins s 1 , s 2 as with ρ = [2 (1 + s µ · s 1 ) (1 + s 1 · s 2 ) (1 + s 2 · s µ )] Similarly for the µ-th site Ω µ,3,4 . For a square and hexagonal arrangement, the triangles specified per lattice site are depicted in Supplementary Fig. S7 (it would be equivalent to use the opposite triangles that also cover a unit cell). It is worth noticing that Berg and Lüscher S5 originally introduced a sign for the exponential in equation 3  where s µ · (s 1 · s 2 ) = 0 and 1 + s µ · s 1 + s 1 · s 2 + s 2 · s µ . For these cases, the exponential falls in a branch cut in the complex plane S5,S6 , which implies a change in the topological charge that can be related to the condition where a skyrmion is created or collapses.