A magnetic domain wall Mackey-Glass oscillator J ér ôme Williame [START_REF] Grollier | Switching a spin valve back and forth by current-induced domain wall motion[END_REF] and Joo-Von Kim 1, a) Centre de Nanosciences et de Nanotechnologies, CNRS, Université Paris-Saclay, 91120 Palaiseau, France (Dated: 19 April 2021) We propose a time-delay oscillator with Mackey-Glass nonlinearity based on a pinned magnetic domain wall in a thin film nanostrip. Through spin transfer torques, electric currents applied along the strip cause the domain wall to deform and displace away from a geometrical pinning site, which can be converted into a nonlinear transfer function through a suitable choice of a readout. This readout serves as a delay signal, which is subsequently fed back into the applied current with amplification. With micromagnetics simulations, we study the role of the readout position, time delay, and feedback gain on the dynamics of this domain wall. In particular, we highlight regimes in which self-sustained oscillations and complex transients are possible.

Magnetic domain walls in thin films possess a number of useful properties for information technologies. Domain walls are topological solitons whose spatial profile is determined by competing exchange interactions and anisotropies, which confer them particle-like properties that make them useful for binary storage or logic gates. For example, the possibility to move sequences of domain walls back and forth under applied currents [START_REF] Grollier | Switching a spin valve back and forth by current-induced domain wall motion[END_REF][START_REF] Vernier | Domain wall propagation in magnetic nanowires by spin-polarized current injection[END_REF][START_REF] Yamaguchi | Real-Space Observation of Current-Driven Domain Wall Motion in Submicron Magnetic Wires[END_REF][START_REF] Yamanouchi | Current-induced domain-wall switching in a ferromagnetic semiconductor structure[END_REF][START_REF] Kläui | Controlled and Reproducible Domain Wall Displacement by Current Pulses Injected into Ferromagnetic Ring Structures[END_REF] allows for information to be shuttled back and forth in magnetic shift registers [START_REF] Parkin | Magnetic Domain-Wall Racetrack Memory[END_REF] . Domain wall devices have been touted as a promising avenue to reduce the energy consumption in information technologies.

Beyond applications in storage and logic, attention in spintronics has recently turned towards other forms of non-Boolean information processing such as neuro-inspired computing. Such efforts are motivated by the goal of identifying ways to mimic the capacity of the human brain to efficiently perform cognitive tasks, such as pattern recognition, with noisy or incomplete data. The inherent nonlinearity in magnetization dynamics is useful in this respect, since neuronal dynamics exhibit complex transients, self-oscillations, and chaotic behavior. For example, it has been shown recently that spin-torque nano-oscillators can be used for classification tasks [START_REF] Torrejon | Neuromorphic computing with nanoscale spintronic oscillators[END_REF][START_REF] Romera | Vowel recognition with four coupled spin-torque nanooscillators[END_REF][START_REF] Riou | Temporal Pattern Recognition with Delayed-Feedback Spin-Torque Nano-Oscillators[END_REF] .

Here, we examine a feature of neuronal dynamics that has been explored sporadically in spintronics, namely that of magnetization dynamics with delayed feedback effects. In general terms, this entails feeding the output of the dynamical system, after amplification and a time delay, back into its input as an additional component of the external drive (e.g., fields or currents). Delayed feedback has been explored in terms of selfinjection in spin-torque nano-oscillators, where their sensitivity to external perturbations [START_REF] Tiberkevich | Sensitivity to external signals and synchronization properties of a non-isochronous auto-oscillator with delayed feedback[END_REF] and changes to characteristics such as the spectral linewidth have been examined [START_REF] Khalsa | Critical current and linewidth reduction in spin-torque nano-oscillators by delayed self-injection[END_REF][START_REF] Tamaru | Analysis of phase noise in a spin torque oscillator stabilized by phase locked loop[END_REF][START_REF] Tsunegi | Self-Injection Locking of a Vortex Spin Torque Oscillator by Delayed Feedback[END_REF][START_REF] Singh | In-teger, Fractional, and Sideband Injection Locking of a Spintronic Feedback Nano-Oscillator to a Microwave Signal[END_REF][START_REF] Singh | Self-Injection Locking of a Spin Torque Nano-Oscillator to Magnetic Field Feedback[END_REF][START_REF] Williame | Effects of delayed feedback on the power spectrum of spin-torque nano-oscillators[END_REF] . Feedback in low-dimensional systems such as macrospin models can result, for example, in nontrivial behavior such as temporal chaos [START_REF] Williame | Chaotic dynamics in a macrospin spin-torque nano-oscillator with delayed feedback[END_REF][START_REF] Taniguchi | Chaos in nanomagnet via feedback current[END_REF] . In other systems, the onset of self-sustained oscillations can appear through a variety of different feedback mechanisms [START_REF] Dixit | Spintronic oscillator based on magnetic field feedback[END_REF][START_REF] Kumar | Coherent microwave generation by spintronic feedback oscillator[END_REF][START_REF] Bhuktare | Spintronic Oscillator Based on Spin-Current Feedback Using the Spin Hall Effect[END_REF] . More broadly, delayed-feedback has been employed in nonlinear elements as building blocks of time-delay architecture for reservoir computing [START_REF] Appeltant | Information processing using a single dynamical node as complex system[END_REF][START_REF] Larger | High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification[END_REF] , a neuroa) Electronic mail: joo-von.kim@c2n.upsaclay.fr inspired scheme that has garnered growing interest in magnetism and spintronics [START_REF] Riou | Temporal Pattern Recognition with Delayed-Feedback Spin-Torque Nano-Oscillators[END_REF][START_REF] Bourianoff | Potential implementation of reservoir computing models based on magnetic skyrmions[END_REF][START_REF] Prychynenko | Magnetic Skyrmion as a Nonlinear Resistive Element: A Potential Building Block for Reservoir Computing[END_REF][START_REF] Nakane | Reservoir Computing With Spin Waves Excited in a Garnet Film[END_REF][START_REF] Marković | Reservoir computing with the frequency, phase, and amplitude of spin-torque nanooscillators[END_REF][START_REF] Araujo | Role of non-linear data processing on speech recognition task in the framework of reservoir computing[END_REF][START_REF] Yamaguchi | Periodic structure of memory function in spintronics reservoir with feedback current[END_REF][START_REF] Pinna | Reservoir Computing with Random Skyrmion Textures[END_REF] .

A prototypical example of delayed feedback dynamics is the Mackey-Glass oscillator, which is described by the onedimensional delay differential equation [START_REF] Mackey | Oscillation and chaos in physiological control systems[END_REF] ,

