Path sampling for lifetimes of metastable magnetic skyrmions and direct comparison with Kramers' method

TL;DR

This paper compares Kramers' theory and path sampling methods for calculating the collapse rates of magnetic skyrmions, highlighting the importance of entropic contributions and rate prefactors in complex systems.

Contribution

It demonstrates the agreement between Kramers' method and forward flux sampling in skyrmion lifetime estimation and emphasizes the significance of entropic effects.

Findings

Good agreement between the two methods.

Collapse rate varies by three to four orders of magnitude with magnetic field changes.

Entropic contributions significantly affect skyrmion lifetimes.

Abstract

We perform a direct comparison between Kramers' method in many dimensions -- i.e., Langer's theory -- adapted to magnetic spin systems, and a path sampling method in the form of forward flux sampling, as a means to compute collapse rates of metastable magnetic skyrmions. We show that a good agreement is obtained between the two methods. We report variations of the attempt frequency associated with skyrmion collapse by three to four orders of magnitude when varying the applied magnetic field by 5$\%$ of the exchange strength, which confirms the existence of a strong entropic contribution to the lifetime of skyrmions. This demonstrates that in complex systems, the knowledge of the rate prefactor, in addition to the internal energy barrier, is essential in order to properly estimate a lifetime.

Equations

$$k={\tau }^{-1}=f_0{e}^{-\beta \Delta E},$$

$$\begin{split}&E=-J_{\text{ex}}\sum_{<ij>}\mathbf{m}_{i}\cdot\mathbf{m}_{j}-\sum_{<ij>}\mathbf{D}_{ij}\cdot\big{(}\mathbf{m}_{i}\times\mathbf{m}_{j}\big{)}\\ &-K\sum_{i}m_{z,i}^{2}-B_{z}M_{s}\sum_{i}m_{z,i},\\ \end{split}$$

$$k=f_0^{\prime }{e}^{-\beta \Delta F},$$

$$k_{\text{FFS}}={\Phi }_{A,0}\prod _{i=0}^{n-1}P({\lambda }_{i+1}|{\lambda }_i),$$

Full text

The rate of decay of metastable states is an ubiquitous problem in physics. Thermal activation processes across an energy barrier are found within fields as diverse as solid state physics (Josephson junctions), chemical reactions, electrical circuit theory (phase-locked loops), laser physics, and magnetization switching in ferromagnets Coffey et al. (2001); Risken (1996). In the context of magnetic data storage, information is stored in the form of 0 and 1 bits, corresponding to uniformly magnetized grains pointing along opposite directions. New challenges arise in the necessity to design small magnetic structures capable of retaining a given state against fluctuations for an average lifetime of 10 years at room temperature Wild et al. (2017). The ability to precisely predict that lifetime is therefore crucial. The rate of such thermally activated processes can be described by an Arrhenius law Hänggi et al. (1990),

[TABLE]

in which $\beta=(k_{B}T)^{-1}$, $\Delta E$ is the internal energy barrier, and the prefactor $f_{0}$, commonly referred to as the attempt frequency, corresponds to a fundamental fluctuation rate. Estimating the stability of magnetic structures is often synonymous with accessing internal energy barriers, while assuming a typical value of the prefactor in the GHz range Weller and Moser (1999); Chen et al. (2010); Lederman et al. (1994); Cortés-Ortuño et al. (2017). Hence, it is generally accepted that $\beta\Delta E\sim 50$ at room temperature is a sufficient and necessary condition in order to achieve the desired stability.

