Skyrmion Gas Manipulation for Probabilistic Computing

TL;DR

This paper introduces a skyrmion gas-based device that leverages thermally driven dynamics for probabilistic computing, demonstrating energy efficiency and potential as a bio-inspired computational element.

Contribution

It proposes a novel skyrmion gas device for probabilistic computing, modeling ensemble dynamics and demonstrating a proof-of-concept with low energy and area requirements.

Findings

Energy-efficient (~μW) device demonstrated

Reshuffles random signals into uncorrelated copies

Potential application in stochastic computing and bio-inspired systems

Abstract

The topologically protected magnetic spin configurations known as skyrmions offer promising applications due to their stability, mobility and localization. In this work, we emphasize how to leverage the thermally driven dynamics of an ensemble of such particles to perform computing tasks. We propose a device employing a skyrmion gas to reshuffle a random signal into an uncorrelated copy of itself. This is demonstrated by modelling the ensemble dynamics in a collective coordinate approach where skyrmion-skyrmion and skyrmion-boundary interactions are accounted for phenomenologically. Our numerical results are used to develop a proof-of-concept for an energy efficient ($\sim\mu\mathrm{W}$) device with a low area imprint ($\sim\mu\mathrm{m}^2$). Whereas its immediate application to stochastic computing circuit designs will be made apparent, we argue that its basic functionality, reminiscent of an integrate-and-fire neuron, qualifies it as a novel bio-inspired building block.

Figures

Tables

Table 1. Table 1: Material Parameters

	$M_S$ ($MA/m$)	$J_{ex}$ ($fJ/m$)	$D$ ($mJ/m^2$)	$K$ ($MJ/m^3$)	$B$ ($mT$)	$R_{\mathrm{skx}}$ ($nm$)
S1	1.4	27.5	2.05	1.45	0	23.2
S2	1.45	27.5	2.05	1.45	20	19.1
S3	1.4	27.5	2.6	1.45	20	30.3
S4	1.4	27.5	2.6	1.60	-15	14.9
S5	1.5	25.0	2.4	1.60	30	27.9
S6	1.0	5.0	2.05	0.9	150	11.9
S7	1.2	8.0	2.1	0.9	150	14.9

Table 2. Table 2: Skx-Skx Repulsion

	$\left|F\right|\propto \exp \left\{-\frac{a_2{d}^2+a_{1}d+a_0}{1+d}\right\}$
	$a_2$	$a_1$	$a_0$
S1	1.245	1.249	-3.449
S2	2.709	-5.643	0.964
S3	2.488	-1.899	-2.676
S4	1.13	2.171	-6.623
S5	2.219	-1.159	-3.358
S6	0.059	10.757	-16.738
S7	0.766	5.452	-11.894

Equations

$${\hat{\mathbf{G}}}^{\alpha }\cdot \dot{\mathbf{x}}=\mu {\hat{\mathbf{G}}}^{\eta }\cdot {\mathbf{v}}_s+{\mathbf{F}}_{\mathrm{Th}}+{\mathbf{F}}_{\mathrm{Ext}}$$

$${\hat{G}}_{ij}^{\sigma }=G{ϵ}_{ij}+\sigma D{\delta }_{ij}$$

$$G$$

$$D$$

$$⟨{\mathbf{F}}_{\mathrm{Th}}⟩$$

$$⟨F_{\mathrm{Th},i}(t)F_{\mathrm{Th},j}(t^{\prime })⟩$$

$$⟨|\mathbf{x}{|}^2⟩(t)=2k_{B}T\frac{\alpha D}{G^2+(\alpha D{)}^2}t.$$

$${\mathbf{F}}_{\mathrm{Ext}}={\mathbf{F}}_S+{\mathbf{F}}_{\mathrm{B}},$$

$${\mathbf{F}}_S$$

$${\mathbf{F}}_{S-S}(\mathbf{d})$$

$$|\mathbf{F}{|}_{\mathrm{B}}\simeq \exp \left[-\sum _0^6{b}_j{d}^j\right],$$

$${\rho }_{X,Y}=\frac{\mathrm{cov}(X,Y)}{{\sigma }_X{\sigma }_Y},$$

$$⟨|\mathbf{x}{|}^2⟩(t)=2\frac{k_{B}T}{\gamma M}\left(\frac{\alpha D}{G^2+(\alpha D{)}^2}\right)t,$$

$$\frac{\rho }{1-\rho }\propto \frac{l_0}{\sqrt{⟨|\mathbf{x}{|}^2⟩(t_{\mathrm{exit}})}}\propto \nu \sqrt{\frac{\gamma M}{k_B}\left(\frac{G^2+\alpha D^2}{\alpha D}\right)\frac{J^3}{LT}},$$

$$E_{\mathbf{M}}$$

$$+$$

$$E_{\mathbf{M}}={\int }_{-\infty }^{\infty }\mathrm{d}x\left[|{\partial }_x\theta {|}^2+Q^2{\sin }^2\theta \right]-\kappa \pi ,$$

$$\theta (x)=2\mathrm{atan}\left[e^{-Qx}\right].$$

$$\theta (x)={\theta }_{+}(x)+{\theta }_{-}(x)=2\mathrm{atan}\left[\eta e^{Qx}\right]+2\mathrm{atan}\left[\eta e^{-Qx}\right],$$

$$E_{\theta }=E_{{\theta }_{+}}+E_{{\theta }_{-}}+E_{{\theta }_{+},{\theta }_{-}}^{\mathrm{int}}(\delta )$$

$$E_{{\theta }_{+},{\theta }_{-}}^{\mathrm{int}}(\delta )={\int }_{-\infty }^{\infty }\mathrm{d}x[{\partial }_x{\theta }_{+}{\partial }_x{\theta }_{-}+\frac{1}{2}\sin 2{\theta }_{+}\sin 2{\theta }_{-}$$

$$-2{\sin }^2{\theta }_{+}{\sin }^2{\theta }_{-}]$$

$$F_{{\theta }_{+},{\theta }_{-}}(\delta )$$

$$I_{Q,R,S,T}$$

Full text

1 Introduction

