Introduction

Magnetic skyrmions are whirls of magnetization stabilized by their nontrivial topology1,2,3,4,5,6,7. They have been observed in bulk materials and magnetic thin films at temperatures ranging from a few Kelvin to far above room temperature8,9,10,11. Skyrmions can be stabilized even without the application of an external field12,13, can be manipulated efficiently with spin torques9,14,15 and exhibit thermally activated dynamics16,17,18. This turns them into attractive candidates for applications in novel devices relying on skyrmion-based logic18,19, Brownian computing20,21, and racetrack memory7,22,23,24,25. High-density lattice states16,26,27,28,29, in which skyrmions strongly interact with each other, have recently moved into the focus of attention. In contrast to colloids30,31,32, sizes and density of skyrmions can be adjusted on the fly which makes micrometer-sized skyrmions in ultrathin films act as ideal model systems to study phase behavior in two dimensions33 as shown by the demonstration of hexatic and solid phases16,28.

To theoretically describe skyrmion ensembles that for the experimentally relevant systems span length scales over microns to hundreds of microns, conventional micromagnetic approaches are not applicable due to the prohibitive computational cost. Therefore coarse-grained particle-based descriptions have been introduced34,35,36,37,38,39,40,41,42, where skyrmions are represented by repulsive soft disks whose dynamics are governed by the Thiele equation43. While this approach is successful in describing individual skyrmion dynamics28,44, modeling of skyrmion ensembles or of skyrmion movement in tight confinements21,25,45 requires knowledge of the interaction potentials. So far, interaction potentials between skyrmions are mostly rationalized by theoretical considerations. If one considers skyrmions as quasiparticles, then the description in the simplest approximation requires interaction potentials between a skyrmion and a sample edge as well as between a skyrmion and another skyrmion. While possible other interactions can be envisaged, one would start with these two simple descriptions and then compare to experiments to see if these suffice to describe the dynamics of multiple skyrmions in confined geometries. So far, the skyrmion interaction potentials have been analyzed based on a micromagnetic energy functional. By determining the spin structure of a skyrmion, interaction potentials between two such structures and between a skyrmion and the magnetic material boundary have been approximated by a potential which is exponentially decreasing with distance34,38,39,46,47,48. Concrete parameters for coarse-grained particle-based potentials at the nanometer scale are obtained by fitting micromagnetic simulation data to this functional form34,37. Reference 28 followed a different approach by directly matching structural information from experimental high-density states with corresponding particle-based simulations. Instead of an exponential function, a repulsive power law was used to enable comparisons with previous simulations of phase transitions33. This approach, however, entails the severe limitation that the interaction potentials are purely repulsive and therefore unable to capture possible attractive parts in the interaction as have recently been theoretically predicted49,50,51 and experimentally shown for 3D skyrmion tubes in the conical phase52,53 as well as skyrmions in frustrated magnets54,55. Methods such as Reverse Monte Carlo56 or Iterative Boltzmann Inversion (IBI)57,58, which are rooted in computational soft matter physics, do not require such assumptions and are inherently better suited for the analysis of the full interaction potentials. In these approaches, potentials are adjusted in successive simulations to match a target pair-correlation function, and employed to construct coarse-grained computational models from more detailed atomistic descriptions. Even though the possibility of experimental reference functions is often discussed in this context59, they have so far been rarely employed.

In this work, we use Iterative Boltzmann Inversion to derive skyrmion–skyrmion and skyrmion-boundary potentials for particle-based simulation models directly from experimental data for skyrmions where the chiral spin structure is stabilized by the Dzyaloshinskii–Moriya interaction and without any assumptions on the shape of the potential. We demonstrate that based on these interactions we can reproduce static properties of the experimental system with particle-based computer simulations. In addition, we find that the two interactions can be described by purely repulsive exponential potentials.

Results and discussion

Description of the system

We propose a method with which we can derive system-specific coarse-grained potentials for skyrmion–skyrmion and skyrmion-boundary interactions from the analysis of a single measurement without applying assumptions about the shape of the interaction potentials using the Iterative Boltzmann Inversion method. This method is exemplarily applied to a system of μm-sized skyrmions in a multilayer stack with a confining stripe geometry as depicted in Fig. 1. More details on the sample and experimental conditions are found in “Methods”. By determining potentials this way, we find that the structural properties of a skyrmion lattice with μm-sized skyrmions are consistent with soft repulsive potentials. The IBI method has two important advantages over similar methods used in previous publications28. First, no assumptions on the shape of the potential must be provided for the simulations. Instead, the shape is obtained directly from the experimental data. Second, IBI allows for determining skyrmion–skyrmion and skyrmion-boundary potentials simultaneously by updating the two in every step. This method can be easily extended to include interactions with a nonflat inhomogeneous energy landscape60. In addition, spin–orbit/transfer torques due to effective spin-polarized currents can be modeled in a straightforward manner by applying forces to all quasiparticles in the system.

