Spin Torque Efficiency and Analytic Error Rate Estimates of Skyrmion Racetrack Memory

In this paper, the thermal stability of skyrmion bubbles and the critical currents to move them over pinning sites were investigated. For the used pinning geometries and the used parameters, the unexpected behavior is reported that the energy barrier to overcome the pinning site is larger than the energy barrier of the annihilation of a skyrmion. The annihilation takes place at boundaries by current driven motion, as well as due to the excitation over energy barriers, in the absence of currents, without forming Bloch points. It is reported that the pinning sites, which are required to allow thermally stable bits, significantly increase the critical current densities to move the bits in skyrmion-like structures to about jcrit = 0.62 TA/m². The simulation shows that the applied spin transfer model predicts experimentally obtained critical currents to move stable skyrmions at room temperature well, which is in contrast to simulations based on spin orbit torque that predict significantly too low critical currents. By calculating the thermal stability, as well as the critical current, we can derive the spin torque efficiency η = ΔE/Ic = 0.19 kBT300/μA, which is in a similar range to the simulated spin torque efficiency of MRAM structures. Finally, it is shown that the stochastic depinning process of any racetrack-like device requires an extremely narrow depinning time distribution smaller than ~6% of the current pulse length to reach bit error rates smaller than 10−9.

. Geometry of the used structure. The front and the back leads have dimensions of 30 nm × 90 nm × 3 nm. The dimensions of the DMI wire are l = 600 nm, w = 90 nm, and t = 1.8 nm. The diameter of the notches, as well of the spacer layer (t = 1.5 nm), the pinned layer (t = 9.0 nm), and the pinned lead (t = 6.0) is d = 60 nm.
www.nature.com/scientificreports www.nature.com/scientificreports/ later in detail. The nucleated skyrmion has magnetization parallel to the pinned layer, which leads to pinning of the skyrmion due to the strayfield of the pinned layer. In the following, a current pulse is applied to move the skyrmion into the center of the DMI wire. At this position, no significant strayfield due to the pinned layer is acting on the skyrmion, and the critical currents to move the skyrmion over the pinned site can be studied accurately.
In Fig. 2, the critical current in order to move the skyrmion over the pinning sites is studied. Two different geometrical realizations of the pinning sites are investigated. The diameter of the pinning sites is d = 120 nm (upper structure) and d = 60 nm (lower structure), respectively. The current density rise time is r = 11 GA/(m²ns). The actual current density is given in Fig. 2, and it is measured in the DMI wire between two notches. The current is applied at the back lead, and the potential is fixed to zero at the front lead. The smallest current density that was applied that moves the skyrmion toward and slightly into the pinning site is j crit = 0.40 TA/m². The minimum current density in experiments to move skyrmions in Pt/Co/Ta stacks is reported to be j crit = 0.20 TA/m² ref. 14 . However, in the simulations of Woo et al. 14 in the supplementary material, assuming even disorder in the film, the critical currents are reported to be j = 0.03 TA/m 2 . Hence, the simulations predict a critical current that is smaller by a factor of about 10 as the experimentally obtained value. As a consequence, one has to conclude that the considered SOT in the simulations of ref. 14 does not describe the experiment well. We attribute this discrepancy to the fact that the SOT of several multilayer stacks is partly compensated, as explained in the beginning of this section, and might not be the dominated driving force for skyrmions. Figure 2. Motion of the skyrmion as a function of the applied current density. The input current is applied at the front lead. The voltage is zero at the back lead. The dimensions of the wire are: l = 600 nm, w = 90 nm, and t = 1.8 nm. The diameter of the notches is d = 60 nm and d = 120 nm, respectively. The given current density is the average density within the wire holding the skyrmion at the center between notches. (left column) A positive current density is applied at the right lead; (right column) a negative current is applied at the right lead.
www.nature.com/scientificreports www.nature.com/scientificreports/ Much better agreement of the critical current can be obtained using the presented model. It does not have any free parameter and predicts a critical current that is only larger by about a factor of 2-3 compared to the experimental values of skyrmion motion in multilayers that have finite but in detail unknown pinning sites. For weaker pinning sites, the simulated current to move skyrmions can be decreased. In addition, it has to be noted that the simulations are performed at zero temperature, and the experiments are performed at room temperature. Hence, the simulated critical current should be overestimated compared to room temperature measurements 29,30 .
In contrast in ref. 13 , the Slonczewski term for spin-orbit torque predicts current densities of 0.001 TA/m² to overcome minor pinning sites and a current to overcome the largest pinning sites of 0.05 TA/m² 13 . For even larger pinning sites, no passing of the skyrmion was possible. Both values are significantly smaller than the experimental data to move skyrmions in multilayers 14 . Good agreement with experimentally obtained critical currents are obtained in disordered films 31 if a spin Hall angle of 0.34 is assumed 32 . However, in multilayers with repeated layer structure 32 , it is not obvious if the spin Hall angle can be used that is obtained for one repetition of the multilayer 32 due to diffusion processes between the layers, as discussed earlier.
The left column in Fig. 2 shows simulations where a positive current is applied at the back lead. Here, the skyrmion moves to the left. Interestingly, the critical current density to overcome the pinning site is smaller for the narrow pinning site (j crit = 0.62 TA/m² for d = 60 nm) compared to the pinning site with d = 120 nm (j crit = 0.67 TA/m²). The origin of larger pinning currents for a larger diameter of the pinning sites may be found in the larger extension of the pinning sites. As a consequence, for the larger radius, the skyrmion can not move freely between the pinning sites and does not gain velocity during the time the current is ramped up. Detailed analysis of this effect will be an interesting study for future work.
The used model, which couples the magnetization dynamics with the spin accumulation, self-consistently solves, besides the magnetization as a function of time, also for the electrical potential as a function of time, as well as the spin accumulation s. The additional spin torque field that is generated due to the spin accumulation s that follows from the solution of a differential equation (Eqs 9, 10 of ref. 22 ) is given by: ST s  and added to the effective field in the Gilbert equation of motion of the magnetization. Here, M s is the saturation magnetization, γ the gyromagnetic ratio, and  the Planck constant. The spin accumulation s, which depends on the magnetization, in turn, acts via the exchange strength J between the conducting electrons and magnetization as a torque term on the magnetization. The three components of the spin accumulation are shown in Fig. 3 for a current strength just before switching over the pinning center, as shown by the position of the skyrmion at the bottom of Fig. 3. It is important to note that all three components of the spin accumulation are unequal zero and in the same order of magnitude. This is in contrast to the simplified but often used model of Zhang and Li 33 . A similar model is the model of Thiaville et al. 34 that is equivalent to the model of Zhang and Li if the magnitude of the magnetization is constrained to remain constant. The used self-consistent spin drift-diffusion/micromagnetic model is equivalent to the model of Zhang and Li for the limiting case of a vanishing diffusion constant, D 0 = 0. The advantage of the Zhang and Li model compared to the used self-consistent spin drift-diffusion/micromagnetic model is the simplified calculation of the spin-transfer torque, which can be expressed explicitly as a function of the magnetization. However, the question must be asked if the approximation leads to different results, and the combined solution of both magnetization and spin accumulation is required for accurate results. Hence, in the following, we study the critical current densities in the limit of D 0 = 0, as shown in Fig. 4. The analysis of the simulation data shows that the skyrmion precesses on an elliptical orbit, as well as changes its size periodically as a function of time. In Fig. 4, a current density is applied, which is slightly smaller than the critical current density to overcome the pinning site. The critical current to overcome the pinning site with D 0 = 0 is j crit = 1 TA/m², which is about a factor of 1.6 larger than for the realistic value of D 0 = 10 −3 m²/s. This clearly shows the importance to solve for the full self-consistent model and the expected error in the prediction of the critical current for the limit of D 0 = 0 (Zhang Li model) in this system.