Magnetic skyrmions promise unique opportunities for the processing, storage and transfer of information at the intersection of both spintronics and nanoelectronics [1, 2, 3, 4, 5, 6, 7]. The experimental verification of skyrmion spin configurations has been a forefront research topic in magnetism over the past decade. Initially observed in bulk non-centrosymmetric crystal lattices [8, 9, 10, 11, 12], they have been more recently stabilized in ultrathin ferromagnetic nanostructures strongly affected by the Dzyaloshinski-Moriya interaction (DMI) resulting from coupling to a heavy metal substrate [13, 14]. In such materials, skyrmions are exceptionally stable spin textures capable of enduring room-temperature environments [15, 16, 17, 18] and manipulable at small current densities ($10^{6}-10^{11}\,\mathrm{A}/\mathrm{m}^{2}$) [19, 20, 17, 21, 22, 23]. Furthermore, their topological resilience [24, 2] allows them to avoid pinning onto defects [25, 26] and guarantees particle number conservation under a wide range of operating conditions [27, 28].

The appeal of skyrmions is so wide that it has defined a field of its own, skyrmionics, which refers to the emerging technologies based on magnetic skyrmions as information carriers. In an effort to push for skyrmion-based electronics, challenges ranging from their creation and annihilation [2, 29, 30, 31], the conversion of their topological properties [5, 32], as well as their efficient transmission and read-out [2, 4, 3, 33, 6] are being tackled and solved. Recently, proposals have exploited their nano-metric size and high mobility [34, 26, 20, 1] in the context of skyrmion-based racetrack memories [3, 29, 35], logic gates [5], voltage-gated transistors [36] and synaptic devices [37]. All such applications typically constrain the skyrmion motion onto 1-dimensional tracks reminiscent of previous device designs seeking to employ domain wall motion. Very little has been proposed to take advantage of the analogy between skyrmions and a generic free moving particle. Particularly, the full 2-dimensional freedom of skyrmion motion has not been leveraged either at a functional level or as a technique to reduce the surface imprint of proposed devices.

All such applications effectively employ geometries where one or more skyrmions are injected, manipulated and read out (see Fig. 1). Currently, the standard approach for studying their behavior is to employ micromagnetic simulations capable of modeling the nonlinear evolution of the structure’s magnetic texture over time [38]. Whereas micromagnetic techniques can be very precise, the execution time of the finite-difference method they employ scales rapidly with the volume of the structure itself [39]. As such, resolving the behavior of large structures and circuits over long time periods becomes impractical even on the nanosecond scale. It is therefore important to have an efficient method for simulating the behavior of skyrmion ensembles which scales with the number of skyrmions instead. Such a technique should be capable of capturing both the influence of magnetic dipole effects as well as boundary interactions on the ensemble dynamics.

At a device level, skyrmionic applications are not currently expected to compete significantly with present-day CMOS implementations of boolean logic circuits whose strengths lie in the accuracy and precision of bitwise operations. However, limits in the exponential miniaturization of transistor-based architectures exposes bottlenecks both with regards to speed as well as their power consumption. This is particularly relevant when requiring for hardware to enable software capable of performing fast inference over large amounts of data. As an example, much attention has been directed towards the development of artificial neural networks for the execution of deep-learning techniques [40]. Such solutions typically seek to draw inspiration from the probabilistic and massively parallel nature of biological information processing [41] whose virtues have been explored for the past half century [42]. In this context, we propose that greater focus should instead be given to the use of skyrmions as basic elements in the design of probabilistic device architectures. By leveraging skyrmion number conservation, thermal susceptibility and ultra-low power transport properties, we believe that a fertile route towards disruptive applications can be achieved. Along this line, our work presents two novel devices implementing basic building blocks for probabilistic computing.

Starting from the collective coordinate dynamical theory introduced by Thiele [43], this paper develops an N-body framework capable of modelling the dynamics of skyrmion ensembles within their geometries and successively demonstrate how they may be employed to perform probabilistic computing. To this aim we first, in Section 2, motivate the deep technological challenges inspiring this work by giving a general background on stochastic computing and justifying the need for a compact and energy efficient stochastic signal reshuffler which we propose employing skyrmions as information carriers. In Section 3, we use micromagnetic simulations to characterize the repulsive skyrmion-skyrmion and skyrmion-boundary interactions. These then enter phenomenologically into the N-body Thiele model where the ensemble dynamics become numerically solvable out to long timescales. These advantages permit us in Section 4 to model the Skyrmion Reshuffler leveraging the full freedom of skyrmionic motion. We characterize the device functioning performance and expected energy consumption. Finally, in Section 5, we argue how our skyrmion reshuffler, subject to minor modifications, can serve a much more general purpose by effectively working as an analog integrate-and-fire neuron [44]. These results lay the groundwork for the application of skyrmionic devices as bio-inspired building blocks for non-conventional computing strategies.

2 Motivation: Stochastic Computing

2.1 Motivation

As feature sizes in semiconductor systems are scaled down to the point where their deterministic behavior cannot be guaranteed, the necessity for error-correcting codes [45] and hardware redundancies [46, 47] have become indispensable for reliable data processing. In this context, stochastic computing seeks to embrace probabilistic computing elements to bypass these issues by achieving tremendous gains in signal processing efficiency at the cost of degradation in computational precision [48, 49, 50, 51, 52]. Instead of defining operations as precise bitwise manipulations of stored finite binary numbers, one instead aims to encode numerical values as the probability of seeing a $1$ or [math] in a random sequence of binary bits (bitstream, as in Fig. 2 [53]). This probability, known as the p-value, can be also thought as the statistical ratio of up-time to signal length in a random telegraph noise signal. CMOS circuits consisting of cascaded logic gates can then be shown to perform mathematical operations on the p-values of distint random streams 111We will refer to telegraph noise signals and bitstreams interchangeably throughout this manuscript in real-time. Complex circuits have already been proposed to perform operations such as square-rooting [54], polynomial arithmetic [55, 56, 57], matrix operations [58] as well as the "tansig" transform function employed in neural-networks [59]. The output of such operations will result in an output bitstream whose p-value can be read out and possibly operated on successively. Last but not least, stochastic computing circuitry is inherently tolerant to faults such as soft bit-flips [60, 61]. Consider the act of randomly flipping one binary bit from a sequence of N bits (ex: $010001101\to 010011101$): the binary number encoded changes drastically while the p-value encoded is altered only by a factor of $o(1/N)$.

