Reversal time of jump-noise magnetization dynamics in nanomagnets via Monte Carlo simulations

Thermal fluctuation of magnetization in nanomagnets is important in the context of superparamagnetism, which is analyzed profoundly in Brown's seminal work,^1,2 and more recently reviewed by Coffey³ to estimate the reversal time of magnetization. Thermal fluctuation is more prominent in smaller magnetic volumes, which leads to spontaneous jumps from one stable state to another. From a practical perspective, this hinders the steady miniaturization of patterned spintronic device elements such as magnetic tunnel junction (MTJ) based magnetic random access memories (MRAMs).^4,5

The conventional way of modeling magnetization dynamics with thermal effects is by including two distinct and disjoint terms in the effective field of the macroscopic equation of motion:¹ (a) a dissipative field and (b) a random thermal field. For a uniformly magnetized particle with magnetization M of magnitude M_s, the classical equation of motion is given as^6,7

where γ₀ is the gyromagnetic constant, α is the damping constant, and H_T denotes an isotropic white Gaussian thermal field. The effective field H_eff includes the applied field and contributions from exchange, anisotropy, and magnetoelastic effects.⁸ The damping term, introduced by Gilbert,⁷ accounts for rapid relaxation toward the equilibrium state, while the thermal field is responsible for the Brownian motion.