Thermal Stability and Spin Torque Efficiency
As mentioned in the introduction, a prerequisite for any potential storage application is the thermal stability of the stored information. The average lifetime τ of magnetic states can be expressed within the framework of the transition state theory (TST) as: where ΔE is the energy barrier between two stable states and f 0 is the attempt frequency. A wide range of attempt frequencies from several Mhz to 10 THz were reported 19,[35][36][37][38][39][40] . In ref. 39  was reported as a result of a simulation based on TST. In order to obtain thermally stable bits in hard disk recording, an energy barrier of about ∆ ≥ E k T to k T 40 50 B B 300 300 is a common requirement. In the following, we present simulations of the energy barrier using the string method, taking into account all energy terms, including the DMI interaction 41,42 . As an input of the string method, we use the magnetization states of current driven skyrmion motion over one energy barrier, as shown in Fig. 2. This path leads to the two stable magnetic states on each side of the pinning site, as shown in Fig. 5, by the images indexed with i = 0 (initial state 1) with i = 19 (initial state 2). The simulation in the right column of Fig. 5 corresponds to a geometry with l = 500 nm, w = 75 nm, t = 1.5 nm, and d = 50 nm. Here, all dimensions of the previous simulations are scaled by a factor 0.83. The reason is that we aim to investigate and present a structure as small as possible, which still (2019) 9:4827 | https://doi.org/10.1038/s41598-019-41062-y www.nature.com/scientificreports www.nature.com/scientificreports/ supports pinning and depinning of a skyrmion without its annihilation. The images indexed i = 2 to i = 18 show magnetic states along the minimum energy path (MEP). This path is the most probable path that is triggered by thermal fluctuation if no external field or current is applied. The energy along this minimum energy path (MEP) is shown in Fig. 6 (black line) as a function of path index i. The reason why the stable state i = 0 has lower energy than the state with i = 19 is due to the strayfield of the pinned layer, which is magnetized in the direction of the skyrmion and stabilizes the skyrmion. The smaller energy barrier (barrier 2) of the two possible barriers (barrier 1: moving the skyrmion from left to right, barrier 2: moving the skyrmion from right to left) is ΔE = 16.4 k B T 300 , which does not provide sufficient thermal stability for stable bits.
The situation becomes even more severe if the dimensions of the structure are scaled to smaller dimensions: l = 450 nm, w = 67.5 nm, t = 1.35 nm, and d = 45 nm. Here, all spatial coordinates are scaled by a factor scale size = 0.75 compared to the original structure of Fig. 1. For this structure, the MEP is shown in Fig. 5 (left column). It can be seen that the skyrmion is annihilated and nucleated again at the boundary along the MEP, which has recently been reported independently 16,17 . It is reasonable that in principle, for finite-size systems, topologically charged structures may simply be driven out of the sample, leading to a homogeneous state with zero topological charge 43,44 . The annihilation of a skyrmion at a boundary was already shown during the current driven skyrmion motion in Fig. 2 (bottom right), indicating that this structure can change its topological charge due to current driven motion without forming Bloch points. Different definitions of topological protection are used in the community. According to ref. 17 , topologically protected means "there is an energy barrier separating the transition of a system from one topological state to another. " According to ref. 45 , protection means that "the spin configuration can not be twisted continuously to result in a magnetic configuration with different S (for example, a uniformly magnetized one), " where S is the topological charge. According to the wide definition of topological protection of ref. 17 , the boundary annihilation is still a topologically protected process, since an energy barrier must be overcome. According to the definition of ref. 45 , the boundary annihilation is not topologically protected, since the skyrmion can be moved continuously to a magnetic configuration with different S, as also reported in ref. 16 . A further very interesting observation is that the energy barriers in the two investigated structures are significantly different. The smaller structure showing annihilation of the skyrmion via the boundary displays a smaller energy barrier for the same initial path. In order to investigate this effect in more detail, the MEP with annihilation, obtained from the simulation with scale size = 0.75, is used as the input path for the simulation with scale size = 0.83, which originally did not lead to skyrmion annihilation. Hence, we aim to calculate the energy barrier of annihilation for the larger system and compare it with the energy barrier to overcome the pinning site. www.nature.com/scientificreports www.nature.com/scientificreports/ As a surprising result, the energy barrier for annihilation (ΔE = 13 k B T 300 ) is smaller than the minimum energy path over the pinning site (ΔE = 17 k B T 300 ). Interestingly, the initial path of the skyrmion over the pinning site is a stable local minimum energy path. But, there exists at least one local minimum energy path (the MEP that shows annihilation of the skyrmion) with an even lower barrier. Most striking about this effect is that the energy barrier at the pinning site is larger than the energy barrier for annihilation. For the investigated number of multilayer repetitions, both energy barriers are significantly too small for applications. In order to study the energy barrier over the pinning site for different geometries of the constriction, the lateral dimensions of the wire are scaled, and the barrier is calculated. In Fig. 7, the energy barrier is shown for structures where the lateral spatial coordinates are scaled by a factor scale size compared to the original structure of Fig. 1, with the Co thickness kept constant at t Co = 1.5 nm. Again, a significant reduction of the energy barrier is obtained for scale size < 0.75, since the skyrmion is annihilated for this size, despite the initial path guiding the skyrmion through the constriction. For scale size = 0.83, the energy barrier is significantly higher, since the skyrmion is not annihilated. It can be seen that there exists an optimal lateral dimension of the structure for scale size = 0.83, which corresponds to l = 500 nm, w = 75 nm, and d = 50 nm, yielding the highest energy barrier. Hence, for this given film thickness t, the energy barrier can not be further increased to make it stable at room temperature.
The annihilation path is via boundary annihilation, which is shown in Fig. 8 for a structure with l = 400 nm, w = 60 nm, t = 1.20 nm, and d = 40 nm. It is interesting to note that the saddle point configuration is a configuration where the entire skyrmion is still well located within the wire. If the skyrmion partly annihilates, the energy decreases. As it can be seen in Fig. 8 (i = 21), the spins between the boundary and the skyrmion rotate by 180° to connect the skyrmion core with the boundary. This is a smooth transition of the spins from the saddle point to states in which the skyrmion exits the structure. Clearly, no Bloch points are formed during this annihilation process. Due to the relative large skyrmion size, the magnetization angle between neighboring spins is small, and the micromagnetic continuum approach is well suited. Hence, the energy barrier of ΔE = 7 k B T 300 can be well calculated with micromagnetic simulations. Hence, there is no need to perform atomistic simulations or multiscale simulations for this kind of annihilation process at the boundary.
As we have already shown, the energy barrier can not be further increased by changing the pinning site dimensions. One way to increase the energy barrier of the device is to increase the number of repetitions of the layer stack 46 . www.nature.com/scientificreports www.nature.com/scientificreports/ In the following, the energy barriers for three different paths are investigated as a function of the total thickness of the Co layer: (i) the path over the pinning sites from one bit position to the next one, with the initial path taken from a simulation in which the structure is moved by current (see Fig. 5 in the right column); (ii) the path of annihilation, where the skyrmion is annihilated at a pinning site, with the initial path taken from the simulation of Fig. 5 (left); and (iii) the path of annihilation, where the skyrmion is annihilated at the top flat edge. Here, the initial path is constructed by shifting the position of the skyrmion in the +y direction. It has to be mentioned that a detailed study for the maximum thickness before inhomogeneous magnetization in the z-direction are formed will require simulations where the reduced exchange strength between the Co layers due to the Ta and Pt layer are taken into account 47 .  www.nature.com/scientificreports www.nature.com/scientificreports/ For the wire with l = 500 nm, w = 75 nm, and d = 50 nm, the energy barriers of these three paths are shown by the plots with circles in Fig. 9. For both paths that lead to annihilation, the energy barrier increases approximately linearly as a function of the film thickness (layer repetitions). For the path over the pinning site, the energy barrier does not exactly increase linearly as a function of the thickness, which is attributed to a strayfield effect. Since the skyrmion changes size via moving through the pinning site, the effect of the strayfield, which changes as a function of thickness, becomes important. The simulations indicate that a film thickness of t Co > 4.8 nm (8 stack repetitions) is required in order to obtain energy barriers larger than 50 k B T 300 for the structure with l = 500 nm, w = 75 nm, and d = 50 nm.  www.nature.com/scientificreports www.nature.com/scientificreports/ For comparison, the energy barrier for annihilation at a smaller defect (l = 400 nm, w = 60 nm, and d = 40 nm) is also shown in Fig. 9 (black line, rectangles). It can be seen that the energy barrier for the annihilation processes decreases with the pinning diameter. The largest energy barrier for annihilation is obtained at the flat edge.
Finally, we calculate the spin torque efficiency for the structure with l = 600 nm, w = 90 nm, and t = 1.8 nm, which can be derived from the energy barrier and the critical current: Here, it might be important to note that the spin torque efficiency is approximately independent from the film thickness, since both the energy barrier and the critical current depend approximately linearly on the film thickness. Hence, this estimate is also valid for thicker films. The spin torque efficiency of a metallic MRAM junction simulated with the same spin diffusion approach is in the similar order, η = ΔE/I c = 0.16 k B T 300 /μA 48 . The spin torque efficiency for experimental CoFeB-MgO-based tunnel junctions with diameters below 30 nm are in the range between 1 and 10 k B T/μA, as reported in ref. 49 . Hence, this study indicates that the spin torque efficiency of the investigated skyrmion structures is in the same order of magnitude or smaller compared to MRAM structures if spin transfer torque is the leading torque.