Fig. 1: Magneto-optical Kerr effect microscopy image of an experimental skyrmion conformation.
figure 1

The field of view covers part of a long rectangular stripe of magnetized material along the horizontal direction with a width of 81 μm. The two parallel yellow dashed lines denote the location of the boundaries of the magnetized area. Skyrmions are displayed as bright bubbles on gray background and marked by red circles as detected using the Trackpy module62,68 for python. The diameter of the circles is about 5 μm and smaller than the actual interaction range. The region enclosed by the cyan dashed line denotes the area from which skyrmions are taken as reference particles in the determination of the pair-correlation function g(r). In the region enclosed by the green dashed line, skyrmions are still counted as neighbors for skyrmions from the inner part for the calculation of g(r). The distance between cyan and green dashed lines is 20 μm to keep g(r) properly normalized in the range of interest.

Full size image

We employ particle-based Brownian Dynamics simulations34,35,36 in the Thiele model43 with potentials determined using the IBI method to match the radial distribution functions of our simulation to the experimental function as depicted in Fig. 2. More details on these methods can be found in “Methods”. The method reproduces both distribution functions very closely within the fitted region and also to a good extent beyond this region at least for the skyrmion–skyrmion interaction. The different shape of the skyrmion-boundary interaction beyond the fitted region is most likely a consequence of pinning-induced quenched disorder which is not considered in the simulation.

Fig. 2: Plots of experimental and simulated radial distribution functions.
figure 2

The experimental (red) and simulated (blue) distribution functions are given for skyrmion–skyrmion (a) and skyrmion-boundary interactions (b). The simulated distribution function is generated by a Brownian Dynamics simulation after performing 500 update steps as described in the main text. (b) additionally contains a corrected experimental function which is used to determine the simulation potential (black dashed line). The bin width Δr of of all functions is 0.1 μm.

Full size image

Shape of interaction potentials

Potentials resulting from this method are shown in Fig. 3. The two potentials we determined for this experimental system are purely repulsive and can be approximated by exponential functions. While the agreement with the generic fully repulsive functional form for interactions of chiral skyrmions predicted from theoretical considerations and micromagnetic simulations for skyrmions on the nanometer scale34,39,47 is good for our extracted skyrmion–skyrmion interactions, somewhat larger deviations can be observed for skyrmion-boundary interactions. It is worth noting that while the potentials presented here turn out to be fully repulsive, IBI would in principle be able to capture attractive interactions57. Particle-based simulations may even be adapted on-the-fly to scenarios in which skyrmion potentials change from attractive to repulsive as described for skyrmions occurring in the conical phase by Du et al.53 or in which the skyrmion sizes change as a result of applied fields as long as radial distribution functions for the extreme positions are accessible.

Fig. 3: Interaction potentials determined by Iterative Boltzmann Inversion.
figure 3

Skyrmion–skyrmion (blue) and skyrmion-boundary (red) interactions are shown. Fits of the potentials to an exponential function \(V(r)=a\exp (-r/b)\) with distance r and fit parameters a and b (dashed lines) are given in the region from 4.0 μm up to the first minimum of the respective distribution function. Fit values for skyrmion–skyrmion interaction are aSkSk = 735.1 kBT and bSkSk = 1.079 μm. Values for skyrmion-boundary interaction are aSkBnd = 176.7 kBT and bSkBnd = 1.673 μm. The parameter a is given in units of the Boltzmann constant kB times temperature T.

Full size image

Verification of results