In recent years, magnetic skyrmions Bogdanov and Yablonskii (1989); Bogdanov and Hubert (1994a) have emerged as potential candidates for spintronics applications in data storage and logic devices Fert et al. (2013); Iwasaki et al. (2013); Zhou and Ezawa (2014); Müller (2017); Hsu et al. (2017). Magnetic skyrmions are particle-like spin textures carrying an integer topological charge. They are stabilized by the introduction of a characteristic lengthscale in a system via competing interactions. In particular, the existence of chiral skyrmions in non-centrosymmetric bulk magnets and thin magnetic films with broken inversion symmetry is made possible by the Dzyaloshinkii-Moriya interaction (DMI) Dzyaloshinskii (1958); Moriya (1960); Fert and Levy (1980); Crépieux and Lacroix (1998). Isolated skyrmions typically live on the ferromagnetic (FM) background as metastable excitations, but, under the effect of thermal fluctuations, will eventually collapse back to the uniformly magnetized state. The problem of their thermal stability has so far yielded vastly different theoretical predictions depending on the system of interest Rohart et al. (2016); Desplat et al. (2018); Bessarab et al. (2018), particularly concerning the order of magnitude of the attempt frequency. Experimentally, extreme variations of $f_{0}$ were observed for small variations of the applied magnetic field in the case of the decay of the skyrmion lattice Wild et al. (2017). The apparent lack of consensus between the results is in part due to the difficulty in calculating rate constants of rare events. For structures with lifetimes well above the nanosecond range, direct Langevin simulations García-Palacios and Lázaro (1998), where one integrates the stochastic dynamics of the spin system at each timestep, becomes unrealistic. In that case, a possible approach is the use of a form of reaction rate theory Hänggi et al. (1990); Coffey et al. (2001), which allows a direct calculation of the rate prefactor by considering details of the fluctuations about the metastable state $A$ and the saddle point (SP) $S$ along the reaction coordinate. Numerical implementations of this method Suess et al. (2011); Fiedler et al. (2012) combined with a path finding scheme Dittrich et al. (2002); Bessarab et al. (2015) have previously been used to obtain lifetimes of magnetic skyrmions Bessarab et al. (2018); Desplat et al. (2018); Haldar et al. (2018); von Malottki et al. (2019). While this is a computationally optimal solution, reaction rate theory is based on many assumptions concerning the damping regime, the energy landscape, and the density of states of the system. Additionally, whenever we are faced with several mechanisms for a single process, we can only assume that the mechanisms are independent in order for the rates to add up, which may not hold. We are also faced with the questions of whether higher-order saddle points should contribute to the rate, and whether Eigenmodes with small Eigenvalues should be treated as Goldstone modes. An alternative method is therefore required in order to validate previous results. For that purpose, we turn to forward flux sampling (FFS) Allen et al. (2005, 2006a, 2006b, 2009); Borrero and Escobedo (2008). FFS is a path sampling method that was initially developed to simulate rare switching events in biochemical networks. It has since then been applied to a wide range of rare event problems Allen et al. (2009). In particular, it was used to obtain magnetization switching rates in magnetic microstructures Vogler et al. (2013, 2015). FFS was shown to be significantly more efficient than brute force direct Langevin simulations, while enabling an exploration of phase space free from assumptions. In this article, we demonstrate the application of the FFS method to the computation of collapse rates of metastable magnetic skyrmions far away from the system’s boundaries, and we compare the results with predictions from Langer’s theory, as well as with direct Langevin simulations whenever it can realistically be achieved.

We simulate $N$ magnetic spins of constant amplitude on a two-dimensional square lattice. We use the classical Heisenberg model Hamiltonian,

[TABLE]

where $J_{\text{ex}}$ is the strength of the isotropic exchange coupling, $\mathbf{D}_{ij}$ is the interfacial Dzyaloshinskii vector, $K$ is the perpendicular uniaxial anisotropy constant, $M_{s}$ is the saturation magnetization, and $B_{z}$ is the perpendicular applied magnetic field. Exchange interactions are restricted to first nearest neighbors. We introduce the reduced parameters: $d=\lvert\mathbf{D}_{ij}\lvert/J_{\text{ex}}$; $k=K/J_{\text{ex}}$ ; $b_{z}=B_{z}M_{s}/J_{\text{ex}}$, and we set $(d,k)$ = (0.36, 0.4), which allows the existence of small Néel skyrmions solutions at zero-field that only span over about 7 lattice sites in diameter (state $A$ in Fig. 1a) Heo et al. (2016). At low temperature, the skyrmions do not exhibit translational invariance on the lattice – i.e., no Goldstone modes of zero-energy fluctuations – but instead experience pinning at particular lattice positions. The applied field is oriented opposite to the skyrmion’s core and has a destabilizing effect. The rest of the material parameters correspond to Pt/Co/AlOx samples Miron et al. (2011); Thiaville et al. (2012); Rohart and Thiaville (2013) and are given in the Supplemental Material (SM) sm_ . We simulate an infinite system by setting periodic boundary conditions, which eliminates cases where the skyrmion escapes at the edges Uzdin et al. (2017); Lobanov et al. (2016); Bessarab et al. (2018); Desplat et al. (2018); Cortés-Ortuño et al. (2017).