Bit Error Rate of Racetrack-like Devices
Another important topic for any skyrmion-like device is the requirement for reliable depinning of the magnetic structure from the pinning center. If current pulses are applied in the investigated skyrmion racetrack-like structure, the depinning process is a stochastic process due to finite temperature and will have further distribution in time due to imperfections of different pinning centers. The distribution of the depinning time of field-driven and current driven domain walls is investigated in various studies [50][51][52][53] . For MRAM structures in the dynamical switching regime, the standard deviation of the switching time is in the order of 10% of the switching time 54 . In ref. 55 , it can be seen that the skyrmion positions are statistically distributed, and after applying current pulses, some skyrmions are moving, but others do not move due to larger pinning potentials 55 . By means of computer simulations, the distribution of velocities and skyrmion positions is studied for disordered films, showing significant challenges for applications 31 .
In applications, the current pulse must have a sufficient duration and strength to move the topological structure over the first pinning site, as shown in Fig. 10(c). However, if the current pulse is too long, the probability increases that the structure moves too far over the second pinning site (Fig. 10d). Hence, the optimization of the correct pulse duration is essential to avoid written-in errors. Let us assume a current pulse with a particular strength leading to 50% depinning probability at a pulse length t a . Let us further assume that the depinning time is normally distributed with a standard deviation of σ. The probability that the structure depins within a time t pulse is given by: The probability that depinning occurs can be increased if the current pulse duration is increased, but the probability that the topological object moves too far over the next pinning site increases, too. The probability that the structure moves over two pinning sites at time t′ in the time interval dt′ is The total probability that within the time of the pulse t pulse the skyrmion moves over the second pinning site is given by: a a which leads to: depinning a pulse depinning a pulse ,2 The probability of a successful write process can be expressed as 1 minus the probability that the structure moves across at least two pinning sites minus the probability that the structure does not move across one pinning site. If the racetrack-like device consists of N bits, the before-mentioned process must be successfully repeated for N times to move the skymrion from the input to the output. Hence, the success probability is given by:  Figure 11 shows P success as a function of the pulse duration t pulse for σ = 0.1 ns and t a = 1.0 ns. It can be seen that for too short a t pulse, the success rate is small, since the topological structure can not move over the first pinning site, as shown in Fig. 10(b). For pulses that are too long, the success rate is small again due to the high probability of moving over two pinning sites. If the pulse duration is t pulse = 1.4 × t a for N = 1, the highest success probability P success is obtained. The bit error rate = − BER P 1 success is plotted as a function of σ in Fig. 12 for different N. For N = 1, it shows that increasing sigma from σ = 0.1 ns to σ = 0.2 ns changes the success rate significantly from P success = 3.3 × 10 −5 to P success = 3.7 × 10 −2 . To meet the write and the read bit error rate requirements of STT-MRAM, which is BER < 10 −9 ref. 56 , a standard deviation of the depinning time σ < 0.06 t a is required for N = 10. In principle, the time t a can be increased by increasing the distance of adjacent bits or decreasing the speed of the topological object. Both possibilities lead to unwanted properties for data storage, restricting either data rate or data density. Hence, detailed optimizations and accurate estimates of σ, for example, by Langevin simulations are required to access the feasibility of racetrack-like devices. A possibility to increase P success might be to apply current pulses with different strengths. First, a small strength just to move the structure to the next pinning site could be used. If the structure is located next to the pinning site, a high current pulse to overcome the pinning site can be applied. This clock cycle might lead to an increase of t a but does not increase σ.
Finally, we also aim to investigate the requirements for a BER of 10 −9 ; if not, such well-defined pinning sites are used, as in Fig. 1, but the film consists of a disordered material with random pinning sites 31 . Let us assume the www.nature.com/scientificreports www.nature.com/scientificreports/ racetrack device stores N bits with a center-to-center distance of the skyrmion positions of d. Let us assume that after each current pulse, the skyrmion position can be described by a normal distribution with a standard deviation of the new position σ d . After N pulses, a skyrmion moves through the entire structure. If we assume a Gaussian process, the standard deviation of the final position is σ N d . In order to achieve the required bit error rate of 10 −9 , six times the final standard deviation of the position must be smaller than the bit length, In both examples (disordered film, well-defined pinning sites), concepts that relax the requirements of the distributions, such as structures reported in refs 57,58 , are required.