ẋ(t) = -Γx(t) + Cx(t -τ) 1 + x(t -τ) p . (1)
The right-hand side of the equation comprises two terms. The first represents damping with a relaxation rate Γ. The second is the delayed feedback term, where C is an amplification, and p determines the nonlinearity of the transfer function. Depending on the values of Γ, C, and p, the system in Eq. ( 1) can exhibit complex transients, self-oscillations, and chaos. In the absence of delayed feedback, Eq. ( 1) exhibits only trivial relaxation dynamics. Optoelectronic schemes exploiting the kind of Mackey-Glass nonlinearity in Eq. ( 1) have been recently used in time-delay architectures for reservoir computing [START_REF] Appeltant | Information processing using a single dynamical node as complex system[END_REF][START_REF] Larger | High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification[END_REF] . In the context of neuro-inspired information processing, it is, therefore, interesting to enquire how similar implementations could be realized using magnetic states and spintronic systems.

Our proposal for a domain wall Mackey-Glass oscillator (DWMGO) is presented in Fig. 1. The system comprises an in-plane magnetized ferromagnetic strip with an elliptical protrusion, which acts as a strong pinning site for a transverse domain wall, similar to a magnetic racetrack memory [START_REF] Parkin | Magnetic Domain-Wall Racetrack Memory[END_REF] . An electric current is applied along x, the long axis of the strip, which simultaneously leads to a displacement and a distortion of the domain wall along this direction. A magnetoresistive readout sensor is positioned next to the protrusion, which is sensitive to the y component of the magnetization. The signal from this readout, which represents the spatial average across the sensor at position i, m y (t) i , is amplified by a factor G with a time delay τ and subsequently fed back into the input as an additional time-dependent current, δJ(t). The current-driven domain wall dynamics is then subjected to the time-dependent drive,

J(t) = J 0 + Ga i m y (t -τ) i , (2)
where a i is a normalization constant discussed below. The dynamics of the DWMGO is simulated using the Mu-Max3 code 32 , which performs a numerical time integration arXiv:2104.00080v2 [cond-mat.mes-hall] 16 Apr 2021

dm dt = -γ 0 m×H eff +αm× dm dt -(u • ∇) m+βm×(u • ∇) m. ( 3
)
γ 0 is the gyromagnetic ratio, m(r, t) is a unit vector representing the magnetization field, H eff is the effective magnetic field, and α is the Gilbert damping constant. The Zhang-Li spin transfer torques are parametrized by u = J(r)Pµ B /(eM s ), which represents an effective spin drift velocity where J is the conventional current density, µ B is the Bohr magneton, e is the electron charge, β is the nonadiabaticity, and P is the spin polarization. The effective field is given by the variational derivative of the total magnetic energy U with respect to the magnetization unit vector, H eff = -(1/µ 0 M s )δU/δm, and comprises contributions from the exchange, dipole-dipole, and the Zeeman interactions.

The oscillator studied comprises a strip with dimensions of 3072 × 200 × 10 nm, which is discretized using finite difference cells 3 × 3 × 10 nm in size. The protrusion is taken to be a half ellipse with the long axis perpendicular to the wire, which extends out to a distance of 92 nm from the wire, and the short axis parallel to the wire with an extension of 100 nm. We assume micromagnetic parameters relevant for permalloy, namely a saturation magnetization of M s = 860 kA/m, an exchange constant of A = 13 pJ/m, and a Gilbert damping constant of α = 0.02. The readout region is taken to be rectangular in shape with dimensions of 100×200 nm, which spans the entire with of the magnetic strip. For the current-driven spin torques, we assume P = 1 for simplicity and a nonadiabacity of β = 0.1. Finite element electromagnetics simulations were used to compute the spatial profile of the currents flowing along the wire. Current densities used here refer to the average current density applied through the rectangular cross section at the ends of the wire, far from the protrusion. In Fig. 2, the static domain wall profile under different applied dc currents is shown.

The magnetoresistive signal, which serves as the transfer parison between the different sensor positions, we rescale the transfer functions in the following way. First, we normalize each curve by a factor a i such that the maximum value of the output is unity; we define the transfer function to be f DW = a i m y i . Second, we define an operating point J * 0,i such that the value of the transfer function is 0.5 at J 0 -J * 0 = 0. The values of a and J 0 for each sensor position is given in Table I. The transfer function for positions 1 to 3 is shown in Fig. 3(b). For the purposes of comparison, the transfer function for the Mackey-Glass oscillator [Eq. ( 1)] is shown in Fig. 3(c) for three values of the nonlinearity, p. f DW captures the nonmonotonic behavior of the Mackey-Glass transfer function for p > 1, with the sensor position offering a means to adapt the nonlinearity.

Next, we discuss the possible dynamical states of the DWMGO and the influence of the sensor position. Simulations with the delayed feedback are performed as follows. For a given value of τ under study, time integration of Eq. ( 3) is first performed over the interval t ∈ [-τ, 0] without the feedback term in order to generate the delay term m y (t-τ) , which is stored in memory. Since we employ a time integration algorithm with an adaptive time step (i.e., primarily the Dormand-Prince method implemented in MuMax3 32 ), we use cubic interpolation for the function m y (t -τ) so that it can be evaluated at arbitrary values of t -τ during the simulation. During this initial phase without feedback, a static current J 0 = J * 0 is applied. The initial magnetization state is computed previously for this operating point in order to minimize transient dynamics. For t > 0, Eq. ( 3) is solved together self-consistently with Eq. (2). Some examples of the domain wall dynamics are shown in Fig. 4. The figure shows the simulated magnetoresistance readout for a sensor in position 2. It is instructive to first consider the effect of long delays (τ = 10 ns). At low gain, G = 2, the feedback acts like a driving term formed by a sequence of steps, where each step corresponds to a different equilibrium position of the domain wall [Fig. 4(a)]. As the step is initially applied, a ringdown of the signal can be seen over a few ns, which corresponds to the damped oscillatory motion of the domain wall about the new equilibrium state. Since f DW is designed to be finite at J 0 = J * 0 , the application of the feedback guarantees some transient behavior since J will deviate from J * 0 instantaneously. The stepped ringdown response eventually converges to a steady state value. At higher gain, G = 6, the initial stepped response comprises larger ampli- tude motion, which over the course of a few hundred ns results in a steady state oscillation of the domain wall, as shown in Fig. 4(b). The initial phase is dominated by transients on a time scale give by τ, but this gradually evolves into a selfoscillatory mode for which this modulation is less present. At short delays (τ = 1 ns), the transient dynamics toward the selfoscillation is considerably shorter, which can be observed in Fig. 4(c).