We firstly relax the paths of minimum energy that lead to the skyrmion collapse on the energy landscape, and identify the saddle point along the path via the geodesic nudged elastic band method with a climbing image Bessarab et al. (2015); Henkelman et al. (2000). The prefactor in Eq. (1), $f_{0}$, is then calculated via an extension of Kramers’ method Hänggi et al. (1990) to many dimensions, namely, Langer’s theory for the decay of metastable states Langer (1969) adapted to magnetic spin systems Braun (1994); Coffey et al. (2001); Desplat et al. (2018). The theory is set in the intermediate-to-high damping regime. It yields no temperature dependence of $f_{0}$ if there are no Goldstone modes, or the same number of Goldstone modes at $A$ and $S$. Eq. (1) may be rewritten in terms of the change in Helmholtz free energy $\Delta F$,

[TABLE]

where $f_{0}^{\prime}$ is a new prefactor, and $\Delta F=\Delta E-T\Delta S$. $\Delta S$ corresponds to the change in configurational entropy undergone by the system upon reaching the saddle point. Details on the calculation of $f_{0}$ and $\Delta S$ can be found in sm_ ; Desplat et al. (2018). Note that $\Delta S$ is defined for stable modes of fluctuations, whereas $f_{0}$ takes into account both stable and unstable contributions. We report two distinct collapse mechanisms. In one case, which we refer to as mechanism 1, the skyrmion shrinks in size while its core coincides with a lattice site, and the core-spin flips past the saddle point Desplat et al. (2018). Alternatively, if mechanisms 2 is realized, the skyrmion core may shift to an interstitial position before uniformly shrinking Bessarab et al. (2015). These two mechanisms involve distinct saddle points, that we respectively refer to as $S_{1}$ and $S_{2}$ [Fig. 1]. $S_{2}$ is found above $S_{1}$ on the energy surface, by advancing along the Eigenbasis coordinate associated with a translation mode. If the translational modes at the saddle points are not Goldstone modes, $S_{1}$ and $S_{2}$ should be treated as distinct states associated with different activation rates, namely $k_{1}$ and $k_{2}$. When the metastable skyrmion is pinned at an interstitial position, the realization of mechanism 1 requires the core to firstly shift onto a lattice site [Fig. 1b]. The way the relaxed skyrmion sits on the lattice depends on its equilibrium size and its commensurability with the underlying lattice, such that there exists only one type of stable skyrmion state per field value, although the skyrmion can be indistinguishably located at either of the $N$ possible sites. $S_{1}$ is a first-order saddle point, with a single unstable mode corresponding to the breathing of the skyrmion Desplat et al. (2018). At $S_{2}$, three unstable modes are found – the unstable breathing mode, and two unstable modes of translation – resulting in a third-order SP. As we increase the field, the stable skyrmion size decreases Bogdanov and Hubert (1994b); Wilson et al. (2014); Romming et al. (2015), and so do the internal energy barriers for both mechanisms, which we plot in Fig. 2a. In Fig. 2b, we show the change in configurational entropy upon reaching the SP, which is found to become less negative as the field increases. Since $\Delta S<0$ (entropic narrowing Desplat et al. (2018)), it is a stabilizing effect which lowers the attempt frequency. Lastly, we assume that the collapse processes are independent, so that the total rate of collapse is $\tau^{-1}_{\text{Langer}}(T)=k_{1}(T)+k_{2}(T)$.