Discussion
In this study, the stability of skyrmions, as well as critical currents to move skyrmions, is investigated. Since, for a storage application, the skyrmion positions must be well defined and thermally stable, skyrmions must be pinned at certain positions. Usually, the presence and the absence of a skyrmion represent the stored information. In this work, pinning sites are introduced by geometrical constrictions to define bit positions. Critical currents, as well as energy barriers to overcome these pinning sites, are investigated. It is shown that the required pinning sites for storage applications lead to critical current densities to write bits that are in a similar range as micromagnetically simulated critical currents in MRAM devices. A very interesting effect is found that the energy barrier to overcome these pinning sites is larger than the energy barrier for skyrmion annihilation in the investigated structure. For the given material parameters, the isotropic annihilation showed the largest energy barrier, as discussed in the method section 'Minimum energy path' . For the investigated multilayer stack, more than 7 repetitions are required to obtain sufficiently high barriers. In order to achieve bit error rates as required by MRAM structures, innovative new concepts will be required. Concerning a practical point of view, it is also worth mentioning that bits in hard disks (45 nm × 10 nm) 59 and FLASH storage at 2.77 Tb/in² 60 (~10 nm × 10 nm) are significantly Figure 11. Dependence of the probability for successful depinning over one pinning site and not depinning over two pinning sites as a function of N write pulses t pulse . The average depinning time is t a = 1 ns. Figure 12. Dependence of bit error rate to move a skyrmion from the input to the output for a racetrack device with N bits as a function of σ. For each σ, the pulse length t a was optimized to reach the highest P success . smaller than the thermally stable skyrmion structures 55,61 and the required distances between them reported in this work. As a consequence, applications and markets need to be found in which the size of the magnetic bit is not the leading objective and multiple bits per cell in RAM-like devices lead to an increased performance 62 , such as in racetrack-like devices for MRAM-like structures. The STT-MRAM cell area is dominated by the access transistor (leading to 40 F² scaling, where F is the minimum feature size 63,64 ), and the MTJ is much smaller in size 65 . Hence, multiple bits per cell will lead to enhanced storage density at the cost of latency time. Special architecture level optimization, such as bit-interleaved, is suggested in order to overcome this problem 63 .
The financial support by the Austrian Federal Ministry of Science, Research and Economy and the National Foundation for Research, Technology and Development, as well as the Austrian Science Fund (FWF), under Grant Nos. F4112 SFB ViCoM, I2214-N20 is acknowledged.