A very attractive property of stochastic computing is that it enables low-energy cost implementations of arithmetic operations using standard logic elements. As an example, consider the effect of feeding two indefinitely long random streams with p-values $p_{A}$ and $p_{B}$ into an AND-gate (see Fig. 2a). The AND-gate’s properties will produce an output signal with p-value $p_{C}=p_{A}\cdot p_{B}$, thus implementing the stochastic computing equivalent of multiplication. Whereas the result of the multiplication is produced in real-time as the input signals are fed to the AND-gate, the accuracy of the result depends on how much of the output signal one wishes to sample. This trade-off of speed vs accuracy is central in the development of integrated circuits performing stochastic computing. A detailed description of the initial encoding and tuning/setting of input p-values goes beyond the scope of this work. We note however that, aside from the traditional shift-register based techniques [62, 63], a compact and energy efficient approach for generating suitable stochastic signals can be achieved through current-biased superparamagnetic magnetic tunnel junctions [64, 65, 66, 67] whose individual use in neuromorphic computing architectures has also been recently explored [68, 69, 70]. This ties in well with the spintronics-based toolbox that can already be used for the nucleation and read-out of skyrmions [14, 15], thus allowing for the development of technologically consistent devices [71].

A major issue has however impeded the blossoming of stochastic computing as an industrially viable competitor. Namely, cascading gates operating on non-fully random signals propagate unwanted correlations rapidly even after a few elementary operations [63, 72, 73, 74, 75, 52]. In Figure 2b, one can immediately see that feeding two identical copies of the same input signal produces an output signal whose p-value does not represent a product operation in any way. It is therefore crucial to be able to regularly reshuffle signals so that they stay uncorrelated. We call this element a stochastic reshuffler, capable of copying an input stream into an uncorrelated new one while preserving the original p-value. By prefixing such devices to each gate in a circuit (see Fig. 2c), inconvenient correlations can be washed out effectively.

As of present, no efficient hardware stochastic reshufflers are known to exist. In CMOS, reshuffling signals can be done by combining a pseudo random number generator with a shift register or counter (both requiring long term memories) [76]. This means that each reshuffling operation has a large area imprint (a typical integrated shift-register has a linear size of $\sim 10-100\,\mu\mathrm{m}$), and consumes a lot of energy. It is therefore impossible to insert these reshufflers after each calculation stage in a stochastic computing circuit, preventing the realization of any large scale demonstration of stochastic computing on chip. To overcome these technological obstacles, we propose a novel reshuffler design as the first low-energy, compact device proposal of its kind: the Skyrmion Reshuffler [77].

2.2 Skyrmion Reshuffler

Our original concept converts an input bitstream into a sequence of skyrmions whose order can be thermally reshuffled due to the diffusive 2-dimensional dynamical nature of interacting skyrmions at finite temperature. If the order of this sequence can be altered enough by the thermal noise affecting skyrmion motion inside the device geometry, the new skyrmion sequence can then be read as a new uncorrelated output signal with p-value identical to that of the input signal.

The device (depicted in Figure 3) consists of two circular chambers with input/output conduit tracks capable of ushering skyrmions into and out of the chambers. The net drift of skyrmions is achieved by a static current flowing across the entire structure from one conduit to the other (see Methods). Skyrmions injected can be nucleated either singularly via injection of localized spin-currents with magnetic tunnel junctions (MTJ) [3], using in-plane spin-torques on a notch [25], by rapid proliferation with appropriately constructed input conduit geometries [15, 78] or via periodic shedding off of material defects with homogeneous DC currents [79]. Analogously, the read-out elements could consist of MTJs capable of measuring magnetoresistance changes due to the passing of skyrmions through the underlying magnetic layer [80]. We will summarily assume injection/detection capabilities by providing each track with an injection/detection element.

The up-and-down (binary) states of the bitstream are used to select into which chamber to inject the generated skyrmions (see Fig. 3). We do so by selectively nucleating skyrmions at a constant rate ($r_{\mathrm{inj}}=0.2\;\mathrm{nucleations}/n\mathrm{s}$ in the simulations presented here) onto the input conduit of one of the two chambers ($\mathrm{UP}$ and $\mathrm{DOWN}$) depending on the input signal’s state. Doing so generates a total population of $N$ skyrmions proportional to the length $\tau_{\mathrm{input}}$ of the input signal ($N=r_{\mathrm{inj}}\tau_{\mathrm{length}}$) and chamber populations $N_{\mathrm{UP}}=pN$ and $N_{\mathrm{DOWN}}=(1-p)N$ given by the $p$-value of the input signal. The current-induced drift then pushes the skyrmions across their respective chambers toward the output track where they will be detected. The sequence of outgoing skyrmions can then be decodified into an outgoing signal. Whenever an outgoing skyrmion is read-out from the $\mathrm{UP}$($\mathrm{DOWN}$)-chamber we switch the outgoing signal into an up (down) state. As a result of this procedure, the output signal’s length will be given by the time $\tau_{\mathrm{output}}$ taken to read all the skyrmions and, by construction, its p-value will be identical to that of the input signal as long as the relative skyrmion particle numbers in the chambers are conserved. The basic operating principles of the Skyrmion Reshuffler allow it to also work as a device capable of generating stochastic bitstreams with well-defined p-values. To do so, it suffices to inject skyrmions into the two chambers with a specified ratio and allowing the reshuffling to scramble the input order into a random outgoing skyrmion sequence.

In order to quantify the performance of skyrmion reshuffling as a function of chamber materials and geometries, we first model the dynamics of thermally diffusing assemblies of interacting skyrmions.

3 Modelling the Dynamics of Interacting Skyrmions

3.1 Isolated Skyrmion Dynamics

The motion of an isolated skyrmion in a two dimensional film can be described by a modified version of Newton’s equation which tracks the coordinate evolution of the skyrmion’s center. Assuming translational invariance of the skyrmionic spin texture, the Thiele-equation of motion can be used [81, 82, 43, 4]:

[TABLE]

where $\mathbf{\hat{G}}$ is the gyrotropic matrix, $\epsilon_{ij}$ the Levi-Civita tensor, $\delta_{ij}$ a Krönecker delta, and $\sigma$ in (2) is equal to either the damping constant $\alpha$ (when referring to the left side of (1)) or a parameter $\eta$ characterizing both the net ratio of adiabatic and non-adiabatic spin-torque intensities as well as the ratio of damping-like to field-like torques depending on the geometries of current injection. Doing so, our model is general enough to capture equally well the skyrmion driving motion originating from both spin-torques due to in-plane currents as well as spin-orbit torques arising from the underlying heavy metal layer [83, 26, 84]. Lastly, the parameter $\mu$ models the net intensity due to the combination of these two driving effects. Time is in units of $\gamma M$, where $\gamma$ is the gyromagnetic ratio and $M$ the magnetic texture’s local moment magnitude. Under the assumption of an invariant skyrmion profile, the gyrotropic matrix (2) is composed by the gyrovector $G$, a topological invariant arising from the twist in the spin-texture, and the dissipative diadic $D$ (also known as the gyrodamping) which, together with $\alpha$, characterizes the net friction acting on the skyrmion [83, 82, 85, 86, 43]. Both are typically computed explicitly from the static magnetic profile:

[TABLE]

where $\mathbf{n}$ is the local, unit-normalized, magnetization orientation. A very convenient property of both terms is that they do not scale with the skyrmion’s size as long the skyrmion is larger than a few times the domain wall width $\sqrt{A_{ex}/K_{u}}$ (where $A_{ex}$ and $K_{u}$ are the material exchange constant and perpendicular magnetic anisotropy respectively). This ensures that their values computed in zero-temperature micromagnetic simulations will also be valid when modelling finite-temperature skyrmions. Thermal effects are known in fact to excite breathing modes in the skyrmion profile which imply a periodic fluctuation in the skyrmion size [87]. These are expected to qualitatively alter the dynamics of large skyrmions ($\sim 200n\mathrm{m}$) where second order inertial terms should be included in the respective Thiele model [86, 88]. For the physical parameters chosen in our simulations (see Methods), we find stable isolated skyrmions with a typical diameter of $\sim 60n\mathrm{m}$.

The three force terms appearing on the right side of (1) model various interactions that dynamical skyrmions are subject to. Both the conduction of electrons through the spin-texture and the pumping of spin through spin-orbit-torque effects are known to induce both an adiabatic and non-adiabatic spin-torque [89] (usually refered as damping-like / field-like spin-torques when spin orbit torques are considered) giving rise to well-defined force components acting on the skyrmion profile [90, 81]. These current-induced drift effects are captured by the first term where $\mathbf{v}_{s}$ is the spin-drift velocity (directly proportional to current density $j$). In the special case where $\eta=\alpha$, skyrmions drift along the current flow lines. In all other circumstances, however, the discrepancy between $\alpha$ and $\eta$ leads to current-induced skyrmion flows which proceed at an angle to the current-flow driving them [3, 21].

Similarly, the thermal forces are modelled by the collective action of independent random magnetic fields acting on the texture’s local magnetic moments. The net effect results in an additive, homogeneous zero-mean stochastic term to the Thiele dynamics [91, 92]:

[TABLE]

where $k_{B}T$ is the thermal energy and $\langle\cdot\rangle$ represents averaging over noise realizations. From (1), the mean squared displacement of skyrmions can be expected to scale linearly in time with diffusion constant given by:

[TABLE]

The remaining term, which we now proceed to discuss in detail, seeks to model all repulsive interskyrmion interactions and boundary effects:

[TABLE]

where $\mathbf{F}_{S}$ and $\mathbf{F}_{\mathrm{B}}$ are, respectively, skyrmion-skyrmion and boundary repulsion terms.

As described in Appendix A, we will mainly consider magnetic materials with a saturation magnetization of $M_{S}=1400\,k\mathrm{A/m}$, exchange stiffness $A_{ex}=27.5\,p\mathrm{J/m}$, interface-induced DMI constant $D=2.05\,m\mathrm{J}/\mathrm{m}^{2}$, perpendicular magnetic anisotropy constant $K_{u}=1.45\,M\mathrm{J}/\mathrm{m}^{3}$, and $\alpha=0.1$. The choice of these values is consistent with interfacially stabilized skyrmions on Pt/Co/MgO nanostructures [18].

3.2 Skyrmion Interactions

Short range repulsions between skyrmions due to exchange-dominated deformations to their respective spin textures have already been theoretically and numerically discussed in the literature in the absence of magnetic dipole effects [26]. However, dense skyrmion populations are known to exhibit long range order leading to the formation of regular lattice structures [34, 17, 13] for which dipole interactions may play a relevant role. To capture the net sum of these effects on a single skyrmion’s Thiele dynamics, we introduce a force term $F_{S}$ in (8) which sums all the two-body interactions among skyrmions in a given ensemble (in the vein of [26]):

[TABLE]

where $d_{ij}$ is the distance between particles $i,j$ and $\mathbf{\hat{d}}$ the normalized distance vector. The specific exponential form was chosen such that the repulsion behaves gaussian-like at short range and scales like a simple exponential at distances greater than the typical skyrmion diameter.

Whereas the long range exponential scaling is justified in the thin-film limit where stray field energies can be modeled by a local shape anisotropy (see Appendix B), the short range Gaussian scaling was chosen to fit the micromagnetic behavior observed at small inter-skyrmion distances. In fact, for separation distances smaller than the skyrmion diameter, the notion of a skyrmion as a particle is expected to lose its meaning since the complexity of the magnetic texture voids the assumptions giving rise to the Thiele approximation. As such, even though we include a gaussian correction to the skyrmion repulsion, we will not consider scenarios where this correction is relevant. This typically happens whenever the interparticle spacing is smaller than the average skyrmion diameter.

To fit the phenomenological parameters $a_{k}$ (see Methods for values), we performed micromagnetic simulations far from boundaries where the net two-body skyrmion-skyrmion interaction in the presence of dipole effects can be isolated. Two skyrmions were initially set at a very close distance from each other and allowed to relax over time (see Methods for details).

In Figure 4, we plot the gyrotropic matrix elements obtained by using (3-4) from the evolving magnetic profiles. Both $G$ and $D$ are seen to quickly stabilize onto equilibrium values, with $G=-4\pi$ as expected from the skyrmion’s topological charge in the continuum limit [93]. Employing these values, the net repulsion force between the two particles was extracted by computing $\mathbf{\hat{G}}\cdot\mathbf{\dot{x}}$ (the left hand side of (1)) as the simulation progressed. This in turn allows for the exploration of how the net repulsion force scales with the inter-skyrmion distance. Figure 5 shows the particle trajectories (left) along with the extracted repulsion force as a function of inter-particle distance in units of skyrmion diameter (right). We find that the fit (10) is adequate in modeling the observed two-body repulsion resulting from exchange and dipolar interactions. Particularly, we note that the Gaussian corrections to (10) are only relevant at interparticle distances smaller than the average skyrmion diameter.

3.3 Boundary Effects

Skyrmions are experimentally studied on finite geometries. As such, it is essential to account for how the particle dynamics may be affected by the geometry’s boundaries. The radial size of isolated skyrmions in magnetic dots is known to be confined by an explicit condition on the orientation of magnetic moments at the boundary [94, 3]. This effect, a direct result of DMI in the material, is also responsible for repelling skyrmions from boundaries thus guaranteeing the transport properties which make them so useful for applications [3, 21]. As skyrmions in rarefied ensembles (as opposed to lattices) are capable of moving about freely, their dynamics eventually lead them close enough to the sample boundary where their behavior must be quantified through the Thiele formalism.