Finally, we test the capability of our potential to reproduce independent static properties of the experimental system. To this end, we turn to structural quantifiers which describe phases and ordering of 2D quasiparticles. We start with the local orientational order parameter33 \({\psi }_{6}(j)=1/{n}_{j}\mathop{\sum }\nolimits_{k = 1}^{{n}_{j}}{e}^{6i{\theta }_{jk}}\), where nj represents the nearest neighbors of skyrmion j as determined by Voronoi tessellation, and θjk is the angle between the line connecting skyrmions j and k and the x axis33. Considering skyrmions in the cyan reference region defined in Fig. 1, the absolute value of the local order parameter 〈∣ψ6∣〉 as averaged over all skyrmions and snapshots is 0.486 ± 0.008 for the experimental system and 0.496 ± 0.007 for the simulated system. Similarly, favorable agreement can be obtained by the percentage of particles with six neighbors P6. For the experiment, we obtain an average P6 = (59.2 ± 0.8)% and P6 = (62.0 ± 1.1)% for the simulation. The ordering of skyrmions beyond the first shell can be visualized by the normalized two-dimensional radial distribution function g(x, y) shown in Fig. 4 which is calculated with the same constraints as in 1D. The six peaks around the center show the preferred local sixfold coordination of the system while the nearly isotropic second circle implies the absence of long-range orientational order. We conclude that the system is in the liquid phase28,29 as already indicated in Fig. 1. Excellent agreement between experiment and simulation with respect to structural quantifiers is found (considering that pinning effects have been neglected in the simulations) confirming the capability and accuracy of our approach. It also shows that pinning effects are reduced in skyrmion lattices due to the skyrmion–skyrmion interaction. While averages over the whole system indeed show results very similar to the unpinned case, the local structure at specific positions such as at a given distance to the boundary is still somewhat impacted by pinning.

Fig. 4: Normalized two-dimensional radial distribution functions in experiment and simulation.
figure 4

a The two-dimensional radial distribution function g(x, y) of the experimental system averaged over all frames of the experimental video. b g(x, y) of the simulated system. The greyscale bar depicts the normalization of skyrmion occurrence probabilities in the two plots.

Full size image

Conclusions

In summary, we have described a method for determining coarse-grained skyrmion interaction potentials, thus bridging the gap in time and length scales between micrometer-sized experimental skyrmions and theoretical predictions so far possible only at the nanoscale34,37. Thereby, our work captures the physics on experimental time and length scales where fully atomistic or micromagnetic simulations are typically not feasible, but are relevant for devices based on skyrmions. Our results allow for more accurate simulations of micrometer-scale skyrmion systems without being compelled to use potentials determined for inherently different systems. Our procedure is well-suited to construct coarse-grained simulation models for concrete experimental setups and is a first step towards the quantitative exploration of novel skyrmion devices with computer simulations. The method can also be applied to a broad range of experimental systems. These include not only the multilayer stacks leading to skyrmion lattices as studied here but also much smaller skyrmions that can be studied in bulk systems hosting skyrmions16 making this a universally applicable method as long as pair-correlation functions are accessible. The latter can also be accessed for nanoscale skyrmions, which cannot be easily imaged by microscopy. They can however be characterized by neutron diffraction26,27. During such experiments, the scattering pattern can be easily converted to the pair distribution function by a Fourier transform. Thus, the approach of extracting skyrmion–skyrmion interaction potentials from a pair distribution function has the potential to be applied to both, nanoscale16 and microscale skyrmions28. However, this applicability to different kinds of systems does not guarantee similar results for the shape of the interaction function. As a consequence, one can use our method to quantify interaction models between skyrmions including those with partially attractive interactions found for skyrmions in the conical phase and skyrmions in frustrated magnets.

A comprehensive and quantitative description of experiments will also require a description of the effects of nonflat energy landscapes60, which could be integrated into our approach as well as dynamical properties. Moreover, analyzing skyrmion systems at different densities and comparing results between these systems could provide further insights in the properties of skyrmion interactions that might go beyond the two-body description explored in this work.

Methods

Experimental methods

Skyrmion interactions were analyzed in a Ta(5)/Co20Fe60B20(0.9)/Ta(0.08)/MgO(2)/Ta(5) (layer thickness in nm in parentheses) multilayer stack which had already been used in the observation of skyrmion lattices in confined geometries17. The sample exhibits ferromagnetic properties with interfacial Dzyaloshinskii–Moriya interaction due to the Ta(5)/Co20Fe60B20(0.9) interface for which we determine a value of D = (0.30 ± 0.10) mJ/m2 18,61. The perpendicular magnetic anisotropy arises due to the MgO and is fine-tuned by the Ta dusting layer. It is determined as Keff = (62 ± 11) A/m using a superconducting quantum interference device with a saturation magnetization of MS = (0.71 ± 0.03) MA/m and a Curie temperature of TC = 448 K28. The exchange constant for the used material is reported to be A = 10 pJ/m. The out-of-plane component of magnetization was imaged using polar magneto-optical Kerr effect microscopy. Rectangles of magnetic material were patterned on the sample with electron beam lithography and Argon ion etching (Fig. 1). The width of the rectangle was 81 μm as estimated from the background image. On the two short sides, the boundaries were far enough outside of the analyzed region to have negligible edge effects. The temperature of the sample (341.5 K) was controlled by a Peltier element and tuned by the current source. This temperature allows for creating a number of well-visible and detectable skyrmions with a density that can be controlled by field cycling. Skyrmions were nucleated after an in-plane magnetic field pulse of 30 mT which saturated the magnetization in the plane under a fixed out-of-plane field of 70 μT. Particle tracking was performed with the python module TrackPy62 (see Supplementary Note 1 and Supplementary Fig. 1 for more information on the tracking procedure), and a trajectory video was taken for 20 min at a frame rate of 16 frames per second. As skyrmion density still decreased in the first 5 min and stabilizes afterward, only the remaining 15 min of the measurement were used for analysis.