Our aim is to compare the results of Langer’s theory with collapse rates obtained from forward flux sampling. The FFS method generates trajectories between two (meta)stable states $A$ and $B$ in a ratchet-like manner without imposing any bias on the microscopic dynamics, which makes it well adapted for the simulation of rare events. Compared to other path sampling methods, it does not require prior knowledge of the density of states, which makes it suitable for nonequilibrium systems that do not obey detailed balance. It employs a set of $n(+1)$ nonintersecting interfaces in phase space to sample the transition path ensemble and compute a transition rate. The interfaces $\{\lambda_{A},\lambda_{0}\dots\lambda_{n}=\lambda_{B}\}$ [Fig. 3] are defined as iso-surfaces of a monotonically varying order parameter, $\lambda(x_{i})=\lambda_{i}$, such that $x_{i+1}>x_{i}$ or $x_{i+1}<x_{i}$ for all $i$. Any trajectory going from $A$ to $B$ must cross all the interfaces at least once. The rate constant from $A$ to $B$ may be expressed as

[TABLE]

in which $\Phi_{A,0}$ is the rate at which trajectories starting from region $A$ cross the first interface $\lambda_{0}$, and the conditional probabilities $P(\lambda_{i+1}|\lambda_{i})$ correspond to the probability that a trajectory coming from $A$ that crossed $\lambda_{i}$ for the first time will cross $\lambda_{i+1}$ before returning to $A$. The protocol is illustrated in Fig. 3 and is as follows. First, a single Langevin simulation is started in state $A$ [Fig. 3a]. Each time the system successfully exits region $A$ and crosses $\lambda_{0}$, its configuration is stored. The simulation ends after $N_{0}$ crossing events have been recorded, and the flux of trajectories out of $A$ is obtained by $\Phi_{A,0}=N_{0}/\Delta t$, in which $\Delta t$ is the total simulated time. After that, a configuration stored at $\lambda_{0}$ is selected at random and used as a starting point for a new simulation [Fig. 3b]. That new simulation ends when the trajectory either crosses $\lambda_{1}$, in which case the crossing configuration is stored, or the system returns to $A$. This procedure is repeated $M_{0}$ times. If $N_{0}^{s}$ trajectories successfully crossed $\lambda_{1}$, we obtain $P(\lambda_{1}\mid\lambda_{0})=N_{0}^{s}/M_{0}$. One then proceeds analogously at subsequent interfaces. During the trial runs, Langevin simulations are carried out by integrating the system of stochastic Landau-Lifshitz-Gilbert equations at each timestep, by means of the stochastic Heun scheme García-Palacios and Lázaro (1998), for which details can be found in the SM sm_ . To obtain the rate of collapse of a skyrmion, a natural choice of order parameter is the skyrmion size, where state $A$ is the equilibrium skyrmion size, and state $B$ – corresponding to the FM state – is associated to a zero-size [Fig. 3]. Arbitrarily, we consider that magnetic sites $\mathbf{m}_{i}$ ($i=1\ldots N$) that satisfy $m_{z,i}\leq 0$ are part of the skyrmion, and we define the order-parameter $x$ as the (integer) number of sites inside the skyrmion. For values of the reduced field in $[0,0.05]$, we compute a total collapse rate, $k_{\text{FFS}}$. We give the results from FFS and Langer’s theory for the attempt frequency and the skyrmion lifetime against collapse in Figs. 4a and 4b. Through each FFS run, we set the temperature such that $\beta\Delta E_{1}=$ 10, so that Langer’s theory may hold Coffey et al. (2001). FFS runs are also carried out at $\beta\Delta E_{1}$ = 15 in the lower field region and yield very similar results, which shows that the attempt frequency has no significant $T$-dependence here. Since the translation of the skyrmion costs little energy compared to $k_{B}T$, we also show Langer’s result with a treatment of the translational modes as Goldstone modes Büttiker and Landauer (1981); Braun (1994); Bessarab et al. (2018). Details on the method can be found in the SM sm_ . We find that this treatment results in $f_{0}$ being overestimated in the low field region, and, in this system, the best overall agreement between Langer and FFS is obtained without considering Goldstone modes. In Fig. 5, we show some examples of stored configurations at interfaces at which the order-parameter equates that of a saddle point. We report both $S_{1}$- and $S_{2}$-types of configurations, as well as some other configurations that don’t clearly pertain to either category. This occurs because the crossing configurations correspond to an order parameter which is either equal to, or smaller than that of the SP, which does not imply that the configuration is in fact a SP. Under the effect of thermal fluctuations, the system does not usually cross the barrier exactly at the SP, but deviates from it by a more or less small amount. Lastly, at higher field values where $f_{0}$ is found within the GHz range, we also compute a collapse rate via direct Langevin simulations at $\beta_{300}\Delta E_{1}\approx$ 2 - 3. Following a similar procedure to Ref. Desplat et al. (2019), we compute an average lifetime out of 400 collapses. The results are shown in Fig. 4a for $b_{z}\geq 0.04$ and match the FFS results.