We repeat our previous phenomenological analysis to model boundary effects by extracting the net force experienced by a solitary skyrmion initially placed adjacent to the boundary of a circular geometry. In Figure 6, the force experienced by the skyrmion results in a net drift both along and away from the dot’s boundary (left) whose scaling behavior as a function of distance from the boundary is distinctly different from that observed for skyrmion-skyrmion repulsions. In order to have an explicit phenomenological expression for the boundary repulsion force, we fit the observed data to a sextic gaussian form (see Methods):

[TABLE]

which is found to model the boundary interaction very well.

Both the skyrmion-skyrmion and skyrmion-boundary interactions should not be considered valid when distances become smaller than the skyrmion diameter. At such length scales the topological texture of each skyrmion is radically altered leading to a behavior both not consistent with the Thiele dynamical model and with the notion of a topologically protected particles. When skyrmions are moved too close to each other or boundaries, annihilation typically ensues which, if anything, should be modelled via an attractive (rather than repulsive) force law [95].

3.4 Ensemble Dynamics

Having modeled both two-particle and boundary repulsions, we now proceed to verify the assumptions of our model. We repeat the same numerical procedures for larger skyrmion numbers to verify that effective forces perceived by single skyrmions can be reconduced to an n-body sum of individual two-particle repulsions plus a boundary interaction term. We simulate 47 skyrmions initially placed randomly inside a circular geometry and allowed them to evolve for $150\,n\mathrm{s}$. In Figure 7 we again see how the gyrotropic matrix elements of each skyrmion are consistent (within a $10\%$ margin of error) with those found for the two-particle simulation in Figure 4. In our particle treatment of the ensemble dynamics, we use identical Thiele parameters for every skyrmion as justified by the small error of our micromagnetically derived gyrotropic constants.

As discussed for the two- and one-particle studies of the previous sections, we reconstruct the net force experienced by each skyrmion from its trajectory and the gyrotropic parameters obtained from the two-skyrmion simulation. In Figure 8, the dynamics of three sample skyrmions are tracked and their net force extracted from the simulated dynamics. We compare this to our phenomenological expressions accounting for just the skyrmion-skyrmion interactions (9) (red curves) and the case where boundary interactions are also included (11) (green curves). The additional boundary term $\mathbf{F}_{\mathrm{B}}$ is included to highlight its importance in modeling skyrmions in the vicinity of the boundary, resulting in a much better fit. Particularly, the boundary effects only become important when skyrmions approach the boundary. As an example, in Figure 8-d, a particle whose trajectory stays close to the center of the nanodot is shown. In this circumstance, the particle’s force fit with and without the boundary correction term are found to be virtually identical. The excellent agreement between micromagnetic simulations and our fit demonstrates that our multi-body Thiele model accurately captures dipolar and exchange skyrmion-skyrmion interactions as well as skyrmion-boundary interactions. We now have all the tools necessary to efficiently model large ensembles of skyrmions in a device setting.

4 Skyrmion Reshuffler

Having evaluated the behavior of single skyrmions in a confined system and extracted the forces and topological parameters of interacting skyrmions as well as the influence of sample boudaries, the next step towards the implementation of bio-inspired skyrmion-based devices is to apply the numerical techniques developed to simulate an actual real-world device. We begin with the Skyrmion Reshuffler introduced in Section 2 and then propose a completely different application of the same device as a neuron in the next section. Many other future devices employing the same principles as ours should be possible, as well as being numerically tractable via the modelling techniques just introduced.

We simulate the Skyrmion Reshuffler in Figure 3 by considering two identical circular $1024\,n\mathrm{m}$ chambers into which we inject skyrmions every $5\,n\mathrm{s}$ depending on the state of a $0.5\,m\mathrm{s}$-long telegraph noise signal with fixed p-value. The skyrmions are initially injected from $100\,n\mathrm{m}$-wide input conduits with currents of varying intensity (see Methods). After drifting and diffusing through their respective chambers, they are then read upon arriving at the output track and used to reconstruct an outgoing signal. The output signal’s correlation to the input is then checked by computing the product-moment correlation coefficient of the two [96]:

[TABLE]

where $\mathrm{cov}(X,Y)$ is the covariance of two signals $X$ and $Y$ while $\sigma$ is their respective variance. We show this result in Figure 9 as a function of the drift current intensity for different signal $p$-values.

If the chamber is significantly larger than the transverse size of the input/output conduit tracks, current densities inside the chamber may become small enough, allowing skyrmions to interact and diffuse thermally in the chambers before exiting. In Figure 9-left, the correlation coefficient $\rho$ is in fact seen to decrease rapidly with lower current intensities. For the physical values considered, our device shows that strong decorrelation is achieved for currents intensities of $\sim 10^{10}\;\mathrm{A}/\mathrm{m}^{2}$ 222For the lowest currents sampled no more than $51$ particles were observed in either chamber at any given time. As we will see later on (see Figure 13), this number is well below the saturation limit of $90$ for this chamber size.. From the moment the output signal reconstruction commences, the correlation reaches its steady state value on a $\mu\mathrm{s}$ timescale as can be seen in Fig. 9-right by the flattening of the correlations’ temporal traces. The time required to inject and read-out a statistically significant number of skyrmions determines the precise timescale for relaxation of correlations towards their steady-state equilibrium. This leads to the intuitive understanding that stronger currents allow the device to function faster at the cost of decorrelation efficiency. In turn, what exactly qualifies as ’satisfactory decorrelation’ will ultimately depend on the details and specifications of the stochastic computing circuit one wishes to implement. Furthermore, we don’t see significant deviation in the results as a function of the p-value assigned to the input signal. In Fig. 10 we instead fix the current intensity at $1.7\cdot 10^{11}\;A/\mathrm{m}^{2}$ and allow for the chamber sizes to vary. Larger chambers allow the particles to diffuse for longer thus resulting in a better decorrelation of the shuffler’s output signal. Ultimately, Figures 9 and 10 show how the degree of decorrelation achieved can be tuned by both lowering the current intensity and/or increasing the chamber sizes.

Ideally, one would like the signal scrambler to work as fast as possible. The speed at which a statistically significant output signal can be generated depends directly on the intensity of the driving current moving the injected skyrmions through the chambers (we neglect effects due to finite response times of nucleation and read-out mechanisms). However, as shown in Figure 9, higher currents lead to highly correlated output signals due to the fact that skyrmion particles do not spend enough time diffusing and interacting in the chambers. In Fig. 9-left we find $64\%$ correlation at the highest driving current sampled ($3.8\cdot 10^{11}\;\mathrm{A}/\mathrm{m}^{2}$).