Solution of the coupled spin drift-diffusion equation and micromagnetic equation. The dis-
cretization of the continuum model is performed using a hybrid finite element/boundary element method, as presented in detail in ref. 21 . The total energy of the magnetic system is composed of the exchange energy, the magnetostatic energy, the anisotropy energy, and the Zeeman energy 66 . In order to take into account the interface DMI, the following energy contribution is added to the total energy E tot s eff tot The right-hand side of Eq. (9) is the Gâteaux derivative of the total energy E m ( ) tot in the direction v. Mass lumping is used in order to calculate the effective field H eff on the node points from the right-hand side of Eq. (9). The approach using the Gâteaux derivative of the total energy has the advantage that the boundary conditions for the corresponding field terms are considered in a natural fashion and no explicit treatment is required. In contrast, if the effective field on the node points is calculated by solving where the second term is derived from the functional derivative of the total energy DMI DMI s z z the boundary conditions for the normal derivative of the magnetization at the boundary due to the DMI interaction and the exchange interaction, z must be explicitly considered.
Effective media model. In the following, we will describe the methodology that is used to model the magnetic multilayers. Physically, the multilayers consist of magnetic layers with different magnetic properties. As it is common in micromagnetic simulations, atomic features are only considered in a continuum sense (for example, the layer structures of L 10 layer materials are not resolved in micromagnetic simulations).
In the following, we start from the magnetic properties of two magnetic layers, as shown in Fig. 13 (left). The layers denoted with the label 2 may represent non-magnetic layers, such as Ta or Pt. The layer with label 1 may denote the magnetic Co layer. In the following, we aim to approximate the layered model (left) with an effective media model (right) with appropriate parameters. www.nature.com/scientificreports www.nature.com/scientificreports/ Let us start to obtain an effective saturation magnetization J s,eff for the effective media model. In order to obtain the same magnetic moment μ as in the layered model, one requires: Hence, one gets: s eff s s , , 1 1 ,2 2 Hence, the Zeeman energy is the same for these two models, independent from the direction of the magnetization.
In order to derive an effective anisotropy constant, K eff 1, , it is worth noting that this can not simply be obtained by the same argument as used above. Let us derive K eff 1, by calculating the total energy in parallel and antiparallel direction, which must be the same for the two models. Since the Zeeman energy and the exchange energy are the same for both configurations, these contributions are not included in the following calculation: Here, the second term on the left-hand side considers the shape anisotropy of the effective media model. The corresponding shape anisotropies of the two layers in the layered model are considered by the second and fourth terms on the right-hand side of Eq. (15). Since we assume that the film thickness is significantly smaller than the width of the layers, the strayfield of one layer to the adjacent layer can be neglected. As a consequence, the layers do not interact with each other due to strayfield interactions, and the demagnetizing field must be considered separately for each layer. K eff 1, is then simply obtained by: For the exchange energy, the simple average approach is used as: Here, it should be noted that the approximation of Eq. (17) is well justified for the exchange energy contribution in the direction within the plane. For variations of the magnetization in the z-direction (perpendicular direction), due to the different exchange constants within the layers, this approximation will not be well suited. However, within this work, the film thickness is thinner than the domain wall width in the z-direction, and hence, the approximation will be well suited, since no inhomogeneities in the z-direction are expected. The eqs Eqs (14)(15)(16)(17) agree with the independently developed effective media model of ref. 14 .
In Table 1, the used material parameters are summarized for a magnetic layer (1) with the material parameters given by column (c). The non-magnetic layer (2) has a thickness of 2 nm, which can be 1 nm Pt and 1 nm Ta. In the layer (2), all magnetic properties are assigned to have value of zero. Six repetitions of the layer stack are assumed. Column (a) denotes the resulting material parameters of the effective media model (b) denotes the media model where only the magnetic layer is simulated and stacked above each other.
The validity of the effective media model is presented in Figs 14 and 15. In Fig. 14, it is shown that the hard axis loop of an extended but finite film agrees very well for the three investigated models. The model (a) is the model of the effective media. Model (b) is a model where only the Co layers are simulated and stacked directly in contact above each other. Since no domain walls are expected within the film thickness, this is a good approximation of the structure. Model (c) is a model where 6 Co layers are simulated with air in between, which represents the non-magnetic Ta and Pt layers.
The energy barriers of these three models are compared in Fig. 15. Whereas the barriers of model (a) and model (b) perfectly agree, the model with the 6 Co layers that are not exchange-coupled shows a significant difference. For model (c), the skyrmion gets annihilated at the boundary. The origin may be in the unrealistic weak coupling between the layers in the model that is only triggered by strayfield interactions and may lead to an individual motion of the skyrmions within the layers. In reality, we expect a significant exchange coupling Effective model (a) 6 con. Co layers (b) 6 sep. Co layers (c) www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 14. Comparison of hard axis loop (field applied in-plane) of an effective media model (effective model), a model where only the Co layer is simulated and directly stacked above each other with full exchange coupling (6. Con Co layers), and a layered resolved model where between the Co layers 2 nm air is assumed, since Pt and Ta are assumed non-magnetic. www.nature.com/scientificreports www.nature.com/scientificreports/ between the Co layers. In ref. 47 , the coupling strength in Co(4 Å)/Pt(t PtÅ)/[Co(4 Å)/Pt(6 Å)] 2 as a function of the Pt coupling layer is investigated. For a spacer layer of 2.5 nm, a coupling strength of 0.5 mJ/m² is reported, which results in an exchange constant of A = 0.6 pJ/m. By comparing the equilibrium sizes of the skyrmions of the different models (Fig. 15), one can observe slightly different sizes. The reason is that due to the layered structure (layered resolved model), an additional induced anisotropy occurs due to dipolar interaction (see, for example, supplementary material of ref. 14 ). This effect is well known in micromagnetics 68 and is included in an additional contribution to the anisotropy constant (for example, in FePt, the L 10 phase consists of one layer Fe and one layer Pt). If the anisotropy constant is measured, this dipolar induced anisotropy term is automatically included. Hence, in the effective media model, no special treatment is required. In contrast, if the measured anisotropy constant is used for the layered resolved model, one has compensated for this effect if the measured value of the anisotropy is used.