We have thus validated the use of Langer’s theory to obtain skyrmion collapse rates by means of forward flux sampling simulations. In the end, we find that due to a higher activation entropy, the path involving the third-order saddle point $S_{2}$ does not contribute significantly to the total rate. Nonetheless, the internal energy barriers associated with $S_{2}$ are almost the same as the ones associated with $S_{1}$, and configurations similar to $S_{2}$ are observed in FFS, so we could not justify neglecting it a priori. We should also note that since the difference in internal energy between $S_{1}$ and $S_{2}$ is quite small (around 0.08 $J_{\mathrm{ex}}$ at zero field), the mechanisms are probably not completely decoupled. Since the skyrmions are coupled to the lattice, we observe lattice effects in the Langer approach, that manifest in non-monotonic variations of the entropic contribution [Fig. 2b], which are in turn found in the attempt frequency and the average lifetime [Fig. 4]. In FFS simulations, lattice effects are likely smoothed out by thermal fluctuations. Nevertheless, FFS shows that translational modes should not be treated as Goldstone modes in this case. FFS is a valuable tool, as it requires no prior assumptions on the system. With appropriate interface design, it could be used to treat problems that have not yet been successfully solved by reaction rate theory, such as the problem of skyrmion nucleation rates. Most notably, FFS and Langer’s theory both yield variations of the collapse rate prefactor by three to four orders of magnitude when the applied magnetic field varies by 0.05 $J_{\mathrm{ex}}$. This effect is due to the important entropic contribution, and implies that reaching the 10-year retention rate necessary for technological applications may require adequately tuning the attempt frequency, in addition to the energy barrier. This result is valid for magnetic skyrmions, but also applies to any (meta)stable state undergoing a consequent change in entropy upon reaching the transition state. Here, a decrease in entropy at the saddle point stabilizes the skyrmion state. This is directly linked to the skyrmion’s internal modes Desplat et al. (2018). Since the skyrmion size decreases with the applied field, we find that the bigger the skyrmion, the higher $\Delta E$, the more negative $\Delta S$, the stronger the stabilizing effect (see also Ref. von Malottki et al. (2019)). On the other hand, for processes with high activation energies, the Meyer-Neldel compensation rule yields a destabilizing, often large entropic contribution Meyer and Neldel (1937); Yelon and Movaghar (1990); Yelon et al. (1992) (e.g. biological death rates, transport in semiconductors, decay of the skyrmion lattice, etc Rosenberg et al. (1971); Peacock-Lopez and Suhl (1982); Wild et al. (2017); Kamiya et al. (2010); Cooper et al. (2001)). These considerations underline the fact that, when estimating transition rates, one should not a priori assume a characteristic value of $f_{0}$, and special care needs to be taken in its evaluation.

Acknowledgements.

This work was supported by the Horizon 2020 Framework Programme of the European Commission, under Grant agreement No. 665095 (MAGicSky), and the Agence Nationale de la Recherche under Contract No. ANR-17-CE24-0025 (TOPSKY). The support from CD-Laboratory AMSEN (financed by the Austrian Federal Ministry of Economy, Family and Youth, the National Foundation for Research, Technology and Development) is acknowledged. FFS simulations were performed on the Vienna Scientific Cluster (VSC).