Stochastic dynamics and pattern formation of geometrically confined skyrmions

Ensembles of magnetic skyrmions in confined geometries are shown to exhibit thermally driven motion on two different time scales. The intrinsic fluctuating dynamics (t ∼ 1 ps) are governed by short-range symmetric and antisymmetric exchange interactions, whereas the long-time limit (t ≳ 10 ns) is determined by the coaction of skyrmion–skyrmion-repulsion and the system’s geometry. Micromagnetic simulations for realistic island shapes and sizes are performed and analyzed, indicating the special importance of skyrmion dynamics at finite temperatures. We demonstrate how the competition between skyrmion mobility and observation time directly affects the addressability of skyrmionic bits, which is a key challenge on the path of developing skyrmion-based room-temperature applications. The presented quasiparticle Monte Carlo approach offers a computationally efficient description of the diffusive motion of skyrmion ensembles in confined geometries, like racetrack memory setups. The anticipated role of skyrmions as information carriers in spintronic devices has, so far, been hampered by difficulties in controlling their motion. Here, the authors use micromagnetic simulations to investigate the temperature-dependent motion of skyrmions, revealing that their magnetic texture reacts on two different time scales.

M agnetic skyrmions [1][2][3] are quasiparticles that are considered as possible carriers of information for future storage devices. Their specific chirality is determined by the antisymmetric Dzyaloshinskii-Moriya exchange interaction (DMI) 4,5 . The DMI can be induced by a broken inversion symmetry in a crystal itself (e.g., MnSi) 2 or by interfacing a heavy metal layer (e.g., Pt, W, Ir) with a ferromagnetic material (e.g., Fe, Co) 6,7 .
Conceptually the utilization of these topologically non-trivial 8 quasiparticles in so-called racetrack setups is of great interest [9][10][11] . Skyrmions can be manipulated (written and deleted) 12,13 and addressed individually 14-16 on a magnetic stripe, allowing a memory device extension from a pure surface density into the third dimension by pushing the quasiparticle-hole-train back and forth, e.g., by applying electrical currents 17,18 . The low threshold driving current 19 along with the small size and high stability of the skyrmions are key features of this concept.
To connect experimental and theoretical model systems with technological applications, investigating the influence of finite temperatures is of crucial importance. The bits on a racetrackbased memory device need to fulfill two main features, stability against external perturbations and addressability. Both are affected by thermal fluctuations, as shown below.
The stability of skyrmions was examined in several publications over the last years [20][21][22][23][24] . Rózsa et al. 21 investigated theoretically periodic two-dimensional Pd/Fe double-layers on Ir(111) and determined the phase diagram as a function of external field and temperature, which includes field-polarized, skyrmion lattice, spin spiral, fluctuation-disordered, and paramagnetic regions. Skyrmion lifetimes in the fluctuation-disordered regime were calculated. The lifetimes of isolated skyrmions in racetrack geometries for Pd/Fe/Ir(111) and Co/Pt(111) systems were investigated in refs. 22,23 , and different mechanisms were revealed for the collapse of a skyrmion inside the track and at the boundary.
The diffusive motion of skyrmions at finite temperature has also attracted significant research attention lately [25][26][27][28] . Previous studies concentrated on diffusion in infinite or extended geometries, but a clarification of the role of the sample shape still seems to be missing in the case where the system size becomes comparable to that of the skyrmions. Effects of this kind directly impede the addressability, which is indispensable when using skyrmions for storage technology, e.g., in racetrack memory devices. As the number of skyrmions during the diffusive motion remains constant, quasiparticle models have been developed for their description in this limit 25,26,29,30 , which are primarily based on the Thiele equation 31 . The advantage of such a collectivecoordinate description over micromagnetic or atomistic spin dynamics simulations is the significantly lower computation cost.
Owing to the finite temporal resolution, imaging techniques on the atomic length scale like spin-polarized scanning tunneling microscopy (SP-STM) 32 or magnetic force microscopy (MFM) 33 cannot access temporally the diffusive motion of the skyrmions. Instead, the time-averaged skyrmion probability distribution is imaged that may well be different from the zero-temperature configuration or a snapshot of a simulation. Here, we discuss the formation, stability, and addressability of diffusive skyrmion configurations. Using micromagnetic simulations, it is shown that the complex interplay of the repulsive interaction between the skyrmions 29,34 along with the confinement effect of the nanoisland and the thermally induced skyrmion diffusion leads to a pattern formation of the skyrmion probability distribution on the nanosecond time scale. A computationally efficient quasiparticle Monte Carlo method is introduced, which is found to yield comparable results to time-averaged micromagnetic simulations regarding the skyrmion probability distribution. The simple implementation and high speed of such a method may make it advantageous over micromagnetic simulations when the actual number of skyrmions has to be determined based on a timeaveraged experimental image.