Alternatively, one may try to increase the temperature in an attempt to amplify thermal diffusion. This, however, leads to both a higher energetic consumption as well as the potential for an increased skyrmion annihilation rate. Alternatively, one could attempt to leverage the presence of random material inhomogeneities to further drive the diffusive processes. This has been recently discussed and characterized in the literature [97] as a possible technique to tune the skyrmion dynamical properties. Ultimately, the condition affecting this trade-off lies in the relative intensity of dynamical drift and thermal noise. From (7) one has that the skyrmion’s diffusivity over time results in a dispersion given by:

[TABLE]

where $t$ is the elapsed time and angle brackets $\langle\cdot\rangle$ refer to averaging over noise.

Significant decorrelation between input and output signals takes place whenever this diffusion is large enough to allow scrambling in the sequence of skyrmions being injected into the chamber. Denote then by $l_{0}$ the relative spacing between successively injected skyrmions. If thermal effects manage to displace skyrmions more than their inter-spacing ($\langle|\mathbf{x}|\rangle\gg l_{0}$) over the time it takes them to cross the chamber, it is reasonable to expect the ingoing and outgoing sequences to become randomized relative to each other. As such, a qualitative estimate of the correlation between input and outgoing signals can be found in the ratio between these two length scales $\rho\propto 1/(1+\sqrt{\langle|\mathbf{x}|^{2}\rangle}/l_{0})$. Recalling then from (1) that the particle drift speed scales with the current intensity ($|\dot{\mathbf{x}}|\simeq\mathbf{v}_{s}\propto J$) one can expect that correlations will scale relative to system parameters as:

[TABLE]

where, denoting the chamber diameter as $L$ and the injection frequency as $\nu$, the time taken to cross the chamber is $t_{\mathrm{exit}}\simeq L/|\mathbf{v_{s}}|\propto L/J$ and the inter particle spacing is $l_{0}\propto\nu\,J$. The fitting curves in Figures 9 and 10 indeed verify that the scalings $\rho/(1-\rho)\propto J^{3/2}/\sqrt{L}$ match our numerical data, thus justifying how correlations will be lower for smaller currents and larger chamber sizes.

A potential improvement to the Skyrmion Reshuffler can theoretically be found in the use of anti-ferromagnetic (AFM) materials where the possibility of stabilizing AFM skyrmions has recently been suggested [98, 99, 100]. At the core, the main differences in the material parameters between AFM and FM devices lies in the sign of the exchange interaction and the magnitude of the magnetostatic interaction. These differences have potentially large repercussions on the Thiele modeling of skyrmions in the AFM systems. The negative exchange interaction leads to a vanishing gyrovector $G=0$, implying two important consequences. First, the AFM skyrmion always has zero transverse velocity relative to the direction of current flow. This, in FM skyrmions, is true exclusively in the singular scenario $\alpha=\eta$. Secondly, and most importantly, a vanishing gyrovector implies a larger diffusion constant (7) and, consequently, smaller correlations (14). Depending on the value of the damping constant $\alpha$ in the materials employed ($\alpha\sim 0.01-1$), the difference in magnitude of the diffusion constant in FM and AFM skyrmions can easily span an order of magnitude difference. The magnitude of the magnetostatic interaction also complements the AFM skyrmion as it may lead to the stabilization of much smaller, and more thermally sensitive, particles. For all these reasons, one expects that AFM skyrmions would allow operation of our proposed device at higher currents and faster speeds.

Any accidental particle nucleation or annihilation taking place in the chambers will not influence the outgoing signal’s p-value as long as its rate remains proportional to the total number of skyrmions present in each chamber. This proportionality is guaranteed for accidental particle annihilations as it will depend on how many skyrmions are present in the chamber at any given time. Nucleations on the other hand can be potentially problematic for the device’s proper functioning. Luckily, the experimental observations of skyrmions at room temperature show no signs of thermally influenced skyrmion instabilities. This leaves the interactions with material defects and background inhomogeneities [101] as the only other potential source of skyrmion anihilations. As long as both chambers are similarly constructed, however, this will lead to a similar rate of annihilations in both chambers thus preserving input-to-output p-values modulo a possible decrease in the output signal’s frequency.

Overall, given the DC-current driving the device (and possibly nucleating the necessary skyrmions as in [79]), we expect the skyrmion reshuffler to operate at very low energies ($\sim\mu\mathrm{W}$). A more accurate metric would involve calculating the energy cost per reshuffled bit of information being transmitted through the telegraph signal. Typical CMOS-based solutions reshuffle clocked stochastic bitstreams [76] whose energy consumption can be estimated at $\sim 10\;p\mathrm{J}/\mathrm{bit}$. Placing their large area imprint problems aside, they are expected to perform efficiently only for slow bitstreams as their energy consumption will scale linearly with the signal’s frequency. On the other hand, our Skyrmion Reshuffler constant energy costs will prove ideal for reshuffling high frequency streams instead. As already mentioned, the device modeled in this section considered a constant skyrmion injection rate of $0.2\;\mathrm{nucleations}/n\mathrm{s}$ thus allowing us to potentially reshuffle signals with frequencies as high as $5\;G\mathrm{Hz}$, resulting in an equivalent energy-per-bit cost of $5\;f\mathrm{J}/\mathrm{bit}$ which is orders of magnitude lower.

5 Skyrmion Neuron

We now turn to demonstrate how a similar design can be used to emulate the behavior of a neuron. Neural functioning relies on the ability to transmit and receive electrochemical signals to/from a large number of interconnected neurons by means of its dendritic connections (see Fig. 11). As coupled neurons spike, the receiving neuron registers the event with a voltage increase (or decrease) across its cellular membrane. The underlying physical mechanism for this behavior is the presence of ion channels scattered over the cell’s membrane responsible for these potential differences. Such physical channels are not ideal however and leak over time resulting in a steady decay of the neuron’s voltage towards its resting potential in the absence of incoming electrochemical pulses. As a result of this behavior, a neuron can be thought of as a leaky capacitor.

If incoming neuronal activity is registered at a high enough frequency, however, the receiving neuron overcomes its leaking voltage and steadily charge. A threshold potential exists beyond which the neuron suddenly discharges and resets. This basic model of neuronal functioning is known as the leaky integrate-and-fire neuron model [102, 103]. The neuron effectively serves as a basic memory integrating the sum of past voltage spikes received with the added property that this memory is lost over time if a potential threshold is not reached.

