Switching probability of all-perpendicular spin valve nanopillars

In all-perpendicular spin valve nanopillars the probability density of the free-layer magnetization is independent of the azimuthal angle and its evolution equation simplifies considerably compared to the general, nonaxisymmetric geometry. Expansion of the time-dependent probability density to Legendre polynomials enables analytical integration of the evolution equation and yields a compact expression for the practically relevant switching probability. This approach is valid when the free layer behaves as a single-domain magnetic particle and it can be readily applied to fitting experimental data.

In all-perpendicular spin valve nanopillars the probability density of the free-layer magnetization is independent of the azimuthal angle and its evolution equation simplifies considerably compared to the general, nonaxisymmetric geometry. Expansion of the time-dependent probability density to Legendre polynomials enables analytical integration of the evolution equation and yields a compact expression for the practically relevant switching probability. This approach is valid when the free layer behaves as a single-domain magnetic particle and it can be readily applied to fitting experimental data.
With the recent advent of all-perpendicular Spin Transfer Torque Magnetoresistive Random-Access Memory (STT-MRAM) [1,2], a significant industrial and academic effort has materialized to better understand the physics of these devices and to improve their performance. STT-MRAM applications demand high thermal stability, switching speed, and energy efficiency. Assessing STT-MRAM experiments requires knowledge of the theoretically expected performance and there are a number of interrelated models that can be used for this purpose. The simplest and most common approach is to assume that the free layer is a single-domain magnetic particle and simulate the dynamics of the free-layer magnetization with the Landau-Lifshitz-Gilbert (LLG) equation [3], augmented by a term for the spin-polarized current [4,5]. In this formulation, stochastic fields are employed to introduce finite temperature into the LLG equation [6,7]. The resulting stochastic LLG equation must be simulated repeatedly for different realizations of the thermal fields to yield good statistics.
Alternatively, the advection-diffusion (Fokker-Planck) equation [6,7] for the probability density of the free-layer magnetization can be discretized and modelled using a standard numerical scheme. Only one copy of such a simulation is needed because the formalism does not involve any stochastic fields; it encompasses finite-temperature effects as a diffusion term in the Fokker-Planck equation. When high fidelity statistics are required, this approach can be faster than the stochastic LLG. Nonetheless, having to represent a distribution on a grid and then update it at each time-step is also computationally intensive.
Both simulation methods have been applied to study the physics of STT-MRAMs [8][9][10][11][12] but they are too cumbersome for fitting experimental data, especially when the simulation time is long. To circumvent this limitation, one usually identifies two regimes of distinct physical behaviors [6,13,14]: (a) a thermal regime, where the spin current is low and its effect can be approximated by modifying the energy barrier between the two stable states, and (b) a dynamic regime, where the current is high, the spin transfer torque dominates the thermal * michail.tzoufras@spintransfer.com fields, and the switching process is nearly deterministic. In this latter case the finite-temperature effects are only considered for initializing the system before the currentpulse arrives and are subsequently ignored during the dynamic switching process [10]. These approximations can be effective in predicting outlier events, principally in the limits of very low and very high current, but they do not cover the entire range of interesting parameters.
An analytical solution of the full thermal problem for the general geometry was derived in Refs. [15,16] by employing the spherical harmonic expansion. Thus, the derivatives in angle were reduced to algebraic expressions and yielded a hierarchy of differential-recurrence relations. These were then recast in matrix form and finally the solution was expressed in terms of Laplace transforms.
For all-perpendicular spin valves specifically, where the azimuthal angle drops out and the magnetization distribution ρ becomes a function of the polar angle θ and the time τ only [6], the spherical harmonic expansion reduces to a Legendre polynomial expansion. The Fokker-Planck equation may then be integrated for any current-pulse amplitude and duration, and the solution for ρ(θ, τ ) can be expressed in straightforward matrix notation. Obtaining the switching probability from ρ(θ, τ ) is a matter of evaluating the dot product between two vectors. Below, I start from the Fokker-Planck equation for the free-layer magnetization probability density ρ(θ, τ ) of an all-perpendicular spin valve as expressed in Ref. [17]: The four normalized quantities in Eq. (1) are: the energy barrier ∆ = µ0H k MsV , the external field h, the spin current i, and the time τ , where h, i, and τ may be denormalized following Ref. [17]: U is the bulk magnetocrystalline anisotropy, K S U /t is the surface anisotropy with t the film thickness, and N zz the demagnetizing factor. The parameters α and η correspond to Gilbert damping and spin polarization respectively.
The expansion of the distribution ρ(θ, τ ) to Legendre polynomials may be written as: Where the Legendre polynomials P n (x) are solutions to Legendre's differential equation d dx 1 − x 2 d dx P n + n (n + 1) P n = 0 and obey the recurrence relations: The zeroth order coefficient, r 0 , is proportional to the number density of the system as can be seen by integrating ρ on the spherical surface: 2π 0 π 0 ρ(θ, τ ) sin θdθdφ = 4πr 0 . Similarly, the first order term corresponds to the expected magnetization Substitution of the expansion (2) into the Fokker-Planck equation (1) and application of the recurrence relations (3)-(4) yields: with the coefficients a ı, defined as: a n,n = − n (n + 1) a n+2,n = (n + 1) (n + 2) (n + 3) (2n + 1) (2n + 3) Rearranging Eq. (5) and using the orthogonality of the Legendre polynomials enables integration of the Fokker-Planck equation: where r = (r n ) is the vector of the expansion coefficients and A = (a ı, ) is the pentadiagonal matrix with elements in Eqs. (6)- (10). All elements in the first row of A are identically equal to 0, i.e. a 0,n ≡ 0, ∀n. This ensures that-since there are no sources or sinks of spin valves in Eq.
In deriving Eq. (11), the spin current i and the external field h were assumed to be independent of time. In general, when ∂ ∂τ (i − h) = 0, the matrix A is timedependent with A(τ 1 )A(τ 2 ) = A(τ 2 )A(τ 1 ). The solution for r(τ ) should then be expressed using the Magnus expansion [18] instead of the matrix exponential from Eq. (11). Below, only the case of constant driver, With the probability density in terms of the expansion coefficients r(τ ) from Eq. (11), the probability that the magnetization is in the upper hemisphere, P(m z > 0), can be calculated as: where r s is the dot product between the vectors r and s = (s n ) with elements given by (see Ref. [19]): Typically, one is interested in the probability that a device with a given energy barrier, ∆, at some initial r(0), switches orientation for certain external field and current-pulse amplitude and duration. Using Eqs. Of the four operations listed above, only matrix exponentiation can be of non-trivial computational cost, and this only happens when the number of terms N max that are retained in the Legendre polynomial expansion is large.
However, for problems of practical interest N max ∼ 100 usually suffices and the resulting 100 × 100 matrix A can be immediately exponentiated by standard linear algebra software packages. The reason behind the rapid convergence of the expansion rests with Eqs. (6)-(10): while the off-diagonal terms scale as O(n) the diagonal ones scale as O(n 2 ), therefore ∂ ∂τ r n 1 − n(n+1) 2∆ r n 1 ⇒ r n 1 ∼ e − n(n+1) 2∆ τ 0. In other words, since ∆ is finite for T > 0K, highly peaked components of the distribution tend to diffuse away. In Fig. 1, I present a number of distributions generated by applying Eq. (11), where I vary N max . This figure confirms that a relatively modest increase in N max leads to a significant improvement in accuracy for the probability density ρ(θ, τ ).
In the limit n N max , the second term on the right hand side of Eq. (8), which derives from the anisotropy inherent in the magnetic material, becomes comparable with the first and counteracts diffusion, that is, it prevents r n>0 from vanishing. To see this, one can consider a system without a net external driver (i − h = 0) after a very long time (τ → ∞). In this scenario, it is known that the magnetization orientation is not completely randomized, but it lies along the anisotropy axis, and therefore r n>0 (τ → ∞) = 0. Still, for sufficiently large values of n the effect of magnetic anisotropy on r n is overwhelmed by thermal diffusion and r n diminishes as discussed above. The critical value n = n c for which the diffusion term overtakes the magnetic anisotropy and a n,n becomes negative is: n c = ∆/2 + 1−1/2. Hence, higher ∆ is associated with higher n c , and increasing the thermal stability produces magnetization distributions with sharper peaks on the anisotropy axis that require more terms in the Legendre polynomial expansion.
The off-diagonal terms a n±1,n are proportional to i − h. When these terms vanish, both directions on the anisotropy axis are equivalent; without a spin current or an external field there is nothing to break the symmetry of the bistable system. If one starts from a symmetric system without a net external driver and lets the system relax, the final distribution will also be symmetric.
In practice, one usually attempts repeated switches of a device, or an ensemble of devices, for a range of pulse amplitudes and durations and counts the number of failed switches for each set of parameters. This yields the error rate as a function of amplitude or pulse-length. The same approach can be taken using Eqs. (11)-(12): Starting from the trivial case r n>0 = 0, one generates a distribution as in Fig. 1 and then lets it relax by multiplying it with a suitable matrix exponential following Eq. (11). This relaxed distribution, r(0), becomes the initial condition from which new distributions are generated by applying Eq. (11). The error rate is calculated by substituting each new distribution, r(τ ), into Eq. (12).
For example, the "read disturbance" method for evaluating thermal stability via the parameter ∆ involves: attempting to switch a device with low current amplitude i, generating the error rate as a function of i, and fitting the measured data with an approximate expression such as: In this expression f 0 1GHz is the attempt frequency and the exponent β was calculated analytically in Ref. [10] for all-perpendicular geometry: β = 2. It is assumed that h = 0. The inner exponential derives from the Néel-Arrhenius model [20] and the outer from the approximation that the switching time is much longer than 1/f 0 , which also implies |i − h| 1. Unfortunately, it can be difficult to determine the experimental P |i−h| 1 (m z > 0) because very few switches happen under conditions of weak disturbance and a massive number of measurements needs to be taken to ensure adequate statistics. It is more convenient to obtain reliable data in the intermediate range, where |i−h| ∼ 1, which is beyond the range of validity of Eq. (14). In Fig. 2(a), the Read Error Rate [RER ≡ P(m z < 0)] is shown for pulse-length 100ns from two separate calculations: the solid blue line derives from Eqs. (11)-(12) and the dashed red line from Eq. (14). Good agreement between the two approaches is found in the limit |i − h| 1. However, Eq. (14) tends to overestimate the RER when |i − h| ∼ 1.
In Fig. 2 is shown as a function of the normalized spin current i = J/J c0 for various pulse-lengths. The parameter τ D = 2ns is chosen to facilitate comparison with Fig. 13 from Ref. [21], which was generated using Eq. (3.15) from Ref. [17]. Although the expressions in Refs. [17,21] recover the correct asymptotic behavior of the WER(J/J c0 ) in the limit J J c0 , they fail in the region J J c0 , especially for relatively long pulses, τ > 10. The lines shown in Fig.  2(b) were calculated by directly applying Eqs. (11)- (12).
Another useful plot is that of the switching current J WERt in terms of the pulse-length, where J WERt is defined as the spin current density required to achieve a specific WER target, P(m z > 0) = WER t . In Fig.  2(c), J WERt is shown as a function of the pulse-length for WER t = 0.5 as well as for smaller WER t values. These lines were calculated by numerically solving the transcendental equation 2π e Aτ r(0) s = WER t .
In summary, from the Fokker-Planck equation for the probability density of the free-layer magnetization of an all-perpendicular spin valve nanopillar, ρ(θ, τ ), a compact analytical solution, Eq. (11), was derived by employing the expansion of ρ(θ, τ ) to Legendre polynomials. From this expression the expected switching probability can be straightforwardly calculated using Eq. (12) under the assumption that the free layer acts as a single-domain magnetic particle.