Results
Skyrmion stability. We study metastable skyrmions and their motion and mobility at finite temperatures. At first, the phase space has to be explored with respect to the external magnetic field and temperature. A detailed description of the phase diagram for extended systems of ultrathin biatomic layers of Pd/Fe on an Ir(111) substrate hosting magnetic skyrmions can be found, e.g., in refs. 21,35 . Here, simulations are performed for the same system in the micromagnetic framework, by solving the Landau-Lifshitz-Gilbert equation (LLG) 36 , as described in the Methods section.
In Fig. 1 possible equilibrium spin configurations are investigated in a 21-nm-diameter nanodisk of 0.4 nm thickness at zero temperature. For each fixed external magnetic field (B = B z ), 20 different randomly generated initial states are considered and relaxed towards the closest local minimum of the total energy hypersurface. Subsequently, the total skyrmion number N Sk of the relaxed states is calculated and plotted against the external magnetic field in Fig. 1 to gain an overview of the possible metastable states. Below the strip-out instability field 37 , skyrmions as localized cylindrical spin structures are not stable and the configuration consists of spin spiral segments. For these the topological number is not a well-defined quantity, as spin spiral structures can extend over the whole island, and therefore the boundaries have a significant effect. This is reflected by the continuous distribution of the topological charge values found below 0.75 T in the red regime in Fig. 1, whereas discrete skyrmion numbers correspond to discrete steps in the total topological charge. The upper boundary of this region is reasonably close to the strip-out field of B = 0.65 T 37,38 determined for extended systems.
With increasing the external field, individual skyrmions may be stabilized in the system, where they will coexist with spin spiral segments. This corresponds to the darker green area in Fig. 1~B ≈ 1 T, where besides the continuous distribution also discrete steps can be observed in the topological charge.  In the following bright green region, stronger magnetic fields lead to shrunken skyrmions, which in turn enables the magnetic island to host a larger number of quasiparticles. In this regime the topological charge is always close to an integer, with the small deviation caused by the tilting of the spins at the edge of the sample. Reaching field values above B ≈ 3.5 T, the system becomes completely field-polarized and skyrmions are not stabilized anymore starting from a random configuration, even at zero temperature. Previous calculations [37][38][39] indicated that skyrmions on the lattice collapse at around B = 4.5 T in the system.
In order to prevent the appearance of spin spiral states and to be able to consider the skyrmions as well-defined quasiparticles, we choose B = 1.5 T for the following calculations, unless mentioned otherwise. Figure 1 shows that several metastable configurations associated with different topological charges are accessible starting from randomly generated initial states. However, in this approach not all stable structures are covered, as the included ones are obtained from relaxing random initial configurations. Different configurations may also be generated in a controlled way by adding skyrmions to the field-polarized state one-by-one, and relaxing the state at zero temperature. This is shown in Fig. 2 for the disk-shaped (Fig. 2a) and a triangular ( Fig. 2b) geometry for B = 1.5 T.
On top of the transitions owing to a variation of the magnetic field, the impact of finite temperature is also of crucial importance. In Fig. 3a temperature effects on the topological charge are displayed. Starting from a relaxed configuration at T = 1 mK, the temperature is increased every 500 ps by ΔT = 2 K and the topological charge is averaged over time at each temperature. The results show a first discontinuity at T 1 ≈ 40 K. Up to this temperature the thermal fluctuations are relatively weak, and the number of well-defined skyrmions inside the sample remains constant. The small standard deviation of the topological number in this regime can be attributed to the thermal motion of the spins at the edge. However, the combination of thermal fluctuations and the repulsion between the skyrmions triggers the escape of a single skyrmion out of the sample at T 1 ≈ 40 K. Hence, the total topological number is changed from N Sk ≈ −5 to N Sk ≈ −4, see Supplementary Movie 1. This change is also shown in the plot of the topological number susceptibility in Fig. 3b, defined as Only a minor deviation from the otherwise smooth function is visible. This means that the thermal fluctuations are still not strong enough to perturb the quasiparticles drastically.
A first characteristic change in the slope of the topological susceptibility can be seen at T 2 ≈ 50 K. Comparing Fig. 3a, b it is obvious that the fluctuations of topological number are increasing with temperature, and around this point the lifetime of skyrmions is reduced below the averaging window. The average topological charge continuously approaches zero above this temperature. In this disordered regime, the lifetime of skyrmions is shorter than the observation time, as the fluctuations allow a collapse inside the sample. Hence, this temperature range is not favorable for information storage applications.
Finally, at about T 3 ≈ 100 K the paramagnetic regime is entered, which is characterized by a completely disordered time-averaged magnetic configuration. Here the average topological charge is very close to zero. This behavior is also indicated by a change of sign in the first derivative dχ N ðTÞ=dTj T¼T 3 .
In the following, the external parameters are fixed to T = T 0 = 15 K and B = B 0 = 1.5 T, shown by the solid red lines in Figs. 1 and 3. This ensures that the thermal fluctuations do not lead to a collapse or escape for the skyrmions, only to a diffusive motion, and that the number of particle-like skyrmions will be constant during the simulation time.
Skyrmion diffusion. With these first results, giving information on the segment of the parameter space we deal with, we will focus on the influence of thermal fluctuations on the dynamics of the skyrmions. In Figs. 4a and 5, results for different initial numbers of skyrmions presented in Fig. 2 are shown for the circular and triangular geometries, respectively. The first row shows the magnetic configuration after 20 ps of simulation time and the second row the final state, which means the magnetization state after our total simulation time of 20 ns. Only the out-of-plane z component is shown in grayscale. The magnetic configuration after only 20 ps is qualitatively very similar to the final one, displaying slightly deformed skyrmion configurations owing to the The saturation magnetization is locally set to zero in the dark gray areas. Scale bar corresponds to 10 nm thermal fluctuations compared to the zero-temperature initial states in Fig. 2. This short timescale, on which the short-range Heisenberg and DMI dominate the dynamics and cause shape deformations, is examined more extensively below. The third row corresponds to the time-averaged z component, calculated over the simulation time divided into 1000 snapshots. In our understanding this result comes as close as possible to realspace scanning-probe measurements, due to the limited time resolution and finite measurement duration in such techniques. According to ref. 37 , typical limits of the time resolution in SP-STM are t res > 5 ms. With the present simulation parameters, we found that increasing the simulation time further does not lead to a significant change in the time-averaged images, indicating that similar observations may be expected on the much longer experimental time scales as well. The observed two characteristic time scales can be attributed to different interaction types of clearly distinguishable energy scales: local deformations on the picosecond time scale of the magnetic texture are dominated by thermal excitations competing with the short-range symmetric and antisymmetric exchange interactions, whereas the translational motion leading to the complex pattern formation over tens of nanoseconds is caused by the repulsive interactions between pairs of skyrmions and between skyrmions and the boundaries.
In case of the disk, it is extremely challenging to make a clear statement about the number of skyrmions in the system based on the time-averaged images only. As the symmetry of the system is a cylindrical one, also the resulting picture is cylindrically symmetric, consisting of concentric bright and dark circles and rings. The small contrast differences along the angular direction within a single diffusive ring-like area (third row in Fig. 4a) arise because of the finite simulation time, but also as an artifact of the finite grid, leading to a deviation from the ideal radial symmetry. For all initial configurations the total skyrmion number is conserved. As the skyrmions repulse each other, the radius of the resulting gray ring increases for higher skyrmion densities, before one of them becomes localized quite strongly in the center. This behavior is analyzed in Fig. 4b, where the out-of-plane magnetization component depending on the distance from the center of the disk is calculated by integration over the angular coordinate. For comparison, the static profile of a single skyrmion in the center of the disk is plotted as well (black curve).  Fig. 5, the lower symmetry of the geometry is reflected in the resulting time-averaged pictures. The symmetry of the spin configuration coincides with that of the sample for one, three, and six skyrmions without thermal fluctuations, as shown in Fig. 2b. In the time-averaged images in the third row of Fig. 5, this results in well-localized skyrmions with a strong dark contrast for these numbers of quasiparticles in the system. Interestingly, also the adjacent configurations show almost the same time-averaged pattern as the highly symmetric configuration, e.g., two or four skyrmions compared with the case of three skyrmions. For the case where one skyrmion is missing from the highly ordered configuration, the remaining skyrmions imitate the absent skyrmion and effectively show up as an additional "phantom skyrmion". In contrast, a surplus skyrmion smears out in the time-averaged picture, so that mostly the highsymmetry points remain observable, even though they can also be slightly broadened-see the 6th and 7th columns in Fig. 5. For an even larger numbers of skyrmions, the repulsive interaction between them in combination with temperature fluctuations is strong enough for the skyrmions to escape at the boundary during the simulation time, as shown in the decreased number of quasiparticles in the second row compared to the first row in the 8th and 9th columns in Fig. 5.
To identify the different time regimes more clearly, the convergence behavior of the time-averaged pictures is investigated. The time-integrated pictures are calculated up to each snapshot time and compared with the averaged image after 20 ns via a matching parameter, where a complete agreement is denoted by a value of 1. The mathematical definition of the matching parameter is given in the Methods. The results are plotted in Fig. 6. Panels a and b show the convergence behavior for different skyrmion occupations for the disk and for the triangular geometry. For the circle, there are no major differences in terms of the convergence speed for different skyrmion numbers. Only the ensembles of three and seven skyrmions exhibit a slightly faster convergence, which is explicable by the resemblance of those states to the close-packed arrangement. As expected, the close-packed array possesses a lower mobility. In the case of the triangle, for three and six skyrmions the value 1 is reached much faster than for the other configurations. This effect can again be explained by the strong localization of the skyrmions in these cases, preventing an exchange of their positions. A different trend is visible for the eight-skyrmion case (purple line in Fig. 6b), which shows a much slower convergence. This delayed behavior arises due to the non-conserved skyrmion number, which can be seen in Fig. 5. Consequently, after the escape of the skyrmions it takes more time until the averaged picture has adapted to this change.
For exploring the stochastic dynamics on the short time scale, which we suppose to be governed by the skyrmion deformations in contrast to the slow skyrmion displacement, multiple simulations are performed for the same initial configuration, two skyrmions inside the disk relaxed at zero temperature, but with different pseudorandom number seeds. The simulations cover a time of 100 ps each with 50 different initializations. After every 0.5 ps, a snapshot is taken and the positions of the two skyrmions extracted. On the scale of the simulation time no noteworthy displacement of each of the skyrmions takes place, however, a deformation of the magnetic textures is visible. The mean skyrmion radius relative to the zero temperature radius calculated as described in the Methods, is averaged over the different simulations, shown in Fig. 7 for each of the skyrmions (Fig. 7a, b). Here, two features are remarkable. First, the radius increases, approaching an enhancement of~10% after 100 ps. Second, the radius is oscillating with a high frequency of approximately 60 GHz, estimated from the average time between the peaks. This oscillation can be identified with internal skyrmion modes, like the so-called breathing modes (see ref. 38 : at T = 0 K, B = 1.5 T and the triangular atomic lattice the  All of the discussed examples show the strong influence of the interaction between skyrmions as well as the role of the geometrical aspects on their mobility, that can vary between delocalization and a strong localization, e.g., in the center of the disk in Fig. 4a. This feature of variable mobility, depending on the geometry and packing density of skyrmions, is of general nature. In real experimental systems, further contributions may affect the mobility, and ultimately the measured position of the skyrmion as well. As an example of such a feature, defects in the grown thin film nanoislands will be discussed. In Fig. 8a, a nanoisland geometry inspired by the experimental results in refs. 7,12 is considered. In addition to the previous simulation parameters, a magnetic defect is inserted in the top-left corner of the island, treated as a simulation cell with the magnetization frozen along the in-plane +x direction. The initial configuration including five skyrmions is evolved in time over a span of 20 ns, as shown in the snapshots in panel a. During the time evolution one skyrmion is pinned at the defect, whereas the others are moving rather freely. The resulting time average of the images reflects this behavior, as only a single spot is visible with a welldefined position in the upper left corner, and the rest of the quasiparticle features form a blurred trace mostly following the geometry of the island. In addition, the profile of the out-ofplane spin component is shown in Fig. 8b along the dashed arrow in the lower right panel of Fig. 8a. In this representation, the discrepancy between the measured profiles of the different skyrmions is even more pronounced. Where the localized skyrmion appears as a sharp dip almost reaching m z = −1, the superposition of the other four skyrmions leads to broader and shallower features, and therefore no clear indication of their positions.
Quasiparticle model. The localization of the quasiparticles discussed above may also be interpreted as the probability density of finding a skyrmion at a selected position over the complete simulation time. Such a probabilistic interpretation is capable of describing the different contrast levels observed for nominally equivalent skyrmions, and it proves to be valuable in the interpretation of scanning-probe experimental results where the characteristic measurement times are significantly longer than the time scale of the diffusive motion. In this section, a simplified quasiparticle model is introduced, motivated by the time-averaged data of the micromagnetic simulations, offering a time-efficient opportunity in predicting the complex pattern formation of the skyrmion probability distribution discussed above.
For the purpose of using a quasiparticle approach, the repulsive potential between pairs of skyrmions and between the skyrmion and the boundary needs to be quantified. Due to the Neumann boundary condition 40 at the edge of nanoislands, the DMI leads to a noncollinear spin texture, similar to a skyrmion 41 . Therefore, the interaction mechanism is the same between the skyrmions and between a skyrmion and the boundary, as has been studied in previous publications 29 Supplementary Fig. 1.
Starting from the modeled Sk-Sk and Sk-Bnd potentials quasiparticle simulations are executed, in order to compare the calculated probability densities with the time-averaged configurations obtained from full-fledged micromagnetic simulations. The main assumption for this endeavor is that neither the thermal fluctuations lead to the collapse, creation or escape of skyrmions, nor does a high skyrmion density lead to a merging of the skyrmions on the island as discussed in ref. 39 . A triangular setup containing two quasiparticles serves as the model system. The side of the triangle is chosen to be 30 nm, the same as in the case of the micromagnetic calculations shown in Figs. 2b and 5. The motion of the skyrmions inside the potential is calculated by modeling them as interacting random walkers, using an implementation of a Markov-chain Monte Carlo algorithm with Metropolis transition probabilities to study their diffusion, described in detail in the Methods section. This method offers a fast way of obtaining the cumulative probability density function of the positions of the two particles, which in turn can be compared to the time-averaged images in Figs. 4 and 5 obtained from micromagnetic simulations. In the following, a quantitative comparison between full micromagnetic simulations and Monte Carlo calculations will be presented to review the validity of the simplified quasiparticle approach. For this endeavor, we tracked the positions of two skyrmions in the triangular system from the micromagnetic simulation snapshots after every 20 ps and calculated the probability density distribution of the center coordinates. The resulting distribution is shown in Fig. 9a, along with the Monte Carlo result after 10 5 Monte Carlo steps in Fig. 9b. The largest difference between the two results arises because of the substantially less data for the micromagnetic calculation. As the histogram is calculated with only 10 3 snapshots compared with 10 5 Monte Carlo steps, this behavior is expected. Nevertheless, the characteristic features are similar. In both results the three peak densities are not radially symmetric, but are resembling the triangular boundary shape. Hence, this deformation is not an artifact of the Monte Carlo simulation.
In addition, the average distance between the density maxima may serve as a good quantity for comparing both methods. As Fig. 9a lacks the resolution for calculating the peak positions reliably, we start instead from the time-integrated magnetization pictures in Fig. 5. Based on the averaged result for two skyrmions, we extract again the peak positions and calculate subsequently the mean distance. The positions for the Monte Carlo calculations are at (11.0, 6.5) nm, (15.0, 13.5) nm, (19.0, 6.5) nm, the micromagnetic simulations deliver (11.1, 6.6) nm, (15.0, 13.5) nm, (19.0, 6.5) nm. This leads to an average distance of 8.0 nm for both the Monte Carlo calculations and the micromagnetic simulations, supporting the validity of the simplified quasiparticle approach in the investigated parameter space. The model is capable of predicting the skyrmion distribution in the long-time limit and thereby serves as a method for computationally efficient calculations for larger or more complex systems.
As an example, the probability density functions of the Monte Carlo simulations are shown in Fig. 10 after 10 5 Monte Carlo steps for different temperatures and sample sizes. Owing to the repulsion between the two skyrmions, equilibrium positions are found close to the vertices of the triangle, with energy barriers between these local energy minima. At low temperature (T = 1 K) and a small island size (30 nm edge length) in the upper left panel in Fig. 10, no transitions between the minima occur during the duration of the simulation, and two skyrmions may be observed in the obtained probability densities. As the temperature is increased to T = 15 K or T = 30 K in the first row of Fig. 10, the transition rate between the preferential positions is enhanced, and an additional "phantom skyrmion" appears in the probability density function, in remarkable agreement with the micromagnetic results as analyzed before. By enlarging the size of the island (second and third rows in Fig. 10 for 45 nm and 60 nm edge length), the energy surface becomes flatter, which in turn leads to higher transition probabilities between the energy minima at the same temperature. Therefore, raising the temperature or increasing the size of the system have similar effects on the resulting probability density functions for the skyrmions. Consequently, large island sizes or high temperatures result in completely delocalized skyrmions as indicated in the bottom right panel in Fig. 10.

Discussion
In this paper, we studied the temperature-driven diffusive motion of ensembles of magnetic skyrmions in finite magnetic islands. Based on experimentally determined system parameters for Furthermore, we proposed and developed a Monte Carlo quasiparticle model, which can be of great use to efficiently calculate the stochastic motion of larger systems of skyrmions in arbitrarily shaped geometries. According to the results, the time-averaged images of the micromagnetic simulations can be qualitatively well reproduced by the probability density distributions obtained from the simple quasiparticle model, as long as it is possible to treat the skyrmions as stable objects with lifetimes significantly longer than the simulation time at the given temperature. This method can be used to describe larger systems containing more skyrmions at various temperatures in a computationally efficient way, including model applications in storage technology like racetrack memory devices. Here, not only the stability of the information bit is important, but also its addressability, which is immediately affected by the diffusive motion of the skyrmions.

Methods
Micromagnetic simulations. Finite-temperature micromagnetic calculations, using the open-source, GPU-accelerated software package mumax3 40 , were performed to solve the LLG equation, for every simulation cell m i of the discretized magnetization vector field. Here, γ 0 = 1.76 × 10 11 (T −1 s −1 ) is the gyromagnetic ratio of an electron and α is the Gilbert damping parameter. The time-and space-dependent effective magnetic field, is composed of the external field B ext i ; exchange interaction field B exch i ¼ 2A exch =M sat Δm i , with A exch the exchange stiffness and M sat the saturation magnetization; demagnetizing field B d i ¼ M sat K ij Ã m j , where details on the calculation of the demagnetizing kernel K can be found in ref. 40 ; uniaxial magnetocrystalline anisotropy field B a i ¼ 2 K u =M sat m z e z , with K u the anisotropy constant; and the field generated by the DMI B dmi i ¼ 2D=M sat ð∂m z =∂x; ∂m z =∂y; À∂m x =∂x À ∂m y =∂yÞ T , with D the strength of the interfacial DMI. The simulation parameters were determined experimentally in ref. 7  In addition, an effective thermal field is included in Eq. (3) as where η is a random vector generated according to a standard normal distribution independently for each simulation cell and time step. k B is Boltzmann's constant, T is the temperature, ΔV is the size of the simulation cell and Δt is the time step of simulation. During the simulations, the cell size was set to ΔV = 0.3 × 0.3 × 0.4 nm 3 . The linear size of the cell is comparable to the lattice constant a = 0.271 nm of the Pd/Fe bilayer on Ir(111). This means that the micromagnetic simulations performed here should closely resemble the results of atomistic simulations where the skyrmions collapse when their size becomes comparable to the lattice constant 37,39 , and where the use of temperature-independent model parameters is justified.
The total skyrmion number N Sk in the simulations was calculated via Processing of the micromagnetic results. For the results shown in Fig. 6, the z component of the magnetization averaged up to time t is stored in grayscale pictures as a matrix, A ij (t) = (m z (x i , y j , t) + 1)/2, where x i and y j denote the position of the micromagnetic simulation cell. The distance induced by the Frobenius norm between the image averaged up to time t, A ij (t), and the image averaged over the whole simulation length τ = 20 ns, A ij (τ) is divided by the square root of the number of matrix elements and subtracted from 1, yielding the matching parameter This procedure gives a matching parameter of 1 for a perfect agreement of the compared pictures.
In order to determine the skyrmion radius as a function of time in Fig. 7, first the contour lines where the z component of the magnetization is zero are calculated for the skyrmion at each time step. The contour lines are discretized on N points in space, L i (t) = (x i (t), y i (t)), i ∈ (1, ..., N). Subsequently, the center of the skyrmion is calculated via which in turn is used to obtain the average skyrmion radius Finally, the radius is normalized with respect to the zero temperature radius r 0 obtained from the relaxed initial configuration of choice.  Quasiparticle simulations. For the calculation of the stochastic motion of skyrmions following the quasiparticle approach, the Metropolis algorithm 42 was utilized. With this method, the probability density function converges to the Boltzmann distribution determined by the energy functional of the system. For the triangular geometry, the potential surface was computed by taking the superposition of the potentials shown in Supplementary Fig. 1 from the three boundaries for every grid point of the chosen finite, rectangular mesh. The cell size was ΔV = 0.5 × 0.5 × 0.4 nm. The initial positions of the skyrmions were randomly chosen inside the confined structure, excluding the case in which both quasiparticles start from the same simulation cell. At each time step, a possible adjacent position is selected for each skyrmion simultaneously via pseudorandom numbers. If the energy of the attempted new state, including interaction between the skyrmions, is lower than the energy of the initial one, the skyrmion will move there. If not, the transition into this state happens with a probability of p(ΔE) = exp(−ΔE/k B T) following the Boltzmann distribution. Here ΔT = E 2 − E 1 is the energy difference between the final and the initial states, k B is Boltzmann's constant and T is the temperature. Subsequently, new attempted positions are generated and accepted or rejected as before until a fixed number of simulation steps is reached.

Data availability
Data that support the findings of this work are available from the corresponding author on request.

Code availability
For obtaining the results, we used the open-source, GPU-accelerated software package mumax3 40 for the micromagnetic simulations. The quasiparticle simulations were carried out via a Mathematica 43 code developed for this purpose. Input files are available from the corresponding author on request.