The analogies to our skyrmion reshuffler are already apparent as the injection of particles into it can be thought of as a response to a voltage spike train (shown in Fig. 3). In this sense, the collecting of skyrmions inside the chambers serves as a history of the input spike train. If the passage of skyrmions through the output conduit were to be precluded (such as with a voltage gate [35, 104]), skyrmions progressively accumulate inside the geometry where their density can be estimated via a reading element placed inside each chamber (see Fig. 12). Each chamber can then be thought of as a reservoir collecting the memory of whatever input signal was employed to populate it.

Furthermore, an upper limit exists to how many skyrmions may be crammed into each such chamber. In fact, as skyrmions become compressed with their neighbors to length scales comparable to the typical skyrmion diameter, the combination of dipolar strain among particles along with thermal effects breaking topological stability may cause them to annihilate. The topological stability of the skyrmion profile is only guaranteed under continuous deformations of the magnetic texture [24, 105]. The introduction of a sufficiently large thermal field in the magnetization dynamics can make a skyrmion collapse with some probability [14, 106, 107, 108, 109, 110]. Alternatively, and possibly more importantly, skyrmion annihilation can be biased by the presence of randomly scattered material defects in the magnetic sample populated by the skyrmions. As such, small room temperature skyrmions can potentially annihilate (and nucleate) simply due to thermal fluctuations and material disorders. The maximal particle density inside a chamber can hence be thought of as a saturation threshold beyond which memory of the input signal will certainly be destroyed.

In Fig. 13, we perform micromagnetic simulations at $300\,K$ to establish the maximal number of $\sim 39\,n\mathrm{m}$ skyrmions which can be stabilized inside chambers of varying sizes. We nucleate the skyrmions randomly in the chamber in ever increasing numbers and allow the magnetic texture to stabilize for $2\,n\mathrm{s}$ before counting how many skyrmions are left. This process was then repeated ten times before computing the average. For each chamber diameter explored, a saturation value is found beyond which no more skyrmions can be added reliably into the geometry.

The combination of memory due to particle conservation and memory-loss from annihilations allows one to view our proposed skyrmionic chamber as a lossy memory storage device which integrates an incoming signal while allowing a certain decay of the accumulated information. Upon switching off the voltage gate, the collected skyrmions can be released into the output conduit whenever a certain threshold skyrmion density is reached (as measured by a reading element inside the chamber). Following these principles, the functioning of our device behaves as a leaky integrate-and-fire neuron where, instead of integrating discrete current pulses through electric potential accumulation before discharging, our device collects skyrmions instead and successively dumps them into the output conduit.

6 Conclusion

The impracticality of performing complex micromagnetic simulations out to long timescales currently poses an obstacle to the design and proof-of-concept of novel spintronic-based devices with a $\sim\mu\mathrm{m}$ area imprint. Under suitable simplifying assumptions however, the behavior of solitonic magnetic spin-textures can be radically simplified and reduced to classical n-body problems which can be tackled through a plethora of numerical techniques. We have shown how one such particle model (the Thiele equation) can be phenomenologically tuned in such a way to account for complex exchange and dipole mediated interactions among stable magnetic skyrmions. By means of the tools developed, we proposed and modeled two skyrmion devices whose behavior would be computationally intractable via standard micromagnetic techniques.

The Skyrmion Reshuffler encodes a telegraph noise signal into a sequence of skyrmions injected and current-driven through two separate chambers. The dominance of thermal diffusion in each chamber scrambles the particle order in such a way that the output sequence can be used to reconstruct an output signal with negligible correlations to the input one. Through simple scaling analyses, we have shown how the system parameters impact these correlations and balance the trade-offs between speed and effectiveness of the device’s operation. Leveraging the interplay between skyrmion particle stability and thermal diffusivity, we have argued that the skyrmion reshuffler stands to solve a long standing problem in the field of stochastic computing where computations defined over stochastic signals require correlations to be absent for proper functioning. The spintronic basis of the skyrmion reshuffler’s design allows for a highly attractive area imprint and energy costs.

The Skyrmion Neuron adapts the skyrmion reshuffler’s design to propose a device capable of emulating the integrate-and-fire characteristics of a neuron. By means of an added voltage gate selectively prohibiting the outflow of skyrmions from the chamber, as well as an extra read-out element inside the chamber, the skyrmion density in the chamber can be measured and the voltage gate released upon a certain threshold being reached. In doing so, the chamber integrates the signal intensity (or length) driving the input sequence of skyrmions injected into the chamber. To complete the analogy with the animal neuron, we argue that skyrmions perform the role of both voltage accumulation as well as neurotransmitters in this specific layout. The devices we propose are compatible with standard CMOS and can also be the building block of future fully spintronic probabilistic computers.

Overall, the sensitivity of magnetic soliton structures to small driving currents and thermal noise makes them optimal candidates for the development of novel probabilistic devices. Furthermore, the nano-scale size of these particles implicitly guarantees that devices employing them will be both energy efficient and scalable. Our results suggest that the basic research currently being done on these exotic magnetic solitons fertilizes the field for future disruptive applications to come.

Acknowledgments

This work was funded by the European FET OPEN ‘Bambi’ project 618024 and the European Research Council ERC grant ‘bioSPINspired’ 682955. J.V.K and V.C acknowledge FET OPEN ‘MAGicSky’ project 665095. The authors would like to thank A. Thiaville for his very useful comments during the revision of the manuscript.

Appendix A Methods

The three-dimensional micromagnetic simulations were performed using the MuMax3 GPU-accelerated micromagnetic simulation program [111]. The average energy density contains the exchange energy, the anisotropy energy, the applied field (Zeeman) energy, the magnetostatic (demagnetization) energy and the DMI energy terms. In all simulations, the thickness of the magnetic nanotracks was $2\,n\mathrm{m}$. The length of the input/output conduits was $50\,n\mathrm{m}$, while the width was set to $100\,n\mathrm{m}$. Magnetic parameters used in the simulations: saturation magnetization $M_{S}=1200\,k\mathrm{A/m}$, exchange stiffness $A_{ex}=8\,p\mathrm{J/m}$, interface-induced DMI constant $D=2.1\,m\mathrm{J}/\mathrm{m}^{2}$, perpendicular magnetic anisotropy constant $K_{u}=0.9\,M\mathrm{J}/\mathrm{m}^{3}$ and gyromagnetic ratio $\gamma=-2.211\cdot 10^{5}\,\mathrm{m}/\mathrm{A}\mathrm{s}$. The Gilbert damping coefficient $\alpha$ was set to $0.1$. The choice of these values is consistent with interfacially stabilized skyrmions on Pt/Co/MgO nanostructures [18].

