{'fraction_non_alphanumeric': 0.04224052250461451, 'fraction_numerical': 0.016505750390458613, 'mean_word_length': 4.411640414905878, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Brownian motion of magnetic domain walls and skyrmions, and their diffusion constants

Extended numerical simulations enable to ascertain the diffusive behavior at finite temperatures of chiral walls and skyrmions in ultra-thin model Co layers exhibiting symmetric - Heisenberg - as well as antisymmetric - Dzyaloshinskii-Moriya - exchange interactions. The Brownian motion of walls and skyrmions is shown to obey markedly different diffusion laws as a function of the damping parameter. Topology related skyrmion diffusion suppression with vanishing damping parameter, albeit already documented, is shown to be restricted to ultra-small skyrmion sizes or, equivalently, to ultra-low damping coefficients, possibly hampering observation.


I. INTRODUCTION
The prospect of ultra-small stable information bits in magnetic layers in presence of the Dzyaloshinskii-Moriya (DM) interaction [1] combined to the expectation of their minute current propagation [2], notably under spin-orbit torques [3], builds up a new paradigm in information technology. In stacks associating a metal with strong spin-orbit interactions e.g. Pt and a ferromagnetic metal such as Co, that may host isolated skyrmions, large domain wall velocities have also been forecast [4] and observed [5]. The DM interaction induces chiral magnetization textures, walls or skyrmions, that prove little prone to transformations of their internal structure, hence their extended stability and mobility.
In order, however, to achieve low propagation currents, steps will need to be taken towards a reduction of wall-or skyrmion-pinning. Recent experimental studies indicate that skyrmions fail to propagate for currents below a threshold roughly equal to 2 10 11 Am −2 for [Pt/Co/Ta] n and [Pt/CoFeB/MgO] n multilayers [6], or 2.5 10 11 Am −2 for [Pt/(Ni/Co/Ni)/Au/(Ni/Co/Ni)/Pt] symmetrical bilayers [7]. Only in one seldom instance did the threshold current fall down to about 2.5 10 10 Am −2 for a [Ta/CoFeB/TaO] stack, still probably, however, one order of magnitude higher than currents referred to in simulation work applying to perfect samples [8].
In a wall within a Co stripe 50 nm wide, 3 nm thick, the number of spins remains large, typically 2 16 for a 5 nm wide wall. A skyrmion within a Co monolayer (ML) over Pt or Ir, on the other hand, contains a mere 250 spins, say 2 8 . Assuming that a sizeable reduction of pinning might somehow be achieved, then a tiny structure such as a skyrmion is anticipated to become sensitive, if not extremely sensitive, to thermal fluctuations.
In this work, we show, on the basis of extended numerical simulations, that both chiral walls and skyrmions within ferromagnets obey a diffusion law in their Brownian motion at finite temperature [9,10]. The diffusion * jacques.miltat@u-psud.fr law is shown to be valid over a broad range of damping parameter values. The thermal diffusion of domain walls seems to have attracted very little attention, except for walls in 1D, double potential, structurally unstable, lattices [11], a source of direct inspiration for the title of this contribution. Chiral magnetic domain walls are found below to behave classically with a mobility inversely proportional to the damping parameter. As shown earlier [12,13], such is not the case for skyrmions, a behavior shared by magnetic vortices [14]. Vortices and skyrmions in ferromagnetic materials are both characterized by a definite topological signature. In contradistinction, skyrmions in antiferromagnetic compounds are characterized by opposite sign spin textures on each sublattice, with, as a result, a classical, wall-like, dependence of their diffusion constant [15]. Lastly, ferrimagnets do display reduced skyrmion Hall angles [16], most likely conducive to modified diffusion properties.

II. DOMAIN WALL DIFFUSION
We examine here, within the micromagnetic framework, the Langevin dynamics of an isolated domain wall within a ferromagnetic stripe with thickness t S , width w S and finite length L (see Fig. 1). The wall is located at mid-position along the stripe at time t = 0. Thermal noise is introduced via a stochastic field H Rd uncorrelated in space, time and component-wise, with zero mean and variance η proportional to the Gilbert damping parameter α and temperature T [17] : where, k B is Boltzmann constant, µ 0 and γ 0 are the vacuum permability and gyromagnetic ratio, respectively, M S the saturation magnetization. Written as such, the functions δ( r − r ) and δ(t − t ) have the dimension of reciprocal volume and time, respectively. Applied to numerical simulations, the variance of the stochastic field becomes η = 2kBT γ0µ0MSV dt α, where V is the computation cell volume and dt the integration time step.

A. Simulation results
The full set of numerical simulations has been performed by means of an in-house code ported to graphical processing units (GPU's). Double precision has been used throughout and the GPU-specific version of the "Mersenne twister" [18] served as a source of longsequence pseudo-random numbers generator.
Material parameters have been chosen such as to mimic a 3-ML Co layer (thickness t S = 0.6 nm) on top of Pt with an exchange constant equal to A = 10 −11 J/m, a M s = 1.09 10 6 A/m saturation magnetization, a K u = 1.25 10 6 J/m 3 uniaxial anisotropy constant allowing for a perpendicular easy magnetization axis within domains, and a moderate-to-high DM interaction (DMI) constant D DM = 2 mJ/m 2 . In order to temper the neglect of short wavelength excitations [19], the cell size has been kept down to L x = L y = 1 nm, whilst L z = t S = 0.6 nm. The stripe length has been kept fixed at L = 1 µm, a value compatible with wall excursions within the explored temperature range. The latter has, for reasons to be made clear later, been restricted to ≈ 1/3 of the presumed Curie temperature for this model Co layer. Finally, the integration time constant, also the fluctuating field refresh time constant, has been set to dt = 25 fs.
As shown by the snapshot displayed in Fig. 1b, the wall may acquire some (moderate) curvature and/or slanting during its Brownian motion. Because wall diffusion is treated here as a 1D problem, the wall position q is de- fined as the average position owing to : where, i and j are the computation cell indices, N x and N y the number of cells along the length and the width of the stripe, respectively, m z L is the fluctuations averaged value of the z magnetization component far left of the domain wall, m z R the average value of m z far right.
Regardless of sign, m z R and m z L are expected to be equal in the absence of any H z field. Fig. 2 displays the position as a function of time of a wall within a w S = 100 nm wide stripe immersed in a T = 77 K temperature bath. A 2 ns physical time window has been extracted from a simulation set to run for 1.5 µs. The figure shows short term wall position fluctuations superimposed onto longer time diffusion. According to Einstein's theory of Brownian motion [9], the probability P (x, t) of finding a particle at position x at time t obeys the classical diffusion equation So does the raw probability of finding a (stiff) wall in a narrow stripe at position q after a time interval ∆t, as shown in Fig. 3 (see Fig. 2 for variable definition). It ought to be mentioned that the average wall displacement q(∆t) is always equal to 0, with an excellent accuracy, provided the overall computation time is large enough. The fit to a normal distribution proves rather satisfactory, with, however, as seen in Fig. 3, a slightly increasing skewness in the distributions as a function of increasing ∆t. Skewness, however, 1) remains moderate up to ∆t values typically equal to 5 − 10 ns, 2) is seen to reverse sign with time interval (compare Fig. 3b and c), excluding intrinsic biasing. The distributions standard deviation is clearly seen to increase with increasing ∆t. Alternatively, one may represent the variance q 2 ( q = 0) as a function of the time interval ∆t : if diffusion applies, then a linear dependence is expected, with a 2D slope for a one-dimensional diffusion. Fig.4a shows, for various temperatures, that a linear law is indeed observed. Lastly, as shown in Fig.4b, the diffusion constant increases linearly with increasing temperature. The error bars measuring the departure from strict linearity in Fig.4a remain limited in extent. For the stripe width and damping parameter considered here (w S = 100 nm, α = 0.5), the ratio of diffusion constant to temperature is found to amount to D/T = 0.187 nm 2 ns −1 K −1 .

B. Wall diffusion constant (analytical)
Thiele's equation [20] states that a magnetic texture moves at constant velocity v provided the equilibrium of 3 forces be satisfied: where, F is the applied force, F G = G× v is the gyrotropic force, G the gyrovector, F D = αD v the dissipation force, D the dissipation dyadic. For the DMI hardened Néel wall considered here : G = 0. For a 1D wall, the Thiele equation simply reads : The calculation proceeds in two steps, first evaluate the force, hence, according to Eqn.4, the velocity autocorrelation functions, then integrate vs time in order to derive q 2 . The force, per definition, is equal to minus the partial derivative of the energy E w.r.t. the displacement q, As noticed earlier [14], since the random field noise is "multiplicative" [17], moving the magnetization vector out of the average brackets is, strictly speaking, not allowed, unless considering the magnetization vector to only marginally differ from its orientation and modulus in the absence of fluctuations (the so-called "low" noise limit [14]): If due account is being taken of the fully uncorrelated character of the thermal field (Eqn.1), the force autocorrelation function becomes: The velocity auto-correlation function follows from Eqn.4. Lastly, time integration (q(t) = t 0 v x (t )dt ) yields : In order to relate the diffusion constant to a more directly recognizable wall mobility, D xx may be expanded as : where, ∆ T has been called the Thiele wall width (implicitly defined in [21]). D may thus be expressed as : thus, proportional to the wall mobility γ 0 ∆ T /α. A directly comparable result may be obtained after constructing a full Langevin equation from the (q, φ) equations of domain wall motion (Slonczewski's equations [22]), where φ is the azimuthal magnetization angle in the wall mid-plane. In this context, the wall mobility is µ W = γ 0 ∆/α, where ∆ is the usual wall width, incidentally equal to the Thiele wall width in the case of a pure Bloch wall. The Langevin equation [10] here reads: where, m D is Döring's wall mass density (kg/m 2 ): an expression valid in the limit | D DM | K Eff = K u − 1 2 µ 0 M 2 s . Note that the DMI constant D DM explicitly enters the expression of the wall mass, as a consequence of the wall structure stiffening by DMI. In the stationary regime, q 2 is proportional to time t and the wall diffusion constant exactly matches Eqn.10, after substitution of ∆ T by ∆. Finally, the characteristic time for the establishment of stationary motion is: For the parameters of our model 3-ML Co layer on top of Pt, Döring's mass density is equal to ∼ 3 10 −8 kg/m 2 for α = 0.5, and the characteristic time amounts to t 0 25 ps. Still for α = 0.5, w S = 100 nm and t S = 0.6 nm, D/T amounts to 0.153 nm 2 ns −1 K −1 for ∆ T = 4.13 nm, i.e. the value computed from a properly converged wall profile at T = 0. The relative difference between simulation and theoretical values is found to be of the order of ≈ 20%. Owing to Eqn.10, D is expected to prove inversely proportional to both the stripe width w S and the Gilbert damping parameter α, a behavior confirmed by simulations. Fig.5a displays the computed values of the diffusion coefficient as a function of temperature with the stripe width as a parameter, whilst Fig.5b states the linear behavior of D vs w S −1 . The slope proves, however, some 13.5% higher than anticipated from Eqn.10. Lastly, the 1/α dependence is verified in Fig.6 showing the computed variation of D vs temperature with α as a parameter for a narrow stripe (w S = 25 nm) as well as the corresponding α dependence of D/T . The dotted line represents Eqn.10 without any adjusting parameter. The relative difference between simulation data and theoretical expectation is beyond, say α = 0.25, seen to grow with increasing α but also appears to be smaller for a narrow stripe as compared to wider tracks.
Altogether, simulation results only moderately depart from theoretical predictions. The Brownian motion of a DMI-stiffened wall in a track clearly proves diffusive. The diffusion constant is classically proportional to the wall mobility and inversely proportional to the damping parameter. Unsurprisingly, the smaller the track width, the larger the diffusion constant. In order to provide an order of magnitude, the diffusion induced displacement expectation, √ 2D∆t, for a wall sitting in a 100 nm-wide, pinning-free, track for 25 ns at T = 300 K proves essentially equal to ± the stripe width.

A. Simulation results
In order to monitor the Brownian motion of an isolated skyrmion, rather than micromagnetics, it is preferred to simulate the thermal agitation of classical spins, s (| s |= 1), on a triangular lattice. Lattice effects and frequency cutoffs in thermal excitations are thus avoided. Such simulations have already been used e.g. for the determination of the barrier to collapse of an isolated skyrmion [25,26]. The parame- Fast Fourier Transforms implementation ensues from the decomposition of the triangular lattice into two rectangular sublattices, at the expense of a multiplication of the number of dipole-dipole interaction coefficients. Lastly, the base time step, also the stochastic field refresh time, has been given a low value in view of the small atomic volume, namely dt = 2.5 fs for α ≥ 0.1, dt = 1 fs below. Time steps that small may be deemed little compatible with the white thermal noise hypothesis [17]. They are in fact dictated by the requirement for numerical stability, primarily w.r.t. exchange interactions. Fig.7 is a snapshot of an isolated skyrmion in the model Co ML with a temperature raised to 12.5 K. The skyrmion is at the center of a 200 at. u.-i.e. ≈ 50 nm-size square computation window, that contains 46400 spins and is allowed to move with the diffusing skyrmion. Doing so alleviates the computation load without restricting the path followed by the skyrmion. Free boundary conditions (BC's) apply. The window, however, proves sufficiently large to render the confining potential created by BC's ineffective. The skyrmion position as a function of time is defined simply as the (iso)barycenter of the contiguous lattice site positions x(k), y(k), where s z ≥ 0.5: where, k is the lattice site index, K the number of lattice sites satisfying the above condition. Such a definition proves robust vs thermal disorder such as displayed in Fig. 7. Similarly to the case of wall diffusion, we analyze first the distributions of the displacement components q x , q y . The event statistics for each value of the time interval is clearly gaussian (see Fig.9). However, the noise in the distributions appears larger when compared to the wall case. It also increases faster with ∆t. On the other hand, the raw probabilities for q 2 x and q 2 y barely differ as anticipated from a random process. The behavior of q 2 (q 2 = q 2 x + q 2 y ) vs ∆t is displayed in Fig.10a. The range of accessible temperatures is governed by the thermal stability of the tiny skyrmion within a Co ML: with a lifetime of 1 µs at 77 K [25][26][27][28], temperatures have been confined to a ≤ 50 K range. When compared to the wall case (Fig.4a), the linear dependence of q 2 with respect to ∆t appears less satisfactory, although, over all cases examined, the curves do not display a single curvature, but rather meander gently around a straight line. The slope is defined as the slope of the linear regression either for time intervals between 0.25 and 2.5 ns (thick line segments in Fig.10a) or for the full range 0 to 5 ns (dashed lines). Then, the ratio of the diffusion constant to temperature, D/T , for an isolated skyrmion within the model Co ML considered here is equal to 0.250 and 0.249 nm 2 ns −1 K −1 , respectively, for α = 0.5 (see Fig.10b). The difference proves marginal. Lastly, error bars appear even narrower than in the wall case.

B. Skyrmion diffusion constant (analytical)
The gyrovector G in Thiele's equation (Eqn.3) has in the case of a skyrmion or a vortex, and in many other instances such as lines within walls, a single non-zero component, here G z . Thiele's equation, in components form, reads: Because of the revolution symmetry of a skyrmion at rest, D xy or D yx may safely be neglected and D yy = D xx . Accordingly, the velocities may be expressed as: Similarly to the stochastic field, the force components are necessarily uncorrelated. The velocity autocorrelation functions may now be obtained following the same lines as in the wall case, yielding, in the low noise approximation: The average values of the displacements squared, q 2 x and q 2 y follow from time integration: As shown previously [12,13], the diffusion constant for a skyrmion thus reads: The following relations do apply: Relation (19) implies a peculiar damping constant dependence with, assuming for the time being D and G to have comparable values, a gradual drop to zero of the diffusion constant with decreasing α (α ≤ 1), termed "diffusion suppression by G" by C. Schütte et al. [12]. Diffusion suppression is actually not a complete surprise since, for electrons in a magnetic field, a similar effect is leading to the classical magnetoresistance. A similar dependence D(α) is also expected for a vortex. Boundary conditions, however, add complexity to vortex diffusion. What nevertheless remains, is a linear dependence of D vs α [14], namely, diffusion suppression.
The classical expressions for G z and D xx valid for a magnetization continuum need to be adapted when dealing with discrete spins. We obtain: where, µ At is the moment per atom. The dimensionless product γ0S At µ0µ At G z (Eqn.21), where S At is the surface per atom, amounts to 4π, irrespective of the skyrmion size in a perfect material at T = 0. Stated otherwise, the skyrmion number is 1 [29]. In the Belavin-Polyakov profile limit [30], the dimentionless product γ0S At µ0µ At D xx (Eqn.21) also amounts to 4π. In this limit, D is proportional to α/(1 + α 2 ). D xx increases with skyrmion radius beyond the Belavin-Polyakov profile limit (see supplementary material in [7]). For a skyrmion at rest in the model Co ML considered here, D = D xx ≈ 14.5 µ 0 µ At /(γ 0 S At ). For that value of D xx , and for the parameters used in the simulations, D/T , the ratio of the theoretical skyrmion diffusion constant to temperature, is equal 0.234 nm 2 ns −1 K −1 , for α = 0.5 (S At = a 2 √ 3/2), to be compared to the 0.250 value extracted from simulations. More generally, Fig.11 compares numerical D/T values calculated for a broad spectrum of α values with theoretical expectations for D = 14.5 µ 0 µ At /(γ 0 S At ) and in the Belavin-Poliakov limit. The average difference between analytical and simulation results is, in the α = (0, 1) interval, seen to be of the order of 15%.

IV. DISCUSSION
In the present study of thermal diffusion characteristics, satisfactory agreement between simulations and theory has been attained for DMI stiffened magnetic textures, be it walls in narrow tracks or skyrmions. The α dependence of the diffusion constants has been thoroughly investigated, with, as a result, a confirmation of Brownian motion suppression in the presence of a nonzero gyrovector or, equivalently, a topological signature. The theory starts with the Thiele relation applying to a texture moving under rigid translation at constant velocity. Furthermore, the chosen values of the components of the dissipation dyadic, are those valid for textures at rest, at T = 0. The α dependence of the diffusion constants clearly survives these approximations. And, yet, a wall within a narrow stripe or a skyrmion in an ultra-thin magnetic layer are deformable textures, as obvious from Figs.1,7. Simulations, on the other hand, rely on the pioneering analysis of Brownian motion, here meaning magnetization/spin orientation fluctuations [17], within a particle small enough to prove uniformly magnetized and then extend the analysis to ultra-small computation cell volumes down to the single spin. Both approaches rely on the hypothesis of a white -uncorrelated-noise at finite temperature.
The discussion of results is organized in two parts. In the first, results are analyzed in terms of a sole action of structure plasticity on the diagonal elements of the dissipation dyadic. In the second, we envisage, without further justification, how the present results are amended if, in the diffusion constants of walls and skyrmions (Eqns. 8 and 19), the gyrotropic and dissipation terms are re- placed by their time average as deduced from simulations.

A. Size effects
The integral definition of wall position adopted in this work (Eqn.2) allows for a 1D treatment of wall diffusion, thus ignoring any diffusion characteristics potentially associated with wall swelling, tilting, curving or meandering. Additional information is, however, available in the case of skyrmions. We concentrate here on the number, n, of spins within the skyrmion satisfying the condition s z ≥ 0.5, and its fluctuations as a function of time. The surface of the skyrmion is nS At and its equivalent radius, r Eq , is defined by r 2 Eq = nS At /π. The skyrmion radius r Eq is found to fluctuate with time around its average value, according to a gaussian distribution that depends on temperature, but becomes independent of the autocorrelation time interval beyond ≈ 25 ps. The power spectrum of the time series r Eq (t), shown in Fig.12a, excludes the existence of a significant power surge around the fundamental breathing mode frequency of the skyrmion (≈ 39.3 GHz for the present model Co ML) [31]. The skyrmion radius as defined from the discrete n distribution is thus subject to white noise. The average radius r Eq , on the other hand, varies significantly with temperature, increasing from ≈ 1.6 nm to 2.4 nm when the temperature is increased from 4.2 K to 50 K (Fig.12b) and the diagonal element of the dissipation dyadic is expected to increase with increasing skyrmion radius [3,7].
Owing to relations (19,21), the maximum of D(α) is found for α = G z /D xx = G/D. For α < G/D, resp. α > G/D, D increases, resp. decreases, with D, hence the relative positions of the blue and black continuous curves in Fig.11. At maximum, D is independent of D and amounts to k B T γ0S At that the discrepancy between numerical and analytical D values around α = 1 may not be relaxed by a sole variation of D. On the other hand, allowing D to increase with skyrmion radius, itself a function of temperature, leads to an increase (decrease) of the diffusion coefficient for α < G/D (α > G/D). Likely more important is the reduction, as a function of skyrmion size, of the α window where diffusion suppression is expected. If including the (R/∆ + ∆/R) dependence of D xx (see supplementary material in [7]; ∆ is the wall width and R the skyrmion radius), the skyrmion diffusion constant may be expressed as: The general shape of function f (α, R/∆) is shown in Fig.13a. The maximum of f (α, R/∆) is equal to 1 for all values of α and R/∆. The crest line Rα = ∆ is seen to divide the parameter space into two regions (see Fig.13b), a region close to the axes where ∂D/∂α > 0, i.e. the region of diffusion suppression, from the much wider region where ∂D/∂α < 0, that is, the region of wall-like behavior for skyrmion diffusion. Clearly, the α window for diffusion suppression decreases dramatically with increasing skyrmion size R/∆. A first observation of skyrmion Brownian motion at a video recording time scale (25 ms) may be found in the Supplementary Material of Ref. [32]. Skyrmions are here unusually large and most likely escape the diffusion suppression window (α < 0.02 for R/∆ = 50). Combining skyrmion thermal stability with general observability and damping parameter tailoring may, as a matter of fact, well prove extremely challenging for the observation of topology related diffusion suppression.

B. Time averaging
One certainly expects from the simulation model a fair prediction of the average magnetization M z or S z vs temperature T , at least for temperatures substantially lower than the Curie temperature T C . Fig.14 shows the variation of M z /M z (T = 0) or S z /S z (T = 0) with temperature for the two model magnetic layers of this work. Although simulation results do not compare unfavorably with published experimental data [33][34][35], where, typically, the Curie temperature amounts to ≈ 150K for 1 ML, and proves larger than 300K for thicknesses above 2 ML, a more detailed analysis, potentially including disorder, ought to be performed.
Let us now, without further justification, substitute in the expression of the skyrmion diffusion coefficient time averaged values of G and D, owing to relations (23). Keeping in mind the geometrical meaning of G VF z , the dimensionless vector function in G, G z is anticipated to be a sole function of s z . Inversely, D VF xx , the (dimensionless) vector function in D xx , a definite positive quantity, steadily increases with thermal disorder. It is even found to be proportional to temperature (not shown). Its time averaged value for the sole skyrmion may only be obtained by subtraction of values computed in the presence and absence of the skyrmion.
For the skyrmion in our model Co monolayer, D VF xx is found to increase moderately with temperature (see Fig.14), a result also anticipated from an increase with temperature of the skyrmion radius. Besides, both G z and D xx are expected to decrease with temperature due to their proportionality to s z . D xx is thus subject to two competing effects of temperature T . Present evidence, however, points at a dominating influence of s z (T ) .

V. SUMMARY AND OUTLOOK
Summarizing, it has been shown that the Brownian motion of chiral walls and skyrmions in DMI materials obeys diffusion equations with markedly different damping parameter (α) dependence. Although not a new result, skyrmions Brownian motion suppression with decreasing α (α < G/D) is substantiated by a wide exploration of the damping parameter space. The observation of this astonishing topological property might, however, be hampered by the restriction to ultra-small skyrmion sizes or ultra-low α values. The discrepancy (up to 20%) between simulation results and theoretical expectations could be reduced by the introduction of time averaged values for the gyrotropic and dissipation contributions to the analytical diffusion coefficients in the "low" noise limit, at the expense of a tiny upwards curvature in the D(T ) curves. A strong theoretical justification for doing so remains, however, lacking at this stage.
In this work, the sample has been assumed to be perfect, i.e. devoid of spatial variations of the magnetic properties, even though the lifting of such a restriction is anticipated to prove mandatory for a proper description of experiments. Diffusion in the presence of disorder has been theoretically studied for a number of disorder and random walk types [36,37]. Generally, disorder changes the linear growth with time of the position variance into a power law, a behavior called superdiffusion if the exponent is larger than 1 and subdiffusion if smaller. For instance, if the skyrmion motion in a disordered system may be mapped onto a 2D random walk with an onsite residence time τ , probability ∝ τ −(1+µ) (µ < 1), then the diffusion exponent will be µ, meaning subdiffusion. Besides, choosing a physically realistic disorder model for a Co monolayer might well prove equally arduous [38]. Altogether, skyrmion diffusion in the presence of disorder has been left out for future work.