Analysis of distribution functions

From these data, we compute the average of the one-dimensional radial distribution function for each frame (Fig. 2a):

$$g(r)=\frac{1}{\rho {N}}\frac{1}{2\pi r{{\Delta }}r}\mathop{\sum}\limits_{i\ne j}\left[\theta ({r}_{ij}-r)-\theta ({r}_{ij}-r-{{\Delta }}r)\right].$$
(1)

Here, ρ is the skyrmion area density calculated from the region inside the green dashed lines in Fig. 1. N is the number of skyrmions enclosed by the cyan lines and Δr the bin width of the radial distribution function. θ refers to the Heaviside function and the sum over i ≠ j indicates summation over all pairs of particles i within the cyan and particles j within the green dashed line in Fig. 1. This function describes the normalized probability distribution of distances between skyrmions with the first peak corresponding to nearest neighbors. The boundary distribution function (Fig. 2b) was calculated in a similar manner via

$${g}_{{{{{{{{\rm{Bnd}}}}}}}}}(r)=\frac{1}{\rho }\frac{1}{L{{\Delta }}r}\mathop{\sum}\limits_{i}\left[\theta ({r}_{i}-r)-\theta ({r}_{i}-r-{{\Delta }}r)\right],$$
(2)

where r is the distance in the direction perpendicular to the boundary and L is the length of the boundary in the recorded video. The yellow line in Fig. 1 shows the position of the boundary determined from the position of average brightness between darker parts outside and brighter parts inside the rectangle. The boundary distribution function for the two boundaries exhibited a noticeable shift of 0.6 μm, which can be explained by the drifting of the sample compared to the background image used to determine the positions of the boundary. Because the exact amplitude of the drift cannot be determined directly using the Kerr videos, we shift both functions by 0.3 μm so that the first maximum is located at the same position assuming that the two boundaries interact identically with skyrmions. The boundary distribution function determined this way is further corrected by calculating the running average over eight neighbors in both directions for the region between the first peak and the first minimum of gBnd. This is done because the small peak at around 7.5 μm in the experimental function, likely caused by the pinning-induced quenched disorder, does not contribute to the real shape of the skyrmion-boundary interaction of this system. The quenched disorder may prevent skyrmions from leaving certain positions and can thereby introduce a larger number of skyrmions found at one specific distance to the boundary. The applied correction is shown in Fig. 2b and used for all simulations. Note that each skyrmion only contributes once per frame to the skyrmion–wall correlation function, while for skyrmion–skyrmion interactions pairs are considered, which is the main reason for the difference in statistics. To improve statistics for skyrmion-boundary interactions and obtain a more accurate interaction potential, a Kerr microscope with higher resolution and a larger field of view would be required.

Computer simulations

To derive interaction potentials for the given experimental setup we performed particle-based Brownian Dynamics simulations34,35,36 according to the Thiele equation43

$$-\gamma {{{{{{{\bf{v}}}}}}}}-G{{{{{{{\bf{z}}}}}}}}\times {{{{{{{\bf{v}}}}}}}}={{{{{{{{\bf{F}}}}}}}}}_{{{{{{{{\rm{SkSk}}}}}}}}}+{{{{{{{{\bf{F}}}}}}}}}_{{{{{{{{\rm{SkBnd}}}}}}}}}+{{{{{{{{\bf{F}}}}}}}}}_{{{{{{{{\rm{Random}}}}}}}}}$$
(3)