Minimum energy path.
For the calculation of the minimum energy path, the string method is used 41,42 . An advantage of the string method (details of the implementation for magnetic systems can be found in ref. 18 ) over the nudged elastic band 69,70 method is that it does not require any modifications on the existing micromagnetic simulator. The entire method was implemented within the python interface of the micromagnetic software magnum.fe 71 . In order to relax the images of the string method, the Landau-Lifshitz Gilbert equation is integrated without the gyromagnetic precession term for 20 ps. After this time, the images are equally distributed along the string. As stopping criterion, the change of the energy barrier is used, which must be smaller than 10 −23 J. In ref. 72 , two completely independent atomistic implementations of the barrier for isotropic skyrmion annihilation are implemented. The barrier obtained by the GNEB method (ΔE = 4.416 × 10 −20 J) agrees well with the simple string method (ΔE = 4.421 × 10 −20 J).
The energy barrier as a function of the mesh size (denoted with FE) is shown in Fig. 16 for boundary annihilation and isotropic annihilation via a Bloch point for a system of size 90 nm × 90 nm × 0.6 nm, which represents one 0.6-nm-thick Co layer (1 repetition of the layer stack). It can be seen that the boundary annihilation barrier does not depend on the mesh size, since no Bloch point is formed. Hence, it can be reliably calculated with micromagnetics. As expected, for the isotropic annihilation, a mesh size dependence is obtained due to the formation of a Bloch point. For comparison, we also performed atomistic simulation using 3D Heisenberg spins within one layer of Co atoms with a lattice constant of a = 0.5 nm using the approach of ref. 18 , denoted with the label (atom). It can be seen that the results agree well with a finite element method, where the mesh size is chosen in the order of the lattice parameter. This agreement is expected for boundary annihilation, since no Bloch point is formed. However, the agreement is surprisingly good for isotropic annihilation, since for this path, the micromagnetic prerequisite of continuous magnetization as a function of space is not fulfilled. Figure 16. Comparison of energy barrier for boundary annihilation and isotropic annihilation via a Bloch point as a function of the mesh size for a system of size 90 nm × 90 nm × 0.6 nm, which represents one 0.6-nmthick Co layer (1 repetition of the layer stack). The data points denoted with (FE) are obtained from the finite element method, as used within the paper. For comparison, atomistic simulations are also performed with a lattice parameter a = 0.5 nm using the approach of ref. 18 .