All models are discretized into tetragonal cells with the constant cell size of $2\times 2\times 2\,n\mathrm{m}^{3}$ in the simulations, whose linear size is smaller than both the fundamental length scale $A_{ex}/D\simeq 4\,n\mathrm{m}$ and the domain wall width $\sqrt{A_{ex}/K_{u}}\simeq 3\;n\mathrm{m}$. With these parameters, we observed the stabilization of skyrmions with an average of diameter of $\sim 39\,n\mathrm{m}$ both at zero and $300\,K$ temperature. The integration time steps used were $5\cdot 10^{-12}\,\mathrm{s}$ and $5\cdot 10^{-14}\,\mathrm{s}$ at $0\,K$ and $300\,K$ respectively. At finite temperature breathing modes were observed at GHz frequencies.

The tracking of skyrmion particles was performed on dynamical snapshots taken at each time step using the python implementation of the DLIB open-source machine learning library [112]. The classic Histogram of Oriented Gradients (HOG) feature combined with the linear classifier, image pyramid, and sliding window detection scheme were used to train a detector to recognize skyrmions using a training set of $200$ stabilized skyrmion profiles. From the geometry’s total magnetization profile, the tracking scheme allowed us to obtain a set of regions containing each one skyrmion only.

The computation of the skyrmion center was performed by first fitting the magnetic profile of each region with a cubic spline interpolation and then finding the maximum z-component of the magnetization on the interpolated profile. The algorithm allows for a $\sim p\mathrm{m}$ resolution of the skyrmion center. The gyrotropic parameters were computed by performing the relelvant surface integrals on the interpolated profile of each skyrmion particle as discussed in equations (3) and (4).

The n-body Thiele dynamics were solved using a specialized GPU-accelerated CUDA/C++ solver developed by the authors. The stochastic dynamical system of equations was solved employing a Heun scheme ensuring convergence to the proper Stratonovich solution. The current density distributions were first obtained by simulating current flow on COMSOL [113] for the geometry discussed. Room temperature skyrmions are known to exhibit breathing modes [87], these will however not impact the net dynamics due to the gyrotropic terms being effectively scale-free in the particle size [86] when the skyrmion size is much larger than the domain wall width $\sqrt{A_{ex}/K_{u}}\simeq 4.35\;n\mathrm{m}$. Qualitatively different dynamical behavior can however be expected for large skyrmions where an extra mass term is expected to play a role in the Thiele-description. For the material parameters in question, the numerical fitting constants for the skyrmion-skyrmion two-body repulsions were found to be $\mathbf{a}\simeq[2.709,-5.643,0.964]$, whereas the fitting constants for boundary repulsions were $\mathbf{b}\simeq[0.001419,-0.02631,0.24135,-1.1609,3.3547,-2.1786]$. In this work, we have taken $\eta=\alpha$ and $\mu=1$ for convenience.

Appendix B Qualitative influence of dipole effects on skyrmion repulsion in the ultrathin film limit.

A magnetic texture $\mathbf{M}=\mathbf{M}(x,y,z)$ with constant magnetic saturation $|\mathbf{M}|=M_{S}$ extended over a 3-dimensional domain $\Omega$ of thickness $d$ can be described by a micromagnetic free energy which, in the presence of an out-of-plane uniaxial anisotropy (along the $z$-axis) and an interfacial Dzyaloshinskii-Moriya interaction (DMI), takes the form [114, 115]:

[TABLE]

where $A$, $K$ and $D$ are the exchange stiffness, anisotropy and Dzyaloshinskii-Moriya strength respectively. In order of appearance, the terms appearing in the magnetic free energy are: the exchange term, which prefers constant magnetization configurations; the magnetocrystalline anisotropy, which favors out-of-plane magnetization configurations; the magnetostatic term, which prefers divergence-free configurations; and the surface DMI term, which favors chiral symmetry breaking.

In the limit of ultrathin ferromagnetic films where the sample thickness is smaller than the exchange length $l_{ex}$ ($d\leq l_{ex}=\sqrt{2A/(\mu_{0}M_{S}^{2})}$) the magnetic free energy above can be described by a model where the stray field energy is expressed as a local shape anisotropy term [116, 117]. This model can in turn be used to treat the simplified 1-dimensional case of magnetic wires in the presence of DMI with magnetic energy given by [118]:

[TABLE]

where $\kappa=D\sqrt{\frac{2}{\mu_{0}M_{S}^{2}A}}$, $Q=\sqrt{\frac{2K}{\mu_{0}M_{S}^{2}}-1}$, and the magnetic profile has been parametrized as $\mathbf{M}=M_{S}\,(\sin\theta,0,\cos\theta)$. The dimensionless constant $Q$, known as the quality factor, captures both the contributions of the crystalline anisotropy and the stray field in the ultrathin film limit considered. The energy is minimized by domain wall solutions whose explicit form can be verified to be:

[TABLE]

Let us then consider a magnetic texture consisting of two such domain walls separated by a distance $\delta$:

[TABLE]

with $\eta=\exp(-Q\delta/2)$. The total energy of the system can be written as:

[TABLE]

[TABLE]

where the first two terms are simply the minimum energy of each domain wall taken independently while the third measures the interaction potential of the two walls. The repulsive force between the two domain walls will then be straightforwardly given by the derivative of interaction potential with respect to the distance separating the walls ($F_{\theta_{+},\theta_{-}}(\delta)=-\partial_{\delta}E^{\mathrm{int}}_{\theta_{+},\theta_{-}}$). Tedious, but fairly elemental analysis then results in the following repulsion force:

[TABLE]

The repulsion force is effectively dominated by the behavior of the powers of $\eta$ appearing in the prefactors. As such, the repulsive force can be expected to decay exponentially in the distance $\delta$ between the two domain walls. This result is a direct consequence of the ultrathin film limit which allows one to write the magnetic stray field as an effective shape anisotropy and is not expected to qualitatively change for the case of 2-dimensional magnetic skyrmions considered in the main body of the paper.

Appendix C Skyrmion properties under different material conditions.

In this section we list all the material parameters explored throughout our micromagnetic simulations. In order to perform our analysis of the skyrmion-skyrmion repulsive forces in the main text, we required that the micromagnetically simulated magnetic textures relax to their stable values within short enough timescales such that all further dynamical behavior could be considered physically relevant. The seven sets of material values shown below are named S1-S7 (S2 corresponds to results in the main text).