However, there is an important assumption in modeling the dynamics classically: the time scale of the thermal field¹ and the relaxation of angular momentum associated with magnetic moment⁹ should be much shorter than the response times of the system. Simply put, the internal equilibrium of a magnetization orientation should occur much faster than the thermodynamic equilibrium of the ensemble.² The time scale associated with the thermal field is of the order of kT/h = 10⁻¹³ s at room temperature (kT is the thermal energy and h is Planck's constant), that of the angular momentum relaxation is of the order of 10^–15 s, and that of the precessional motion of magnetization is around 10⁻¹⁰ s.^1,9 Therefore, the assumption and the model are reasonable at nominal temperatures. However, the conventional model of magnetization relaxation becomes inadequate for explaining low temperature (10 m K–10 K) phenomena, such as macroscopic tunneling of magnetization,^10,11 in which the transition rate between two magnetic states is independent of temperature.

By modeling the dynamics due to thermal effects using a single random process composed of discontinuous transitions (jumps) with Poisson arrivals, called the jump-noise,^12,13 a unified model applicable over a broad range of temperatures is obtained. The jump-noise provides a justification for quantum tunneling of magnetization without invoking quantum mechanics.¹¹ For small noise, the classical Gilbert damping term emerges from average effects of the jump-noise.^13,14

The reversal time of magnetization extracted from jump-noise dynamics requires primarily the knowledge of the energy landscape of the magnetic body. As an application, it can be used to study the retention time and stability of emerging magnetic memories. An example is magnetoelectric antiferromagnetic (AFM) memory using (0001) Cr₂O₃ (chromia) thin film. In chromia thin films, the energy landscape strongly depends on the product of the applied longitudinal electric field and magnetic field.^15–17 However, experimental measurements of retention time of AFM domain state are elusive due to the absence of a net macroscopic magnetic moment in AFMs. The results obtained here can provide useful insight into the design and optimization of such memory systems.

The direct approach to realize jump-noise dynamics is to solve the Kolmogorov-Fokker-Planck equation for the transition probability density of the stochastic process using a specialized finite element and finite difference method. For small noise, this equation can be reduced to a stochastic process on energy graphs and can be set up and solved much faster.¹⁸ This paper focuses on the algorithmic approach to simulate jump-noise and extract the statistics empirically via Monte Carlo simulations. While Monte Carlo simulations are computationally expensive, they are fairly easy to implement and can be executed in parallel on a multi-core processor for significant speedup.^19,20

The goal of this paper is twofold. First, we present a complete methodology for numerical modeling of jump-noise statistics (Sec. III). Here, we contrast our work against prior works^21,22 in which simulations are performed using a sub-optimal global mesh. Second, we empirically extract the magnetization reversal time in nanomagnets (Sec. IV) with uniaxial and cubic anisotropy and compare the data with Néel-Brown theory's solution to analogous classical dynamics (Secs. V and VI). While Ref. 11 discusses the thermal switching rate, the reversal time of magnetization of jump-noise dynamics is analyzed in our work.

The magnetization dynamics driven by jump-noise is described by the equation¹² 

where $m = M / M s , h eff = H eff / M s , t = ( γ 0 M s ) t$ are dimensionless quantities. If the total magnetic free energy density g is known as a function of the state m, then $h eff = − ∇ m g ( m )$⁠. The jump-noise is represented by the weighted Dirac comb which expresses random jumps $Δ m i$ on the unit sphere $‖ m ‖ = 1$ occurring at random time instants t_i.

The jump-noise is formulated as a Poisson process characterized by a transition probability rate function S between state pairs $( m 1 , m 2 )$ in the phase space $‖ m ‖ = 1$⁠. The transition probability rate is given as¹³ 

where B and σ are the jump-noise parameters, μ₀ is the vacuum permeability, V is the volume of the nanomagnet, and $e b 0$ is the energy barrier parameter. The parameters B and σ can be determined experimentally by measuring the escape rate of the magnetization as a function of temperature. From Eq. (3), the scattering rate λ from a state m follows:

The probability density function f of a jump $Δ m i$ to occur at time t_i given the state $m i$ follows:

The statistic of the jump instants t_i is given by the inter-arrival times in a nonhomogeneous Poisson process

To simulate the jump-noise, we need to generate the correct statistics of the jumps $Δ m i$ and the instants $Δ t i$ from Eqs. (4)–(6) and evolve Eq. (2) in time.

The first step is to realize a discrete state space of $‖ m ‖ = 1$⁠. It is convenient to use spherical coordinates $( θ , ϕ )$ to represent the components of the state m such that $m x = sin θ cos ϕ , m y = sin θ sin ϕ , m z = cos θ$⁠. By uniformly sampling $θ ∈ [ 0 , π ]$ and $ϕ ∈ [ 0 , 2 π ]$⁠, we can form a finite element mesh on the unit sphere. In Ref. 22, the spherical mesh forms a global state space, which remains fixed regardless of the state from where jump occurs. Discretizing the space this way is not ideal, especially when $σ 2 ≪ 1$⁠, because majority of the jumps would occur in a small solid angle centered at the current state. Moreover, the error due to discretization depends on the state as the mesh points are unevenly distributed [see Fig. 1(a)].

A better approach, adopted in this work, is to form a local spherical mesh with half angle Θ as shown in Fig. 1(b). Such mesh inherently keeps the density of states high near the local state and less far away, which is beneficial because the probability density (5) behaves alike. Tessellation of the sphere by a regular convex polyhedron²³ does not exhibit this property. Θ can be estimated by first finding an upper bound on the transition probability rate (3) for $‖ Δ m ‖ = 2 sin ( Θ / 2 )$ as

then setting the right-hand-side to logarithm of the desired suppression. Θ measures the strength of jump-noise given the parameters σ and $e b 0$⁠.

The random samples from the probability density function in Eq. (5) are generated by a fundamental method called inverse transform sampling. First, label the state space by a single index j such that $S = { ( θ j , ϕ j ) : j = 1 , 2 , … , M }$ is the collection of states, where M is the total number of states. Now, find the cumulative distribution function $F$ as

where δ denotes the sample spacing. Generate a random number v from the uniform distribution between 0 and 1, and compute the value of j satisfying

The state $m j ≐ ( θ j , ϕ j )$ is the random jump destination.

The jump statistic given by Eq. (6) cannot be numerically realized in its current implicit form. This is resolved by homogenizing the scattering on the state space via a self-scattering (nil jump) process λ₀.²¹ The total scattering rate after homogenization is given by

From Eqs. (6) and (10), the time between scattering τ can be explicitly expressed as

By replacing the probability with a uniform random number between 0 and 1, the random duration between the jumps is generated. But, the probability of a non-zero jump is only $λ ( m ) / Γ$⁠, which can be discerned by generating another uniform random number u between 0 and 1, and checking if $u ≤ λ ( m ) / Γ$⁠.

Lastly, the state m is evolved in time according to Eq. (2) by repeating the following steps: (a) integrate the deterministic term in the scattering-free duration $t ∈ [ t i , t i + 1 )$ using a finite difference method, (b) at $t = t i + 1$ perform the jump to $m j$ and generate a new duration τ such that the next jump occurs at $t i + 2 = t i + 1 + τ$⁠.

To corroborate the numerical method to emulate jump-noise dynamics, we perform Monte Carlo simulations to study the equilibrium distribution of magnetization for nanomagnets possessing magnetocrystalline anisotropy only, under no applied field. We consider 1000 samples aligned along the same lowest energy state $m z = − 1$⁠, without loss of generality, at time t = 0, and let the ensemble evolve with time until the absolute ensemble mean $| m z ¯ | < = 0.001$⁠. The simulation is implemented in MATLAB with the help of Parallel Computing Toolbox on a server with 20-core CPU @ 2.3 GHz and 512 GB memory. The time step is 0.05, the angular spacing of mesh is $0.72 °$⁠, and the suppression factor for Θ is $10 − 8$⁠. The computation time for these specifications scales $∼ 0.1 × τ$ s, where τ is the reversal time (16).

Analytically from statistical mechanics, the probability density function of the states m in thermal equilibrium is given by the Boltzmann distribution

where Z is the normalization constant given by

From Eq. (12), the equilibrium distribution of m_z follows:

For uniaxial anisotropy, $g ( θ , ϕ ) = 1 2 sin 2 θ$⁠. For cubic anisotropy, $g ( θ , ϕ ) = 1 2 ( sin 4 θ sin 2 2 ϕ + sin 2 2 θ )$⁠.

The simulation results for the equilibrium distribution of m_z for the uniaxial and cubic anisotropy are shown in Fig. 2. The Monte Carlo simulations replicate the Boltzmann distribution correctly for both the energy landscapes. This is expected because the transition probability rate in Eq. (3) satisfies the detailed balance

for all $( m 1 , m 2 ) ∈ ‖ m ‖ = 1$⁠. The small yet notable error in the simulations is affected by one or more of the sample size and mesh spacing.

In Néel-Brown theory,¹ the reversal time of magnetization τ characterizes the longitudinal relaxation of ensemble magnetization as

Even for jump-noise dynamics, Eq. (16) holds true because the Markov chain represented by it satisfies the detailed balance condition, and therefore has an exponential rate of convergence to the stationary distribution (12).²⁴ So, τ can be estimated from the asymptotic value of $τ ( t ) = − t / ln [ | m z ¯ | ]$ from simulations as featured in Fig. 3, by following the same procedure mentioned for extracting the equilibrium distribution.

It is essential to see how the reversal time of magnetization extracted from jump-noise dynamics compares with the reversal time obtained from Néel-Brown theory. This requires establishing an equivalence between the jump-noise and classical magnetization dynamics.

Equation (1) can be rewritten in the dimensionless Landau-Lifshitz form,² excluding the thermal field, as

In the absence of applied field, $h eff = − ∇ g ( θ , ϕ )$ is orthogonal to the radial direction m, so

From Eqs. (17) and (18), the component of $d m / d t$ along $∇ g$⁠, in the absence of applied field, gives the damping term $α / ( 1 + α 2 )$⁠. Indeed for the jump-noise dynamics (2), the average component of jump-noise along $∇ g$ gives rise to an equivalent damping term locally.

The expected value of jump-noise at a given state m is expressed as

For small noise $Θ ≪ 1$⁠, the phase space of jumps $Δ m ′$ can be approximated by the local tangent plane. The energy difference $g ( m + Δ m ′ ) − g ( m )$ in Eq. (3) can be truncated up to second order of Taylor series, which then can be used to evaluate the integral in Eq. (19) definitely by forming a bivariate Gaussian function. The derivation is similar to that for first order approximation of the energy difference in Ref. 14. The inclusion of second-order term is important because the term $‖ Δ m ′ ‖ 2 / σ 2$ in Eq. (3) is quadratic. The final expression of the expected jump-noise is obtained as

where $I 2$ is the 2 × 2 identity matrix and $H$ is the Hessian operator. In local coordinates of the tangent space of $m ( θ , ϕ )$⁠, the gradient and Hessian of g are written as

where each letter of subscript denotes partial derivative with respect to that variable.

The equivalent damping term for jump-noise can now be determined as

Unless $W$ equals identity, the expected value of jump-noise would also yield an orthogonal component along $m × ∇ g = [ − csc θ g ϕ g θ ] T$⁠, which is the precessional motion. The physical significance of this component is a small correction to the gyromagnetic ratio [see Eq. (1)] by a factor of

In Néel-Brown theory, the reversal time in the high energy barrier limit $e b 0 ≫ 1$ is approximated by inverse of the smallest non-vanishing longitudinal eigenmode of the Fokker-Planck equation of the Landau-Lifshitz dynamics.^2,3 The expressions of reversal time for the uniaxial and cubic anisotropy, τ_uni and τ_cub, respectively, in dimensionless units are reproduced from Ref. 3 as

The expression of τ_cub here is multiplied by 2 because the reversal time for cubic anisotropy in the literature^2,3,25 corresponds to the relaxation associated with surmounting one energy barrier from m_z = 1 to $m x , m y = ± 1$ (see Fig. 4). Since the energy landscape is symmetric about m_z = 0 and the barriers are high, the longitudinal relaxation from m_z = 1 to $m z = − 1$ is associated with overcoming two successive energy barriers with intermediate orientations at $m x , m y = ± 1$⁠.

The equivalent damping term α and the gyromagnetic correction factor $γ c$ for jump-noise dynamics are not constants but depend on the state $m ( θ , ϕ )$⁠. To obtain the corresponding Néel-Brown theory's solution of the reversal time, we calculate the numerical bounds of τ_uni (27) and τ_cub (28) by varying α (25) and $γ c$ (26) on the phase space. Upon inspection, when the noise is small, the lower bound occurs close to energy maximum (⁠ $θ = π / 2$⁠) for τ_uni and close to energy saddle points (⁠ $θ = π / 4 , 3 π / 4$⁠; $ϕ = 0 , π / 2 , π , 3 π / 2$⁠) for τ_cub. For large energy barrier, the internal equilibrium within an energy well occurs much faster than the equilibrium between energy wells, leaving the “flow” across energy barrier to be concentrated near the point of least height,² thereby giving the lower bound. The upper bound occurs close to energy minima for both τ_uni and τ_cub because the damping is minimum around it.

The reversal time extracted from Monte Carlo simulations of jump-noise dynamics and the numerical bounds of the reversal time obtained from Néel-Brown theory's solution to equivalent Landau-Lifshitz dynamics are shown in Fig. 5. For both uniaxial and cubic anisotropy, the simulation points coincide with the lower bound for smaller values of σ and deviate from it for larger values while remaining within the bounds. Thus, the reversal time of jump-noise dynamics for small noise is characterized by saddle point values of the equivalent classical dynamics parameters α and $γ c$⁠. The reversal time of jump of dynamics is also always smaller than that of equivalent classical dynamics characterized by energy minima values of α and $γ c$⁠. The qualitative explanation for this is that jump-noise phenomenologically accounts for tunneling unlike classical dynamics, and so yields a lower reversal time.

For large noise, the second order approximation of the energy difference breaks down, and consequently the expressions of α and $γ c$ could yield unphysical results. Additionally, the expected value of jump-noise would produce component along the radial direction which cannot be modeled by the Landau-Lifshitz equation. For large scattering rate, the right hand side of Eq. (25) could exceed unity, thereby giving a complex damping parameter.

The magnitude of noise for desired suppression of probability is determined from the noise strength Θ using Eq. (7). Alternatively, from Eqs. (5) and (20), the value of $σ 2 e b 0 ‖ W − 1 ∇ g ‖$ gives a good measure of noise. However, when comparing different energy landscapes, the magnitude of noise must be measured relative to the distance between the energy minima and maxima. When the noise is large enough to cross the barrier in few jumps rather than gradually, the equilibrium occurs as a simultaneous process within and in between energy wells; therefore, Néel-Brown theory is not valid. Since the angular distance between the energy minima and maxima is larger in uniaxial (⁠ $90 °$⁠) than in cubic (⁠ $55 °$⁠) anisotropy, for a $σ 2 e b 0$ product say, $0.25 2 × 10 = 36 °$⁠, evidently the reversal time of jump-noise dynamics for cubic anisotropy shows larger deviation from Néel-Brown theory's solution in Fig. 5.

In this paper, we develop numerical methods to extract the reversal time of magnetization of jump-noise dynamics in nanomagnets with uniaxial and cubic anisotropy via Monte Carlo simulations. The reversal time of nanomagnets due to thermal fluctuation predicted from jump-noise dynamics is more accurate because it accounts for tunneling effects by its very formulation. Modeling of the jump-noise requires primarily the knowledge of the energy landscape of the system and not the microscopic origin of the underlying phenomena. Results show that the reversal time gathered from Néel-Brown theory's solution to analogous classical dynamics (a) at the energy saddle point characterizes the reversal time of jump-noise dynamics for small noise, (b) at the energy minima is always larger than the reversal time of jump-noise dynamics. For large noise, the magnetization reversal due to jump-noise dynamics has no classical analogue and phenomenologically represents macroscopic tunneling of magnetization.

This work was supported in part by the Semiconductor Research Corporation (SRC) and the National Science Foundation (NSF) through ECCS 1740136. S. Rakheja also acknowledges the funding support from the MRSEC Program of the National Science Foundation under Award No. DMR-1420073.