for each individual skyrmion with temperature kBT = 1, damping γ = 1 and a Magnus force amplitude G = 0.25 as already used in a previous work17 resulting in a skyrmion Hall angle of \({\theta }_{{{{{{{{\rm{Sk}}}}}}}}}=\arctan (G/\gamma )\approx 1{4}^{\circ }\). Note that the addition of a Magnus term in Eq. (3) will not alter the static behavior of the system35. FRandom are forces arising from thermal white noise satisfying the fluctuation-dissipation relation. The width of the simulated system and the skyrmion density are set according to experimentally determined values. Periodic boundary conditions are used in x-direction while the boundary in y-direction cannot be reached by skyrmions due to the repulsive skyrmion-boundary interaction. The timestep used in this simulation is Δt = 10−5 with data taken every 1000 steps. Due to the absence of a mass term in this description of skyrmions, integration of the equations of motion is performed directly using a standard Brownian integrator incorporating the Magnus force

$${{{{{{{\bf{r}}}}}}}}(t+{{\Delta }}t)={{{{{{{\bf{r}}}}}}}}(t)+{{{{{{{\bf{v}}}}}}}}(t){{\Delta }}t$$
(4)

with v(t) being determined by Eq. (3). This particle-based model of skyrmions comes with some limitations as it averages over effects like skyrmion deformations, internal excitation modes, size polydispersity, creation and annihilation of skyrmions, and also long-range and multi-particle interactions. Deformations and excitation modes, however, mostly occur under applied currents63 or fluctuating magnetic fields64,65 and are therefore less relevant to the analysis of static properties of a confined lattice. Furthermore, cases where skyrmions are no longer rotational isotropic but deform into elongated domains at different values of applied fields66 or lower densities67 are not described by this particle-based description. This model is considered a trade-off between the fidelity of treating the skyrmion structures, the complexity of the model and the robustness of the results with limited data typically available experimentally. Therefore, it is used most robustly for analyzing static properties like lattice formation and ordering or dynamics under small applied forces for systems with constant numbers of skyrmions. Its application is very versatile as it can easily be scaled to small skyrmion numbers in confinement17 or the analysis of large lattices28. Moreover, the model can be extended to a plurality of other quasiparticles, provided they have approximately isotropic interactions as, e.g., skyrmions found in different material systems such as frustrated magnets. We note that the model currently assumes circular spin structures and in our experiment, skyrmions are monodisperse in size and show no eccentricity in any preferred direction thus making the model applicable. For skyrmions that are deformed by external in-plane fields as realized previously61, the model would indeed require adaptation.

Iterative Boltzmann Inversion

Interaction potentials for skyrmion–skyrmion and skyrmion-boundary interactions are determined using the Iterative Boltzmann Inversion method57. The first simulation is undertaken with the potential of mean force

$${V}_{0}(r)=-{k}_{{{{{{{{\rm{B}}}}}}}}}T\ln g(r)$$
(5)

for the respective distribution functions g(r) and gBnd(r) where r denotes the distance between two skyrmions or between a skyrmion and the boundary, respectively. The potential of mean force already provides a first estimate for interaction potentials between particles and has already been used to calculate potentials based on experimental videos of colloidal systems68 and superconducting vortices69. Potentials are computed up to the first minimum of the distribution function. This choice leads to faster equilibration and better statistics, particularly for the derivation of skyrmion-boundary interactions. Using the potentials VSkSk,0 and VSkBnd,0, a Brownian Dynamics simulation is performed according to the above model to obtain the distribution functions g0(r) and g0,Bnd(r). Those are used to simultaneously update the interaction potentials by

$${V}_{i+1}(r)={V}_{i}(r)+\alpha {k}_{{{{{{{{\rm{B}}}}}}}}}T\ln \frac{{g}_{i}(r)}{g(r)}$$
(6)

where α ∈ (0, 1] is a parameter determining the speed of updating the potential which is chosen to be 0.2 in this work. Higher values of α lead to a faster convergence but also stronger oscillations around the correct potential and therefore less accurate results while a smaller α would slow down the convergence of the method. The two updated potentials are then used for the next simulation. This process is repeated until the potential converges, and in total 500 updating iterations each including Brownian Dynamics simulations were performed. A plot indicating convergence of the mean squared difference between the experimental and simulated pair distribution is shown in Supplementary Fig. 2.

Mathematically, there is a unique relation between the radial distribution function of a system and the pairwise interaction potential70 even though pair-correlation functions of similar interaction potentials may be hard to distinguish. Beyond this the convergence of results of the IBI method to real interaction potentials has been shown for several systems57,71. Important limitations of the IBI approach include that statistics must be sufficient to get a smooth distribution function and hence a potential without unphysical kinks. In addition, results calculated by this method are only reliable in a density range close to the density used for the determination of the potentials58.