Examples of the resulting power spectrum of feedbackdriven oscillations are presented in Fig. 5. The spectra in Fig. 5(a) and 5(b) correspond to the oscillations shown in Fig. 4(c) and 4(b), respectively. Fig. 5(a) presents an example at short delays (τ = 1 ns), where we can clearly identify the main oscillation peak f 0 at around 1.13 GHz, with harmonics at 2 f 0 and 3 f 0 also visible in the spectrum shown. These main spectral lines are accompanied by a number of modulation sidebands, which result from a self-modulation that is induced by the delayed feedback signal. This is a feature seen in other delayed-feedback systems involving spintronic devices [START_REF] Singh | In-teger, Fractional, and Sideband Injection Locking of a Spintronic Feedback Nano-Oscillator to a Microwave Signal[END_REF][START_REF] Singh | Self-Injection Locking of a Spin Torque Nano-Oscillator to Magnetic Field Feedback[END_REF][START_REF] Williame | Effects of delayed feedback on the power spectrum of spin-torque nano-oscillators[END_REF] . Figure 6 presents a color map of the power spectral density of the DWMGO as a function of time delay at a fixed gain of G = 6, for three sensor positions. For all cases shown, we can observe the primary oscillation peak around 1 GHz, where the strength and number of modulation sidebands vary with the delay. For short delays, the feedback results only in a modulation of the primary oscillation frequency, where a ratchet-like frequency pulling toward lower frequencies is shown as τ increases. This phenomenon is particularly visible in Fig. 6(c). For longer delays, this frequency pulling is accompanied by the appearance of modulation sidebands, which are also visible in the harmonics within the frequency window considered. The onset of these sidebands depend strongly on the sensor position, and therefore, on the feedback nonlinearity. We can see that the sidebands appear at τ 1.5 ns for a sensor centered at the protrusion (position 1), while at positions 2 and 3 the onset of the sidebands occur at τ 4.0 and τ 4.5 ns, respectively. These results show that the rich power spectra can be obtained within a single device by tuning the readout position of the domain wall oscillations.

Figure 7 presents the output power of the DWMGO for different values of the delay time and amplification with the three sensor positions considered in Fig. 6. Similarly to behavior in Fig. 6, we can observe modulations in the output power as the delay is varied. Another notable feature is that a maximum in the oscillator power is observed over a range of feedback gain, which can be clearly seen in Fig. 7(c). This reflects the fact that the oscillatory properties are strongly dependent on the competition between the relaxation dynamics of the confining potential and the form of the feedback signal.

We now discuss two possible issues and their potential impact on the DWMGO functionality described so far. The first concerns the Joule heating related to the current flow in the device. Because the current densities involved are typical for spin-transfer torques, i.e., in the range of 10 12 A/m 2 , we can expect the operating temperature of the device to be significantly above ambient conditions. However, since the crosssectional area at the protrusion is larger than the straight portion of the wire, we do not expect additional complications due to Joule heating near the protrusion or the sensor region, where heating might actually be lower than the baseline value set by the regions far from the protrusion. In terms of thermal noise in the readout, we do not expect any additional issues beyond those faced in designing suitable magnetoresistive readouts, e.g., in three-terminal devices studied in spinorbit torque switching. The second issue concerns the accuracy of the sensor position, where the vagaries of nanofabrication can also result in variations in its size and shape. Our simulation results using different sensor positions show that the qualitative features of the oscillator remain present irrespective of the sensor position. In a given experimental device, there remain other control parameters such as the feedback delay, feedback gain, and the dc current used as the operating point to tailor the desired dynamics. Shifts in the sensor position away from the position, as shown in Fig. 3, result in an overall translation of the nonlinear transfer function along the dc current axis [Fig. 3(b)], so variations in the sensor position due to fabrication can be accommodated by tuning this operating point.

In summary, we have presented a model for a spintronic feedback oscillator with a Mackey-Glass nonlinearity. In contrast to conventional spin-torque nano-oscillators in which nonlinearities are determined by intrinsic micromagnetic energies such as shape anisotropies or spin wave interactions, here the form of the nonlinearity is determined in large part by the position of the readout sensor and on how the pinned domain wall deforms as a function of the applied current. This suggests that multiple nonlinearities, and possibly different functionalities, could be designed on a single device by using different readout sensors.

FIG. 1 .

 1 FIG.1. Geometry of the domain wall Mackey-Glass oscillator. The magnetic strip, 200 nm in width, hosts an in-plane magnetized transverse domain wall, which is pinned by an elliptical protrusion of width 100 nm. A dc conventional current J 0 flows along the track and displaces the wall through spin-transfer torques. A magnetoresistance readout sensor detects the m y component near the protrusion, which is amplified by a factor G and fed back into the input current as a delayed modulation.

FIG. 2 .FIG. 3 .

 23 FIG. 2. Static domain wall profiles under different applied dc currents, J(t) = J 0 . (a) J 0 = -4, (b) J 0 = 0, (c) J 0 = 4, (d) J 0 = 8, and (e) J 0 = 10 TA/m 2 . The critical depinning current is approximately J 0,c = 8.5 TA/m 2 . The functionality of the DWMGO is limited to currents below this critical current.

FIG. 4 .

 4 FIG. 4. Time traces of the magnetization dynamics measured at Position 2. (a) Transient dynamics toward a static configuration under weak amplification, G = 2, and long delays, τ = 10 ns. (b) Selfsustained oscillations at G = 6 and long delays, τ = 10 ns. Selfsustained oscillations at G = 6 and short delays, τ = 1 ns.

FIG. 6 .

 6 FIG. 6. Color map of the power spectral density (PSD) as a function of delay τ for G = 6, with the sensor at (a) position 1, (b) position 2, (c) position 3.

FIG. 7 .

 7 FIG. 7. Phase diagram of the output power of the oscillator as a function of feedback delay, τ, and gain, G, for the readout sensor located at (a) position 1, (b) position 2, and (c) position 3.

TABLE I .

 I Scaling parameter a i and operating point J * 0,i for each sensor position i such that f DW = 0.5 at J 0 = J * 0,i .

	Position, i	a i	J * 0,i (TA/m 2 )
	1	1.145	-4.12
	2	1.201	-0.68
	3	1.307	0.42
	4	1.366	2.20
	5	1.413	4.38
	6	1.366	6.24

ACKNOWLEDGMENTS

This work was supported by the Agence Nationale de la Recherche (France) under contracts no. ANR-14-CE26-0021 (MEMOS) and ANR-17-CE24-0008 (CHIPMuNCS).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.


L Desplat
email: louise.desplat@ipcms.unistra.fr

J.-V Kim

Quantifying the thermal stability in perpendicularly magnetized ferromagnetic nanodisks with forward flux sampling

The thermal stability in nanostructured magnetic systems is an important issue for applications in information storage. From a theoretical and simulation perspective, an accurate prediction of thermally-activated transitions is a challenging problem because desired retention times are on the order of 10 years, while the characteristic time scale for precessional magnetization dynamics is of the order of nanoseconds. Here, we present a theoretical study of the thermal stability of magnetic elements in the form of perpendicularly-magnetized ferromagnetic disks using the forward flux sampling method, which is useful for simulating rare events. We demonstrate how rates of thermally-activated switching between the two uniformly-magnetized "up" and "down" states, which occurs through domain wall nucleation and propagation, vary with the interfacial Dzyaloshinskii-Moriya interaction, which affect the energy barrier separating these states. Moreover, we find that the average lifetimes differ by several orders of magnitude from estimates based on the commonly assumed value of 1 GHz for the attempt frequency.

I. INTRODUCTION

The thermal stability of magnetic states is a challenging problem that underpins the utility of magnetism for information storage. It involves understanding the average retention (or dwell) time of a magnetic bit, which typically comprises regions of uniform magnetization in hard disk media or uniform states in magnetoresistive random access memories. This retention time, τ , can be obtained as the inverse of the transition rate, τ -1 = k, which is governed by an Arrhenius relation of the form [1]

k = f 0 exp - ∆E k B T , (1)
where f 0 is the Arrhenius prefactor often referred to as the "attempt frequency", ∆E is the energy barrier separating the binary '0' ("up") and '1' ("down") states, k B is Boltzmann's constant, and T is the temperature. The vast majority of studies to date have focused on understanding and optimizing the barrier ∆E, whilst assuming a nominal value of f 0 1 GHz which captures the typical time scales of damped precessional magnetization in strong ferromagnets. Based on this assumption, the typical metric of a 10-year retention time demands that ∆E/k B T 50 for operation at room temperature, a rule-of-thumb that has provided guidance both theoretically and experimentally for the feasibility of using various magnetic states for information storage [2][3][4][5][6][7].

One example of current interest in which the issues of thermal stability are intertwined with the complexity of nonuniform magnetic states concerns ferromagnetic nanostructures with large perpendicular magnetic anisotropy (PMA). Such systems are attractive because they offer larger storage densities than in-plane magnetized systems [8]. However, the reversal process in such systems can be nonuniform [5,7,[9][10][11][12][13][14][15][16]. For nm-thick ferromagnetic films with lateral dimensions in the tens to hundreds of nm, magnetization reversal takes place through the nucleation and propagation of magnetic domain walls, whose energies then govern the energy barrier required to transition from one metastable state to the other [9][10][11]. Moreover, PMA often involves coupling to strong spin-orbit materials, which can also induce an antisymmetric exchange in the form of a Dzyaloshinskii-Moriya interaction (DMI) [17,18]. Indeed, the presence of DMI has been shown to be detrimental to the thermal stability in such structures [5,7,12], which results from the fact that the DMI reduces the domain wall energy [19] and therefore the energy barrier.

Recent studies on magnetic skyrmions, however, have shown that the barrier alone does not govern their thermal stability [20][21][22][23][24][25]. It has been found that a strong entropic contribution to the prefactor can arise as a result of the complex topology of the energy surface resulting from the large number of degrees of freedom present. This can be understood from a generalization of Kramers' transition rate theory [1,26] to multidimensional phase spaces, as developed by Langer [27], in which the activation entropy is given by the spectrum of small fluctuations about the initial (meta)stable state and the transition state. Since the DMI also has a strong influence on the spin wave spectrum in the nominally uniformlymagnetized [28] and domain wall states [29,30], the degree to which entropic contributions influence the thermal stability of PMA memory elements with domain-wall mediated reversal remains to be explored.

Here, we revisit the question of thermal stability in PMA disks using the method of forward flux sampling (FFS) [25,[31][32][33][34][35][36]. This method was developed for tackling the problem of rare events, which are unlikely to appear during the course of a conventional simulation run, due to the fact that the mean waiting time between events is much larger than the timescale of the dynamics. In magnetism, FFS has been applied to a similar problem of thermal stability in graded media, where its efficacy with respect to brute-force Langevin simulations was clearly demonstrated [35]. The method has also been applied recently to the study of skyrmion lifetimes, where agreement was found with another approach based on Langer's theory [25].

The remainder of the article is organized as follows. In Section II, we present the system studied and the implementation of the forward flux sampling method. Section III discusses the application of the method to determine the change in average dwell time with the Dzyaloshinskii-Moriya interaction. A discussion and concluding remarks are given in Section IV.

II. GEOMETRY AND METHOD

The system studied comprises a perpendicularlymagnetized ferromagnetic disk, which simulates the free magnetic layer in a magnetic memory device. Following Ref. 5, we consider a CoFeB film with a saturation magnetization of M s = 1.03 MA/m, an exchange constant of A = 10 pJ/m, and a perpendicular anisotropy constant of K u = 0.77 MJ/m 3 . We also take into account the presence of an interfacial Dzyaloshinskii-Moriya interaction, D, whose strength is varied. We model a 1-nm thick disk with a diameter of 32 nm, which is discretized with 32 × 32 × 1 finite difference cells. Let z denote the axis perpendicular to the film, which is defined by the xy-plane.

We use the micromagnetic approximation in which we consider the Langevin dynamics of the magnetization vector, m(r, t) ≡ M(r, t)/M s , such that m(r, t) = 1. The Langevin dynamics is obtained by stochastic time integration of the Landau-Lifshitz equation with a fluctuating thermal field

dm dt = -γ 0 m × (H eff + h th ) + αm × dm dt , ( 2
)
where γ is the gyromagnetic constant, γ 0 = µ 0 γ, and α is the Gilbert damping constant. The dynamics is governed by the deterministic effective field, H eff = -(µ 0 M s ) -1 δU/δm, which is obtained from the variational derivative of the total micromagnetic energy, U , with respect to the magnetization, and a stochastic field, h th , which takes into account finite temperature effects. This thermal field has zero mean, h th = 0, and represents a Gaussian white noise with the spectral properties

h th,i (r, t)h th,j (r , t ) = 2αk B T µ 0 V δ ij δ(r -r )δ(t -t ), (3)
where i, j represent the different Cartesian components of the field vector. We used two methods to compute this stochastic time integration: a homemade code that employs a Heun scheme [23,25,37], and the MuMax3 code [38] with an adaptive time step scheme [39]. Let 'A' denote the initial magnetic state, which comprises a uniformly-magnetized state along the +z direction, and 'B' the second degenerate metastable state, which is uniformly magnetized along -z and separated from 'A' by an energy barrier ∆E. Our goal is to compute the average lifetime τ of the state 'A' at finite temperatures. Estimating τ directly through bruteforce Langevin-dynamics simulations would involve starting with 'A' as the initial configuration and integrating Eq. ( 2) until the state 'B' is reached. This process would then require to be repeated a few hundred times to obtain reasonable statistics. Because of the precessional dynamics, typical time steps for the numerical time integration are in the range of 1 to 100 fs, so it does not appear fruitful to proceed with this program of work where τ is of the order of years for technologically-relevant systems.

Path sampling methods such as forward flux sampling (FFS) [31][32][33][34] can be used instead to estimate the transition rates of rare events such as the thermally-activated escape from 'A' to 'B' [25,35,36]. The basic idea of the FFS approach is illustrated in Fig. 1.

The method involves generating trajectories between 'A' and 'B' in a ratchetlike manner through a series of interfaces λ(ζ) in phase space, {λ A , λ 0 , λ 1 , . . . , λ n-1 , λ n = λ B }, which represent non-intersecting isosurfaces of a monotonically varying order parameter ζ. The basin 'A' is defined for ζ > ζ A , similarly 'B' is defined for ζ < ζ B . For the problem studied here, the order parameter is chosen to be the spatial average of the m z component of the magnetization, ζ ≡ m z (r) = (V ) -1 dV m z (r), and decreases monotonically between 'A' and 'B'. A choice of order parameter that is close to the reaction coordinate -i.e., the minimum energy path (MEP) of highest statistical weight through phase space -improves the computational efficiency of the method. All trajectories from 'A' to 'B' must traverse each interface at least once and the overall transition rate is expressed as

k AB = Φ A,0 P B , (4)
where Φ A,0 represents the flux of trajectories from 'A' to the first interface λ 0 , and P B ≡ P (λ B |λ 0 ) is the probability that a trajectory that crossed λ 0 coming from 'A' reaches λ B before returning to the basin 'A'. The quantity Φ A,0 can be obtained in a straightforward manner since the trajectories emanating from 'A' cross the λ 0 interface with high frequency as a result of its proximity. On the other hand, the probability P B will typically be very small for rare events. However, calculating this probability becomes manageable by decomposing it into a product of partial fluxes at each interface,

P B ≡ P (λ n |λ 0 ) = n-1 i=0 P (λ i+1 |λ i ), (5)
where the conditional probability P (λ i+1 |λ i ) represents the probability that a trajectory starting at λ i reaches λ i+1 before returning to the basin 'A'. Calculation of the FFS rate k AB proceeds in two steps. The first involves initiating the system in 'A' and performing a Langevin dynamics simulation by stochastic time integration of Eq. ( 2). This simulation is used to collect micromagnetic configurations at λ 0 that result from instances in which the system commences in 'A' and crosses the interface λ 0 , as shown in Fig. 1(b). One then waits for the system to return to 'A' before the next crossing configuration at λ 0 can be saved. The simulation proceeds until N 0 such configurations have been obtained, at which the flux of trajectories out of 'A' that cross λ 0 can be estimated as Φ A,0 = N 0 /∆t sim , where ∆t sim is the total simulated time required to obtain the N 0 crossings. Note that ∆t sim should not include the time the system may have spent in the 'B' basin. The second step involves computing the probability P B . This begins by selecting at random a stored configuration at λ 0 and performing a Langevin dynamics simulation until either the system returns to 'A', which counts as a failed crossing, or it reaches the next interface λ 1 , in which case the micromagnetic configuration at λ 1 is stored. This process is repeated M 0 times, and the probability of reaching λ 1 from λ 0 is computed as P (λ 1 |λ 0 ) = N s 0 /M 0 , where N s 0 denotes the number of successful crossings of λ 1 . The procedure then continues in an analogous manner for the subsequent interfaces until λ B is reached, as shown in Fig. 1(c). The overall FFS simulation is successful if at least one trajectory reached λ B .

III. VARIATION OF THERMAL STABILITY WITH DMI

A. Underdamped limit with full dipolar interactions

We first present results using the MuMax3 code with full dipolar interactions, a realistic damping value of α = 0.01 [40], a variable DMI constant D between 0 and 2 mJ/m 2 , and T = 300 K. For this set of simulations, we used N 0 = 50 and M 0 = 1000 with 16 interfaces (n = 15). The interfaces were constructed as follows. The boundary of the basin 'A', λ A , was determined by the median of the value of m z,0 = m z of the thermally equilibrated state, which leads to 'A' being defined as ζ > m z,0 . By symmetry, we define λ B at ζ B = -m z,0 , with 'B' occupying the region of phase space ζ < -m z,0 . We note that m z,0 is a temperature-dependent quantity and also varies strongly as a function of the DMI, since the boundary conditions at the disk edges result in a canting of the magnetic moments away from the z axis [28,41].

Instead of spacing the interfaces equally from ζ = m z,0 to ζ = -m z,0 , we chose instead a weighting function based on the tanh(ζ) function that better mimics the reversal path. For example, for a magnetization profile of the form m z (χ) = tanh(χ), where χ is the reaction coordinate [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF], the choice of interface would represent regular spacing along the χ axis, rather than m z . Ultimately, this is a matter of convenience as the FFS method is not overly reliant on the particular placement of interfaces, which should in principle only affect the efficiency of the method [33], as we show further below.

In Fig. 2, we present the averaged magnetization configurations at different interfaces for three values of D. The configurations represent the stored state at an interface λ i following a successful traversal from λ i-1 . Since wall nucleation is not restricted to a particular edge of the disk, there are a multitude of paths in which the wall can traverse the disk during the reversal. In order to obtain a meaningful average, we first rotate each configuration such that wall displacement takes place from left to right during the reversal, in line with the schematic presented in Fig. 1(a). The overall behavior is similar for all values of D considered and corresponds to the minimum energy path predicted in earlier work [5,10]. The first interface involves the apparition of a nucleation zone at the disk boundary, which is most pronounced for D = 2 mJ/m 2 by virtue of the greater tilt of the magnetization at the disk edges. This proceeds with the propagation of a domain wall that sweeps through the disk, which is also symmetric about the λ 7 interface at which ζ = 0. The domain wall exhibits a stronger curvature at the disk center for the largest value of the DMI considered, because the DMI lowers the energy of the curved wall with respect to a straight wall, and selects a preferred direction of the curvature, in the same way that it selects right-or left-rotating Néel skyrmions depending on the sign of D. The configurations in Fig. 2 show that the trajectories generated by Langevin dynamics resemble the minimum energy path obtained with path finding schemes such as the nudged elastic band method [5,[START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF][START_REF] Jónsson | Classical and Quantum Dynamics in Condensed Phase Simulations[END_REF][START_REF] Bessarab | [END_REF] or the string method [12,45], which typically constitutes the path of largest statistical weight. We note that in FFS simulations, the thermal fluctuations at 300 K consistently excite the wall's curving modes at the saddle point for large DMI, while in MEPs the relaxed wall in the center tends to remain straight for all values of D.

The properties of the conditional probabilities p i-1 ≡ P (λ i |λ i-1 ) are presented in Fig. 3. The variation of p i-1 with the order parameter ζ i is shown in Fig. 3(a). For the first set of interfaces in the vicinity of 'A', the conditional probabilities lie are around 10%, and then progressively increase as the trajectories propagate toward λ B = λ(ζ B ). For D = 0 mJ/m 2 , there is a sharp increase after ζ 0.4 (λ 5 ) and saturation is attained after ζ -0.4 (λ 11 ). The variation is more gradual for the cases with finite values of D, where interestingly for D = 2 mJ/m 2 saturation toward unity does not occur. This suggests that numerous recrossings of the barrier take place even when the state approaches the vicinity of 'B'. This is a general feature of all the cases studied, albeit to different degrees, which is evidenced by the fact that p i-1 approaches 1 well beyond the interface where ζ = 0, i.e., at the top of the energy barrier. In other words, the underdamped dynamics of the magnetization precession means that reversal is not guaranteed even if the energy barrier is crossed [37]. This can also be seen in the evolution of the cumulative product P (λ i |λ 0 ) with ζ, shown in Fig. 3(b), which describes the probability of reaching λ i given the starting point of λ 0 . This function is constructed from the conditional probabilities in Fig. 3(a) as

P (λ i |λ 0 ) = i-1 j=0 P (λ j+1 |λ j ).
In this representation, we observe that most of the evolution in P (λ i |λ 0 ) occurs for ζ > 0, i.e., as the energy barrier is surmounted from 'A', with a further reduction by a factor of ∼ 5 occurring due to barrier recrossings for ζ < 0.

B. Overdamped limit with effective perpendicular anisotropy

The data presented in Figs. 2 and 3 represent a single FFS run. While the computational time required to execute this task was considerably shorter than full Langevin dynamics simulations of the entire reversal process, it was not feasible to obtain statistics of the interface probabilities, i.e., different realisations of the sampling for a given value of D, within a reasonable timeframe. For instance, each data point for the conditional probability in Fig. 3(a) required 1 to 3 days of simulation time on a single NVIDIA GTX 1080 graphics processor unit. Part of this difficulty stems from the calculation of the long-range dipolar interactions, which is computationally intensive. Another difficulty involves the underdamped nature of the precessional dynamics, which results in multiple recrossings at any given interface.

To ascertain whether the trends observed in Fig. 3 are representative of the reversal process, we conducted a different set of simulations in which the long-ranged dipolar interactions are approximated by a rescaling of the perpendicular anisotropy constant, i.e.,

K eff = K u - 1 2 (N z -N x )µ 0 M 2 s = 187 kJ/m 3
, where N x = 0.0418 and N z = 0.916 are demagnetizing factors of the disk [46]. The use of a local approximation for the dipolar interaction is less computationally intense and speeds up the simulation time. We first considered the overdamped limit, α = 0.5, which facilitates comparisons with Langer's theory [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF]. The cumulative probability function P (λ i |λ 0 ) for different values of the DMI is shown in Fig. 4.

The corresponding conditional probabilities P (λ i |λ i-1 ) are shown in the inset of each subfigure. For each value of D, between 7 and 12 distinct FFS runs were performed and are represented by the different colored curves, while the ensemble-averaged result is given by the black curve. While there is some spread in the final probability P (λ B |λ 0 ), the curves for the cumulative probability and the conditional probability share the same qualitative trend for a given value of D. The manner in which the shape of the P (λ i |λ i-1 ) curves change as D is increased also mirrors the evolution seen for the case presented in Fig. 3(b), where the 'S'-shaped variation at D = 0 transitions to more of a square-root like behavior at D = 2 mJ/m 2 .

C. Details of the optimization in the underdamped limit with effective perpendicular anisotropy a. Optimization In the underdamped limit of α = 0.01, similar to Section III A, recrossings are more frequent and the system deviates more significantly from the MEP, resulting in longer waiting times between magnetization reversals compared to the overdamped case. Longer waiting times typically imply longer CPU times. Nevertheless, the FFS simulation times can be reduced through several methods.

First, many trial runs at a given interface can be performed simultaneously on one or multiple CPU units. The only requirement is that the simulation waits for all trials run at one interface to terminate before moving on to the next interface. This can be managed via MPI on multiple nodes, or a simple Bash script on a single node.

Second, the parameters of the FFS run can be optimized in order to improve the efficiency of the method, as derived by Borrero and Escobedo [34]. The variance of the rate yielded by the trial runs, i.e., V [P B ] [32], is minimized under the constraint that the total rate k AB must remain constant. The best efficiency for FFS is obtained for a minimal relative variance of the rate, which is yielded by a constant flux of partial trajectories through all interfaces, i.e., M i p i = cst. We define the relative variance of the rate as [32],

V = n-1 i=0 q i p i k i , (6)
where we have used the simplified notations p i ≡ P (λ i+1 |λ i ), q i ≡ 1 -p i , and k i ≡ M i /N 0 . One may therefore either optimize the number of trial runs per interface, {M i }, at fixed interface placement, or optimize the placement of the interfaces, {ζ i }, at fixed {M i }. In the latter case, the interface placement is used to optimize the conditional probabilities, {p i }. In this work, we chose to set M i = M 0 at all interfaces, and optimize the interface placement. We proceeded as follows. First, a full FFS simulation was carried out with N 0 = 100 and n + 1 = 40 interfaces whose positions in phase space were chosen as described in Section III A. Note that once Φ A,0 in Eq. ( 4) has been determined, then if an estimate of the final rate k is available -for example from the Kramers' method [25,27,[START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF], or from previous runs with similar parameters -it is then possible to estimate a value of M 0 that yields a low relative variance of the rate V, as defined in Eq. (6). For instance, one may aim for V = 1.

For the present set of parameters, we used values of M 0 spanning from 1000 at high D, to 2800 at low D. Once the initial run has terminated, the interface placement is optimized to obtain p i = p opt i = (P B )

1 n ∀i, by following the procedure described in [34], and illustrated in Fig. 5. A monotonic interpolation function is needed in order to establish a one-to-one correspondence between the {p i } and the {ζ i }. A possible choice of such a function is,

f [λ(ζ i )] = i-1 j=0 ln p j n-1 j=0 ln p j , (7)
which reduces to the optimal values f (ζ opt i ) = i/n when p i = p opt i . The initial set of f (ζ init i ) is computed from the {ζ init i } according to Eq. ( 7). The sets of ζ init f (λ opt i ), as represented by pale blue dots. The corresponding values of the order parameter {ζ opt i } that determine the placement of the new interfaces are found by inverting the interpolation function, and are shown as pale blue triangles. Note that, for readability, we do not show all the interfaces in Fig. 5a. In practise, only a single FFS run was typically carried out per value of D, and the computed probabilities were then used to optimize the interfaces for the next value of D, and so on with iteratively decreasing D. as shown by the shaded area in Fig. 5b. For N 0 = 100, M 0 = 2800, the relative variance of the rate yielded by the FFS run is V = 1.37 and the computed lifetime of the uniform state at 300 K is τ = 0.013 ±0.002 s.

Lastly, Fig. 5c shows the evolution of the order parameter, ζ, and of the energy normalized by the thermal energy at T = 300 K, β 300 E, in which β 300 = (k B T 300 ) -1 is Boltzmann's factor at 300 K, as a function of the normalized reaction coordinate, χ, as computed by the geodesic nudged elastic band method (GNEB) [START_REF] Bessarab | [END_REF] in Ref. [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF]. χ goes from 0 in state 'A', to 0.5 at the barrier top, to 1 in state 'B'. The colored background indicates the number of interfaces per small interval of ζ with arbitrary units, with a darker color corresponding to more interfaces. During the initial FFS run, the interfaces were equally spaced out along the reaction coordinate, similarly to the images in the GNEB method, which are shown as blue triangles. The optimized interfaces are not evenly spaced out but are denser about halfway to the saddle point (ζ 0.58 and χ 0.27). In terms of the energy profile, this corresponds to the region preceding the barrier top. The fact that more interfaces are needed in that area is likely caused by the sudden variation in the flux of partial trajectories, similar to the ones shown in Fig. 3.

b. Results The Langevin simulations in the underdamped regime were performed with a homemade code [23,25] with an effective anisotropy and a single FFS run per value of D. Each run was carried out on a single CPU unit (Intel Xeon Gold 6130 processor), and a maximum of 65 simultaneous Langevin trial runs. A complete run at T = 300 K and n = 39 interfaces took between about a week at high D (M 0 1000), to about two weeks for low D (M 0 2800). In Fig. 6, we show, for a set of trajectories starting at λ 0 , the cumulative probability to reach λ i , P (λ i |λ 0 ), as a function of the order parameter ζ, for different values of D. We find that overall, the behavior of the system is the same as in the overdamped case of Section III B, with profiles of the cumulative fluxes similar to the ones shown in Fig. 4. Once more, the global flux of trajectories decreases with decreasing D. The inset in Fig. 6 shows the cumulative fluxes normalized by P B , i.e., P (λ i |λ 0 )/P B . We find that in this new viewpoint, all graphs are superimposed. Since the energy profile is symmetric, the rate of recrossing remains a constant fraction of the transition rate, and the DMI does not fundamentally impact the system in that sense.

D. Comparison of lifetimes

The key result of this study is presented in Fig. 7, where the mean lifetime of the 'A' state, τ , is shown as a function of D for the three cases discussed above. For memory elements, it corresponds to the information retention time. In Fig. 7(a), the results for the simulations with full dipolar interactions are shown. Since these results were obtained from only a single FFS run, we estimate the uncertainty in the lifetime by using the relative variance of the rate V, as defined in Eq. ( 6). We recall that M i = M 0 = 1000 for all interfaces for this case and N 0 = 50. The statistical error in the rate is given by [32],

σ k = k AB V N 0 , (8)
and is thus assumed to arise only from the trial runs. As anticipated from previous work, the lifetime appears to decrease with increasing DMI, which can be attributed to the linear decrease in the energy barrier with D that is directly related to the variation in the domain wall energy (per unit surface area), σ w = 4 √ AK u -πD. The striking result here however is the magnitude of τ , which varies from a few milliseconds at D = 0 to a few tenths of a microsecond at D = 2 mJ/m 2 . These values are in stark contrast with the τ predicted by taking f 0 = 1 GHz with the energy barriers computed elsewhere [5], as shown in the inset of Fig. 7(a). This discrepancy stems from a strong entropic contribution to the prefactor [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF], which results in a slower decrease of the overall lifetime as the barrier is reduced. In the present example, the constant f 0 approximation predicts a decrease by 10 orders of magnitude in τ as D is varied from 0 to 2 mJ/m 2 , whereas the FFS results indicate a change of only 4 orders of magnitude over the same range. As we have shown elsewhere [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF], compensation effects underpin this behavior and the present example serves to highlight the importance of quantifying the Arrhenius prefactor in magnetic nanostructures.

Fig. 7(b) presents the lifetimes calculated with the effective anisotropy approximation for the underdamped (α = 0.01) and overdamped (α = 0.5) limits. In the underdamped case, we observe that the overall lifetimes are two orders of magnitude greater than those shown in Fig. 7(a). This discrepancy is primarily due to the difference in energy barriers when the full dipolar interactions are taken into account, which amounts to a reduction of the barrier by about 4 k B T when dipolar interactions are included [5,[START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF]. As for the Arrhenius prefactor, we have found that, under a full treatment of the DDI, it is lower than in the effective aniostropy treatment for large D, while it is larger for low D [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF]. This leads to the apparent two regimes that can be seen in Fig. 7(a). As in Fig. 7(a), the uncertainties in the lifetime were determined from Eq. ( 8) for α = 0.5, whereas for α = 0.01, they were obtained from several distinct FFS runs. The latter method yields a larger error, because Eq. ( 8) neglects the statistical error that may stem from an incomplete sampling in the 'A' basin, and treats the trial runs as the only source of error. In the overdamped limit, the lifetimes decrease with the same overall slope as in the underdamped case, but are an order of magnitude shorter. This results from fewer barrier recrossings as compared with the underdamped case, where the domain wall exhibits more of a unidirectional motion as it sweeps through the disk. The system also relaxes faster, so the overdamped trajectories should, on average, remain closer to the MEP.

A comparison to direct Langevin simulations is possible for the lowest barriers studied, i.e., at D = 2 mJ/m 2 , where the overall lifetimes remain accessible to direct integration of the Eq. ( 2) for the full thermally-activated switching process. These data are also presented in Fig. 7, where good quantitative agreement with the FFS results can be seen. The statistics and averaged time traces from these direct Langevin simulations are pre- ] and is indicative of recrossing processes that ultimately lead to a longer average dwell times in the 'A' state. In contrast, the overdamped limit exhibits fewer excursions of this nature, where the time window in which large fluctuations occur are more localized to the switching event. This can also be seen in the solid black lines superimposed on the time series data, which represent averages over the 100 instances of the simulations and provide a measure of the reproducibility of the switching event. The averages for the underdamped cases exhibit a gradual decrease in m z toward switching, while the variation is significantly sharper in the overdamped

case.

In light of the previous observations, we show in black solid lines in Fig. 9, the cumulative fluxes of trajectories in FFS simulations normalized by P B , and averaged over all values of the DMI, as a function of the average order parameter, ζ. The transparent colored lines correspond to the individual normalized cumulative fluxes for each value of D, where the color of the line corresponds to the value of the DMI according to the colorscale in Fig. 6. Fig. 9(a) corresponds to the fluxes obtained with the effective anisotropy treatment of the dipole-dipole interactions, with two values of the Gilbert damping, respectively α = 0.01 (underdamped regime) and α = 0.5 (overdamped regime). As one can expect, the magnetization reverses more efficiently in the overdamped regime compared to the underdamped case, with an overall larger flux of forward trajectories along the order parameter and a faster saturation past the saddle point at ζ < 0, which indicates fewer barrier recrossings. This is in line with the behavior of the magnetization in direct Langevin simulations presented in Figs. 8(d) and 8(f). Fig. 9(b) shows the averaged cumulative fluxes of trajectories at a constant damping of α = 0.01, with the effective anisotropy, and the full treatment of dipole-dipole interactions. In the portion of configuration space situated before the saddle point (ζ > 0), the effect of DDI on the average fluxes resembles that of the large Gilbert damping in Fig. 9(a), with larger forward fluxes than in the effective anisotropy treatment. Nevertheless, the Langevin dynamics of the system under DDI, as shown in Fig. 8

IV. DISCUSSION AND CONCLUDING REMARKS

In this work, we have presented an example of how rate constants for magnetic systems may be obtained by means of forward flux sampling simulations. The case detailed here concerns the mean waiting times between magnetization reversals in a nanodisk with parameters similar to a free CoFeB layer, as used in magnetic random access memories, as a function of increasing DMI. In the large DMI case with the largest rate constant, FFS accurately reproduced the result of brute-force direct Langevin simulations, with and without a full treatment of dipole-dipole couplings. When full DDI are taken into account, and under a realistic Gilbert damping of α = 0.01, the lifetimes of the uniform states decrease with increasing DMI, and span from a few milliseconds to tenths of a microsecond, while the stability factor varies from ∆ = 32 to ∆ = 8 [5]. In all the cases studied above, we find that the assumption f 0 = 1 GHz for the Arrhenius prefactor [5,7] is not justified, and does not hold [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF].

In this particular class of ultrathin films, for which magnetization reversals take place via the nucleation and propagation of a domain wall, we have found that the average transition path obtained from forward flux sampling simulations closely reproduces the minimum energy path computed by path finding methods. In more complex systems, e.g., thicker films where DDI modify the transition paths, many paths from one stable state to the other may be statistically relevant [7,47]. Additionally, the path involving the lowest internal energy barrier may not correspond to the dominant transition mechanism, because the activation entropy of the transition accounts for a significant portion of the associated rate constant, and must therefore be taken into account [22,[START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the meyerneldel compensation rule[END_REF]. In such systems, FFS should prove particularly useful, as it directly samples the configuration space to compute rate constants, and the paths of highest statistical weight will naturally also be the ones most followed during the trial runs. Additionally, FFS can be used to obtain stationary distributions as a function of the order parameter, and thus observables allowing for comparison with experiments [33]. FFS is most efficient when the choice of the order parameter is able to closely mimic the reaction coordinate. In the present work, we chose the average z-component of the magnetization, ζ = m z , which relates to the reaction coordinate χ as ζ = tanh χ. While this choice of ζ has proved successful, it is possible that ζ = tanh m z may be an even better choice in this case.

FIG. 1 .

 1 FIG.1. Schematic of the forward flux sampling method. (a) Minimum energy path for magnetization reversal between the up and down states through domain wall nucleation and propagation. Progression along the reaction coordinate is parametrized by the order parameter ζ. (b) Sequence of interfaces {λ} separating the basins 'A' and 'B'. N0 crossings at the interface λ0 are recorded when the system starts in 'A' and successfully reaches λ0. (c) Simulated trajectories between subsequent interfaces. The micromagnetic configuration is stored at λi+1 for each successful crossing from λi to λi+1, which serves as a starting point for Langevin dynamics simulations toward the next interface. Failed crossings involve returns to the basin 'A'.

FIG. 2 .

 2 FIG. 2. Averaged magnetization configurations at different interfaces λ for different runs involving three values of the DMI, D. λ7 corresponds to the interface at which ζ = 0. Each configuration is rotated about the disk center such that wall propagation proceeds from left to right during reversal before the averaging procedure is performed.

FIG. 4 .

 4 FIG. 4. Cumulative probability, P (λi|λ0), as a function of the order parameter, ζ, for different values of D: (a) 0.5 mJ/m 2 , (b) 1.0 mJ/m 2 , (c) 1.5 mJ/m 2 , and (d) 2.0 mJ/m 2 . The insets show the conditional probability P (λi|λi-1) as a function of the interface λi. The colored lines represent single FFS runs, while the solid black lines represent averages over the different runs.

FIG. 5 .

 5 FIG. 5. Example of interface optimization in FFS. (a) The new values of the order parameter, {ζ opt i }, are obtained by inverting the interpolation function in Eq. (7) which is plotted in a dotted line. The {ζ opt i } are designed to approach a constant flux of partial trajectories through each interface. (b) Partial fluxes P (λi|λi-1) as a function of the order parameter ζ resulting from the interface optimization procedure. With the exception of the first and last interfaces, the values are found within the shaded area. (c) Evolution of the order parameter ζ (left), and of the Energy normalized by the thermal energy at 300 K, β300E (right), as a function of the reaction coordinate χ along the minimum energy path. The interfaces are optimized again after the first FFS run with the interfaces shown in (a). The number of interfaces per interval of ζ is represented by the colored background, where a darker color corresponds to more interfaces.

Fig. 5b shows

  FIG. 6. Cumulative probability, P (λi|λ0), as a function of the order parameter, ζ, for different values of D. The inset shows the same probabilities normalized by PB = P (λB|λ0). The normalized probabilities for different values of D are superimposed because the flux of barrier re-crossings remains a constant fraction of the rate constant, independently of D.

FIG. 7 .

 7 FIG. 7. Average lifetime, τ , as a function of the DMI, D, at T = 300 K computed with forward flux sampling (FFS) and direct Langevin dynamics simulations. (a) Lifetimes for α = 0.01 with full dipole-dipole interactions. (b) Lifetimes for α = 0.01 and α = 0.5 with effective perpendicular anisotropy. The average lifetime computed from full Langevin dynamics simulations is given for all three cases at D = 2 mJ/m 2 . The insets give a comparison of the FFS data with the constant f0 = 1 GHz approximation with the energy barriers respectively computed in Refs. 5 and 42.

1 FIG. 8 .

 18 FIG. 8. Direct Langevin dynamics simulations of the dwell time, τ , for D = 2 mJ/m 2 . The results correspond to (a,b) α = 0.01 with full dipolar interactions, (c,d) α = 0.01 with effective perpendicular anisotropy, and (e,f) α = 0.5 with effective perpendicular anisotropy. For each case, 100 simulations were performed with a successful switching event. (a,c,e) Probability distribution of the dwell time. The solid line represents a fit to an exponential function. (b,d,f) mz component of the magnetization as a function of time, where the curves are shifted along the time axis such that the first crossing of -mz,0 occurs at t = 0. The colored curves represent the individual Langevin dynamics simulations, while the solid black line represents an average of these curves.

FIG. 9 .

 9 FIG. 9. Effects of (a) damping and (b) dipole-dipole interactions (DDI) on the cumulative forward fluxes normalized by PB. The transparent colored lines show the cumulative fluxes for all values of D as a function of the order parameter, where the color of the line represents the value of D with the same colorscale as that of Fig. 6. In both graphs, the solid colored lines correspond to the underdamped case with effective anisotropy, while the dashed colored lines corresponds to the other case. The black lines show the average cumulative flux as a function of the average order parameter for each case.

  (b), is drastically different from that of Fig. 8(f). Since DDI result in large fluctuations in the order parameter in the region of the 'A' basin, with 0.8 m z -0.4, crossings of the inter-faces close to 'A', and up to the barrier top, are much more frequent. It follows that the portion of trajectories crossing subsequent interfaces in the ζ > 0 region, as opposed to the ones returning to 'A', is larger when DDI are included. This explains the larger fluxes in Fig. 9(b). Past the barrier (ζ < 0), the ratios of crossing trajectories over the ones returning to 'A' appear more similar with and without DDI, with similar amplitudes of the fluctuations during the reversal process in Figs. 8(b) and 8(d), and the forward fluxes become comparable.

ACKNOWLEDGEMENTS

The authors thank Nicolas Reyren for stimulating discussions and for sharing a dataset for the inset of Fig. 7. This work was supported by the Agence Nationale de la Recherche under Contract No. ANR-17-CE24-0025 (TOPSKY) and the University of Strasbourg Institute for Advanced Study (USIAS) for a Fellowship, within the French national programme "Investment for the Future" (IdEx-Unistra).